library(sf)#"s2" is an enhancement of the "sf" package#It allows one to complete spatial operations with data that is in geographic coordinate system (the coordinates are in degrees)#However, having it enables will take much longer computation timesf_use_s2(FALSE)library(ggplot2)#Step1: Read country shapeusa_cntry <-st_read(dsn="./data/gadm41_USA_shp/gadm41_USA_0.shp", quiet =TRUE)#Step2: Simplify linesusa_cntry<-st_simplify(usa_cntry, dTolerance =0.05)#Define zoommin_lon_x_us<-(-130)max_lon_x_us<-(-57)min_lat_y_us<-(26)max_lat_y_us<-(53)#Mapfig<-ggplot()+geom_sf(data=usa_cntry, linewidth =0.3, fill =NA)+theme_bw()+coord_sf(xlim =c(min_lon_x_us-3, max_lon_x_us+3), ylim =c(min_lat_y_us-3, max_lat_y_us+3))
Mapping the US USA_0
fig
Mapping the US USA_0
Let us examine what the dataframe looks like
head(usa_cntry, 2)
Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -178.9911 ymin: 18.91028 xmax: 179.7485 ymax: 71.33542
Geodetic CRS: WGS 84
GID_0 COUNTRY geometry
1 USA United States MULTIPOLYGON (((179.6461 52...
Mapping the US USA_0
This is how we turn the sf object into a simple dataframe
Geometry objects can be collected into a geometry column - sfc, which is done with st_sfc
A geometry column also contains a Coordinate Reference System
Geometry column - sfc
Four types of CRS definitions are acceptable
An EPSG code (e.g. 4326)
A crs object, as returned by the st_crs function
A PROJ4 string (e.g. '+proj=longlat +datum=WGS84 +no_defs')
A WKT string
If the crs argument is omitted, the default is to set an “undefined” CRS (i.e., NA).
EPSG Codes
When you need precise measurements for things like distances, you should use EPSG codes.
EPSG Geodetic Parameter Dataset is a public registry of geodetic datums, spatial reference systems, Earth ellipsoids, coordinate transformations and related units of measurement.
It came from a member of the European Petroleum Survey Group in 1985.
Different countries have different codes.
EPSG Codes
For example, the standard for the whole world is 4326, the equivalent of WGS84 (also used by Google)
We can combine the attribute table name_df and the geometry column geom (sfc) using st_sf, resulting in a layer
layer <-st_sf(name_df, combined)layer
Simple feature collection with 3 features and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 34.79844 ymin: 31.06862 xmax: 35.01164 ymax: 31.26028
Geodetic CRS: WGS 84
name combined
1 Point 1 POINT (34.79844 31.24329)
2 Point 2 POINT (34.81283 31.26028)
3 Point 3 POINT (35.01163 31.06862)
Mapping the sf layer with plot
plot(layer)
Mapping the sf layer with ggplot
ggplot()+geom_sf(data=layer)+theme_bw()
Mapping the sf layer with ggplot
We can also plot the polygons we made with ggplot
Here is what ab looks like:
ggplot()+geom_sf(data=ab)+theme_bw()
Extracting Layer Components
Sometimes, we are interested in going the other way around and extract:
the geometry
the attribute table
the coordinates
Extracting Layer Components: the Geometry
The following command allows us to extract the geometry
st_geometry(layer)
Geometry set for 3 features
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 34.79844 ymin: 31.06862 xmax: 35.01164 ymax: 31.26028
Geodetic CRS: WGS 84
POINT (34.79844 31.24329)
POINT (34.81283 31.26028)
POINT (35.01163 31.06862)
Extracting Layer Components: the Attribute Table
The following command allows us to drop the geometry
We are then left just with the attribute table:
st_drop_geometry(layer)
name
1 Point 1
2 Point 2
3 Point 3
Extracting Layer Components: the Coordinates
The coordinates of sf, sfc or sfg objects can be obtained with the function st_coordinates.
The coordinates are a matrix:
st_coordinates(layer)
X Y
[1,] 34.79844 31.24329
[2,] 34.81283 31.26028
[3,] 35.01163 31.06862
We can turn the matrix into a dataframe in the following way
data.frame(st_coordinates(layer))
X Y
1 34.79844 31.24329
2 34.81283 31.26028
3 35.01163 31.06862
name addr.street lat lon
1 Pizzeria ai Marmi Viale di Trastevere 41.88826 12.47379
2 Sichuan Haozi Via di San Martino ai Monti 41.89580 12.49948
3 Dar filettaro a Santa Barbara Largo dei Librari 41.89467 12.47370
The restaurants dataframe can be converted to a sf layer using st_as_sf:
Simple feature collection with 2 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 12.47379 ymin: 41.88826 xmax: 12.49948 ymax: 41.8958
Geodetic CRS: WGS 84
name addr.street geometry
1 Pizzeria ai Marmi Viale di Trastevere POINT (12.47379 41.88826)
2 Sichuan Haozi Via di San Martino ai Monti POINT (12.49948 41.8958)
Creating a Point Layer from a table
This is what these points look like, if we map them out:
ggplot()+geom_sf(data=df2, size=0.3)
Creating a Point Layer from a table
Without some context, it is difficult to see what those points
You can see below an updated map with the streets in the background.
sf properties: rows and columns
Some important properties of sf objects include rows and columns.
Number of rows:
nrow(df2)
[1] 2811
Number of columns:
ncol(df2)
[1] 3
We can obtain both rows and column by typing:
dim(df2)
[1] 2811 3
sf properties: bounding box coordinates
The st_bbox function returns the bounding box coordinates
Subsetting columns is very similar to subsetting a dataframe.
The difference is the geometry column that stays.
We turn the spatial dataframe into a regular dataframe.
df3<-st_drop_geometry(df2)
Subsetting based on attributes
This is how the two compare:
head(df2, n=3)
Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 12.4737 ymin: 41.88826 xmax: 12.49948 ymax: 41.8958
Geodetic CRS: WGS 84
name addr.street
1 Pizzeria ai Marmi Viale di Trastevere
2 Sichuan Haozi Via di San Martino ai Monti
3 Dar filettaro a Santa Barbara Largo dei Librari
geometry
1 POINT (12.47379 41.88826)
2 POINT (12.49948 41.8958)
3 POINT (12.4737 41.89467)
head(df3, n=3)
name addr.street
1 Pizzeria ai Marmi Viale di Trastevere
2 Sichuan Haozi Via di San Martino ai Monti
3 Dar filettaro a Santa Barbara Largo dei Librari