3.3 Vector objects in R

In this section we present two packges: sp and sf. The sp package actualy defines classes for both vector and raster objects. Below, we however focus on the vector ones and so we do not detail SpatialGrid and SpatialPixels objects. Also note that sf “[…] aims at succeeding sp in the long term” (Simple Features for R, sf vignette).

Let first create a data frame mydata:

mylon <- -82 + 2 * runif(20)
mylat <- 42 + 2 * runif(20)
mydata <- data.frame(
  lon = mylon,
  lat = mylat,
  var1 = rnorm(20),
  var2 = 10 * runif(20)
)

Let’s have a look at thus data frame:

lon	lat	var1	var2
-81.03876	42.94418	0.1110215	5.591651
-81.39501	42.73904	0.6644820	8.010553
-80.46684	43.70218	2.2764582	1.790977
-81.52099	43.49017	0.5213265	5.067142
-81.29007	42.58602	2.2475670	7.690447
-81.57360	42.57142	0.8358582	2.518603

3.3.1 Package `sp`

3.3.1.1 Classes

The table below includes a description of the classes for points, lines and polygons. Basically all these classes work the same way. For instance, in order to define a SpatialPointsDataPoints object, three elements are required: a set of coordinates, a Coordinate Reference System (CRS) and an attribute table. Intermediate class are also defined (for instance points + CRS = SpatialPoints) and the name of the class is also the name of the function to be called.

Classes and functions	Contents
Points	list of points (set of coordinates)
SpatialPoints	list of points + CRS
SpatialPointsDataPoints	list of points + CRS + attribute table
Line	a line (set of coordinates)
Lines	list of lines
SpatialLines	list of lines + CRS
SpatialLinesDataFrame	list of lines + CRS + attribute table
Polygon	a line (set of coordinates)
Polygons	list of lines
SpatialPolygons	list of lines + CRS
SpatialPolygonsDataFrame	list of lines + CRS + attribute table

3.3.1.2 `SpatialPointsDataFrame`

As an more tangible example, let’s now create a SpatialPointsDataFrame based on our data frame mydata.

library(sp)
mysp <- SpatialPointsDataFrame(
  coords = mydata[, 1:2],
  data = mydata[, 3:4],
  proj4string = CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
  )
)

Note that proj4 is used and we therefore wrote a string that describes the CRS and that proj4 understands. Below are listed some properties of the object we have defined.

isS4(mysp)
class(mysp)
slotNames(mysp)
dim(mysp)

#R>  [1] TRUE
#R>  [1] "SpatialPointsDataFrame"
#R>  attr(,"package")
#R>  [1] "sp"
#R>  [1] "data"        "coords.nrs"  "coords"      "bbox"        "proj4string"
#R>  [1] 20  2

Basically, it is a S4 object of class SpatialPointsDataFrame. All slot names refer to attribute that are accessible via and @:

mysp@proj4string

#R>  Coordinate Reference System:
#R>  Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
#R>  WKT2 2019 representation:
#R>  GEOGCRS["unknown",
#R>      DATUM["World Geodetic System 1984",
#R>          ELLIPSOID["WGS 84",6378137,298.257223563,
#R>              LENGTHUNIT["metre",1]],
#R>          ID["EPSG",6326]],
#R>      PRIMEM["Greenwich",0,
#R>          ANGLEUNIT["degree",0.0174532925199433],
#R>          ID["EPSG",8901]],
#R>      CS[ellipsoidal,2],
#R>          AXIS["longitude",east,
#R>              ORDER[1],
#R>              ANGLEUNIT["degree",0.0174532925199433,
#R>                  ID["EPSG",9122]]],
#R>          AXIS["latitude",north,
#R>              ORDER[2],
#R>              ANGLEUNIT["degree",0.0174532925199433,
#R>                  ID["EPSG",9122]]]]

head(mysp@data)

#R>         var1     var2
#R>  1 0.1110215 5.591650
#R>  2 0.6644820 8.010553
#R>  3 2.2764582 1.790977
#R>  4 0.5213265 5.067142
#R>  5 2.2475670 7.690447
#R>  6 0.8358582 2.518603

In order to change projection, the user must call spTransform(), like so:

(mysp2 <- spTransform(mysp, CRS=CRS("+proj=merc +ellps=GRS80")))

#R>  class       : SpatialPointsDataFrame 
#R>  features    : 20 
#R>  extent      : -9118633, -8957527, 5159576, 5389937  (xmin, xmax, ymin, ymax)
#R>  crs         : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
#R>  variables   : 2
#R>  names       :              var1,              var2 
#R>  min values  : -2.12149445648486, 0.957436570897698 
#R>  max values  :   2.2764581628672,  9.31604944868013

3.3.2 Package `sf`

Below is a very short overview of classes in sf, the reader that requires further explanation would find more details on the first vignette of sf. Basically three classes are defined: sf, sfc and sfg.

3.3.2.1 Class `sf`

library(sf)

#R>  Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE

pts_sf <- st_as_sf(
  x = mydata,
  coords = c("lon", "lat"),
  crs = 4326
)

Let’s examine its class

class(pts_sf)

#R>  [1] "sf"         "data.frame"

3.3.2.2 Class `sfc`

pts_sfc <- st_geometry(pts_sf)
class(pts_sfc)

#R>  [1] "sfc_POINT" "sfc"

3.3.2.3 Class `sfg`

(x <- st_point(c(1, 2)))

#R>  POINT (1 2)

class(x)

#R>  [1] "XY"    "POINT" "sfg"

3.3.2.4 How to import a `sp` object

st_as_sf() can also be used to convert a sp object into a sf one.

st_as_sf(mysp)

#R>  Simple feature collection with 20 features and 2 fields
#R>  Geometry type: POINT
#R>  Dimension:     XY
#R>  Bounding box:  xmin: -81.91407 ymin: 42.18196 xmax: -80.46684 ymax: 43.70218
#R>  CRS:           +proj=longlat +datum=WGS84 +no_defs
#R>  First 10 features:
#R>           var1     var2                   geometry
#R>  1   0.1110215 5.591650 POINT (-81.03876 42.94418)
#R>  2   0.6644820 8.010553 POINT (-81.39501 42.73904)
#R>  3   2.2764582 1.790977 POINT (-80.46684 43.70218)
#R>  4   0.5213265 5.067142 POINT (-81.52099 43.49017)
#R>  5   2.2475670 7.690447 POINT (-81.29007 42.58602)
#R>  6   0.8358582 2.518603  POINT (-81.5736 42.57142)
#R>  7   0.4324725 4.854437 POINT (-81.85505 43.30567)
#R>  8  -1.0001346 5.265450 POINT (-81.32604 42.18196)
#R>  9  -0.7722097 4.248399 POINT (-81.89046 42.54436)
#R>  10  0.6212614 7.869162 POINT (-81.63948 42.46887)