\(~\)

Hello everyone 👋,

we hope you have enjoyed our small presentation, refreshed your knowledge about Coordinate Reference Systems and got some insights on how geographic data they works in R and what potential data operations can be done with it.

Ready to play around with sf and do some simple operations and create a nice map in the end (hopefully) 😍?

In this little exercise we would like you to

  • get a better feeling for CRS and the sf- package
  • understand how to access and read CRS data in R
  • do simple data operations
  • learn how to create cool and interactive maps based on geographic data

Theory

Refreshing you memory from the presentation…

Exercise (1)

In your own words, why do we need a “Coordinate Reference System” (CRS)? 

#
#

Setup

The first step is to download the concerning packages and import the libraries for using sf and the spatial data.

# install.packages("sf", "spData", "spDataLarge")

# The Essentials
library(sf)   # classes and functions for vector data
library(tidyverse)     # for tidy data and ggplot2

# Datasets
library(rgdal) # dataset of city centers (points)
library(spData) # load explanatory geographic data
library(spDataLarge) # load larger geographic data 

#(Intermediary Mapping)
library(tmap) # to create nice and interactive maps
library(jsonlite) # Convert R objects to/from JSON

Be careful: There might be a little bug in installing sf and spData. If you’re running Mac or Linux, installing sf may not work directly. Please have a look here to download the relevant packages. Then use the library() command to load them.

Exploring CRS in R

For these exercises we will work with the same dataset from the video. Our new object, new_vector, will be a polygon representing a world map. (For more info: ?spData::world).

First import the data and make it a readable to sf:

vector_filepath <- system.file("shapes/world.gpkg", package = "spData")
new_vector <- read_sf(vector_filepath)

We’ve also loaded in a dataset of city centers from the rgdal package as the object cities.

Retrieve the CRS from objects new_vector and cities

Exercise (2)

** Using st_crs what can you tell about the data? Ex. the Datums and the retrieved Coordinate Reference System**

#
#

new_vector

st_crs(new_vector)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     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]]

cities

st_crs(cities)
## Coordinate Reference System:
##   User input: EPSG:3138 
##   wkt:
## PROJCRS["ETRS89 / ETRS-GK31FIN",
##     BASEGEOGCRS["ETRS89",
##         DATUM["European Terrestrial Reference System 1989",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["Finland ETRS-GK31",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",31,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Cadastre, engineering survey, topographic mapping (large scale)."],
##         AREA["Finland - east of 30°30'E."],
##         BBOX[62.08,30.5,64.27,31.59]],
##     ID["EPSG",3138]]

\(~\)

Use st_set_crs to Add or Change CRS

CRS_vector <- st_set_crs(world, "EPSG:4326")

To look up different Datums you can use a website for Coordinate Systems Worldwide

Exercise (3)

How can we change the Datum for cities?

#
#

\(~\)

Geographic Data Operations

Geometric measurements

In the short presentation we have shown you the following examples of how to conduct two types of geometric measurements with the sf-package

Example 1 from presentation:

Belgium <- world %>% 
  filter(name_long == "Belgium")
st_area(Belgium) 
## 30019548314 [m^2]
units::set_units(st_area(Belgium), km^2)
## 30019.55 [km^2]

We can also do this calculation for a US-state:

Florida <- us_states %>% 
  filter(NAME == "Florida")
st_area(Florida) 
## 1.51264e+11 [m^2]

Exercise (4)

st_area() is giving us the output in m^2. The output is however really difficult to read, and the US does operate with miles. How can we get an output in squared miles?

#
#

Example 2 from presentation:

Measuring the distance from a cycle hire point at the Palace Gate in London to the center of Redbridge.

cycle_station_palace <- cycle_hire %>% 
  filter(name == "Palace Gate")
redbridge <- lnd %>% 
  filter(NAME == "Redbridge")
redbridge_centroid <- st_centroid(redbridge)
st_distance(cycle_station_palace, redbridge_centroid)
## Units: [m]
##         [,1]
## [1,] 20258.6

Exercise (5)

Have a look at the Code. What are the three steps one has to make, to measure distance between two geometric data points? Can we change the units of the output here, too? How?

#
#
  • Tip: other interesting functions to explore geometric relations and distances with are st_within(),st_is_within_distance() or st_touches(). If you would like to find out what they do have a look here.

Creating maps wit CRS data 🗺

In the presentation we have shown how to use sf in combination with ggplot or the sf and plot() function to create simple maps.

Example with ggplot()

Example from the presentation using ggplot() in combination with geom_sf() from the sf package on data CRS data from the US retrieved from the spData package.

Example with plot()

Example from the presentation using plot() on data CRS data from the US and attributed data on the median income per state, both retrieved from the spData package.

\(~\)

Public Policy implications

Exercise (6)

During this workshop we’ve covered what CRS data is, the functions of the sf package, and how it can be used. Now brainstorm some applications; What are scenarios (how or when) where you can apply these tools in the future?

#
#

Take-home exercise: Introducing tmap()

Exercise (7)

In the presentation we mentioned the tmap package which offers a bit more comprehensive tools, for example interactive mapping. Please do a quick research. Which function from the tmap-package can we use to set an interactive view of a map?

#
#

\(~\)

Sources

The workshop and these exercises are based on Geocompuataion with R Science, section 2.4, tmap: get started!