For the purposes of this chapter, I will consider vector data as consisting of a finite set of points that may or may not be connected. In R, these points can be used directly, as a numeric data type, or as a spatial simple feature object. The sf (Simple Features) library contains a lot of functions that help deal with points, lines, and polygons, and this is going to be a short overview of how they can be derived and manipulated in the pursuit of population genetic studies.
The older S4-based sp package required a rather convoluted set of nested constructions to build spatial objects. The modern standard in the R spatial ecosystem is the sf package, which integrates seamlessly with tidyverse and standard data.frame layouts, making spatial objects much easier to construct, view, and manipulate.
Points are defined by simple geometries wrapped into a collection of coordinates. In sf, we represent spatial geometries as standard data.frame or tibble objects with an additional list-column called geometry that contains the coordinate information.
Let’s start by extracting the centroids of the Arapatus attenuatus populations from our data set:
library(gstudio)
data(arapat)
coords <- strata_coordinates(arapat)
summary(coords) Stratum Longitude Latitude
101 : 1 Min. :-114.3 Min. :23.08
102 : 1 1st Qu.:-112.9 1st Qu.:24.52
12 : 1 Median :-111.5 Median :26.21
153 : 1 Mean :-111.7 Mean :26.14
156 : 1 3rd Qu.:-110.4 3rd Qu.:27.47
157 : 1 Max. :-109.1 Max. :29.33
(Other):33 We can easily turn this standard data.frame of coordinates into a spatial simple features object (sf) using the st_as_sf() function. We must specify which columns contain the coordinates and define the coordinate reference system (CRS). In this case, we use the standard WGS84 geographic reference system, which has the EPSG code 4326.
Simple feature collection with 39 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -114.2935 ymin: 23.0757 xmax: -109.1263 ymax: 29.32541
Geodetic CRS: WGS 84
First 10 features:
Stratum geometry
1 88 POINT (-114.2935 29.32541)
11 9 POINT (-113.9449 29.01457)
20 84 POINT (-113.6679 28.96651)
29 175 POINT (-113.4897 28.72796)
36 177 POINT (-113.9914 28.66056)
46 173 POINT (-112.8698 28.40846)
56 171 POINT (-113.1826 28.22308)
66 89 POINT (-113.3999 28.03661)
76 159 POINT (-113.3161 27.52944)
85 SFr POINT (-112.964 27.3632)Notice that the output shows that it is a standard data frame but with additional spatial attributes: the geometry type (POINT), the bounding box, and the coordinate reference system (CRS: wgs84).
If we ever need to fetch or redefine the CRS of our spatial object, we can use the st_crs() function:
st_crs(pts)Coordinate Reference System:
User input: EPSG:4326
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]Any set of x- and y- coordinates can be turned into a spatial sf object. Since sf objects are standard data frames, associating data with those points is as simple as merging, binding, or adding columns directly to the data frame.
Let’s determine the number of individuals genotyped in each population and add this data to our points:
df <- data.frame( table(arapat$Population) )
names(df) <- c("Population","N")
# We can merge this data directly with our sf object (matching "Stratum" to "Population")
pts.df <- merge(pts, df, by.x="Stratum", by.y="Population")
pts.dfSimple feature collection with 39 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -114.2935 ymin: 23.0757 xmax: -109.1263 ymax: 29.32541
Geodetic CRS: WGS 84
First 10 features:
Stratum N geometry
1 101 9 POINT (-110.5744 27.90509)
2 102 8 POINT (-109.1263 26.38014)
3 12 10 POINT (-112.6655 27.18232)
4 153 10 POINT (-110.4624 24.13389)
5 156 6 POINT (-109.989 24.0438)
6 157 10 POINT (-110.096 24.0195)
7 159 9 POINT (-113.3161 27.52944)
8 160 10 POINT (-112.5296 27.40498)
9 161 10 POINT (-112.986 27.0367)
10 162 10 POINT (-112.408 27.2028)You can manipulate this spatial object exactly like any other data frame, or you can drop the spatial geometry entirely if you just want a plain data frame:
st_drop_geometry(pts.df)[1:5,] Stratum N
1 101 9
2 102 8
3 12 10
4 153 10
5 156 6Since it is a simple feature object, you can easily fetch spatial properties like the bounding box (the coordinates of a box that encloses all the points) using st_bbox():
st_bbox(pts.df) xmin ymin xmax ymax
-114.29353 23.07570 -109.12633 29.32541 Lines are created by linking sequences of coordinate points. In sf, we define a line geometry using st_linestring() on a matrix of coordinate points, and wrap these geometries in a geometry list-column using st_sfc().
# Grab coordinates for first three populations
c1 <- matrix(c(coords$Longitude[1:2], coords$Latitude[1:2]), ncol=2)
c2 <- matrix(c(coords$Longitude[2:3], coords$Latitude[2:3]), ncol=2)
L1 <- st_linestring(c1)
L2 <- st_linestring(c2)
# Create a spatial geometry column with a defined CRS
SLs <- st_sfc(list(L1, L2), crs = 4326)
SLsGeometry set for 2 features
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: -114.2935 ymin: 28.96651 xmax: -113.6679 ymax: 29.32541
Geodetic CRS: WGS 84If we want to add data to the set of lines, we can combine our geometry column with a data.frame of attributes using the st_sf() function:
df <- data.frame( Sequence = c("First","Second"), Like_It = c(TRUE, FALSE) )
SLDF <- st_sf(df, geometry = SLs)
SLDFSimple feature collection with 2 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: -114.2935 ymin: 28.96651 xmax: -113.6679 ymax: 29.32541
Geodetic CRS: WGS 84
Sequence Like_It geometry
1 First TRUE LINESTRING (-114.2935 29.32...
2 Second FALSE LINESTRING (-113.9449 29.01...We can also extract the length of each line in meters (or other geographic units) using the st_length() function:
st_length(SLs)Units: [m]
[1] 48380.9 27465.6A polygon is simply a collection of line segments that closes in on itself. We can use polygons to identify habitat, define boundaries, etc.
To make a polygon, we must close the coordinates, which means taking the first coordinate and appending it to the end of the coordinate matrix, such that the first row is identical to the last row.
c <- matrix(c(coords$Longitude[1:5], coords$Latitude[1:5]), ncol=2)
c <- rbind(c, c[1,])
P <- st_polygon(list(c))
PSimilar to lines, we can bundle this geometry into a spatial list-column and define the CRS:
Geometry set for 1 feature
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -114.2935 ymin: 28.66056 xmax: -113.4897 ymax: 29.32541
Geodetic CRS: WGS 84And attributes/data can be added to it easily using st_sf():
df <- data.frame(Populations=paste(coords$Population[1:5], collapse=", "))
SPDF <- st_sf(df, geometry = SPs)
SPDFSimple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -114.2935 ymin: 28.66056 xmax: -113.4897 ymax: 29.32541
Geodetic CRS: WGS 84
Populations geometry
1 POLYGON ((-114.2935 29.3254...As all of these objects are R objects, they can be saved to disk using the standard save() function, which makes them a .rda object, or we can export them directly to standard spatial file formats (like GeoJSON, Shapefiles, or GeoPackage) using the st_write() function from the sf package.
st_write(pts.df, "populations.geojson")