1 Set up

Install necessary packages

install.packages("sf", dependencies = TRUE)
install.packages("terra", dependencies = TRUE)
install.packages(c("spData", "mapview", "ggplot2", "tmap"))
install.packages("tidyverse", dependencies = TRUE)

2 Working with vectors

Load packages

library(sf)
library(spData)

2.1 Points

Example dataset: Bicycle rental points across London (Source: OSM)

data("cycle_hire_osm", package = "spData")
head(cycle_hire_osm) # plot first few rows
## Simple feature collection with 6 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -0.1293092 ymin: 51.52583 xmax: -0.090836 ymax: 51.53402
## Geodetic CRS:  WGS 84
##      osm_id                       name capacity cyclestreets_id description
## 1    108539            Windsor Terrace       14            <NA>        <NA>
## 2 598093293 Pancras Road, King's Cross       NA            <NA>        <NA>
## 3 772536185 Clerkenwell, Ampton Street       11            <NA>        <NA>
## 4 772541878                       <NA>       NA            <NA>        <NA>
## 5 781506147                       <NA>       NA            <NA>        <NA>
## 6 783824668      Finsbury Library, EC1       NA            <NA>        <NA>
##                      geometry
## 1 POINT (-0.0933878 51.52913)
## 2 POINT (-0.1293092 51.53402)
## 3 POINT (-0.1182352 51.52729)
## 4  POINT (-0.090836 51.52583)
## 5 POINT (-0.1210572 51.53001)
## 6 POINT (-0.1038272 51.52594)

Visualise the points in the geometry column

plot(cycle_hire_osm["geometry"], axes = TRUE)

Examine the CRS

st_crs(cycle_hire_osm)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

Transform the CRS

cycle_hire_osm <- st_transform(cycle_hire_osm, crs = st_crs(32630)) # UTM zone 30N for London
head(cycle_hire_osm)
## Simple feature collection with 6 features and 5 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 699095.4 ymin: 5712291 xmax: 701799.3 ymax: 5713119
## Projected CRS: WGS 84 / UTM zone 30N
##      osm_id                       name capacity cyclestreets_id description
## 1    108539            Windsor Terrace       14            <NA>        <NA>
## 2 598093293 Pancras Road, King's Cross       NA            <NA>        <NA>
## 3 772536185 Clerkenwell, Ampton Street       11            <NA>        <NA>
## 4 772541878                       <NA>       NA            <NA>        <NA>
## 5 781506147                       <NA>       NA            <NA>        <NA>
## 6 783824668      Finsbury Library, EC1       NA            <NA>        <NA>
##                   geometry
## 1 POINT (701607.8 5712674)
## 2 POINT (699095.4 5713119)
## 3 POINT (699892.7 5712401)
## 4 POINT (701799.3 5712314)
## 5 POINT (699685.1 5712696)
## 6 POINT (700897.9 5712291)

Interactive map of the ‘capacity’ of each bicycle rental point using the package mapview

library(mapview)
mapview(cycle_hire_osm["capacity"], layer.name = "Rental capacity")

2.2 Lines

Example dataset: A river network in France

data("seine", package = "spData")
head(seine)
## Simple feature collection with 3 features and 1 field
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: 518344.7 ymin: 6660431 xmax: 879955.3 ymax: 6938864
## Projected CRS: RGF93 / Lambert-93
##    name                       geometry
## 1 Marne MULTILINESTRING ((879955.3 ...
## 2 Seine MULTILINESTRING ((828893.6 ...
## 3 Yonne MULTILINESTRING ((773482.1 ...

Plot the lines using ggplot2

library(ggplot2)

ggplot(seine) +
  geom_sf(aes(color = name)) +
  geom_sf_label(aes(label = name)) # add labels for river names

2.3 Polygons

Example dataset: Country polygons worldwide

data(world, package = "spData")
head(world)
## Simple feature collection with 6 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 11
##   iso_a2 name_long  continent region_un subregion type  area_km2     pop lifeExp
##   <chr>  <chr>      <chr>     <chr>     <chr>     <chr>    <dbl>   <dbl>   <dbl>
## 1 FJ     Fiji       Oceania   Oceania   Melanesia Sove…   1.93e4  8.86e5    70.0
## 2 TZ     Tanzania   Africa    Africa    Eastern … Sove…   9.33e5  5.22e7    64.2
## 3 EH     Western S… Africa    Africa    Northern… Inde…   9.63e4 NA         NA  
## 4 CA     Canada     North Am… Americas  Northern… Sove…   1.00e7  3.55e7    82.0
## 5 US     United St… North Am… Americas  Northern… Coun…   9.51e6  3.19e8    78.8
## 6 KZ     Kazakhstan Asia      Asia      Central … Sove…   2.73e6  1.73e7    71.6
## # ℹ 2 more variables: gdpPercap <dbl>, geom <MULTIPOLYGON [°]>

Plot the per capita GDP across all countries using an appropriate colour scale

breaks <- quantile(world$gdpPercap, probs = c(0.025, 0.25, 0.5, 0.75, 0.975), na.rm = TRUE) # get quantile values

ggplot(world) + geom_sf(aes(fill = gdpPercap)) + 
  scale_fill_gradientn(colours= c("#ffffcc", "#a1dab4" , "#41b6c4", "#2c7fb8", "#253494"), # binned color scale
                       trans = 'log', # log-transformation
                       breaks = breaks) +
  theme_minimal()

Use the package tmap to quickly plot thematic and customisable maps

library(tmap)

tmap_mode("view")
tm_shape(world) +
  tm_polygons("gdpPercap", palette = "YlGnBu", style = "quantile") +
  tm_view(set.view = 1) # set zoom level

2.4 Geometry operations

Locate the center point of geometries

world_pts <- st_centroid(world)

Simplify shapes by removing vertices

world_projected <- st_transform(world, "EPSG: 3857") # need to transform to projected CRS
world_simple <- st_simplify(world_projected, dTolerance = 1e5) # 100km simplification

Cast the simplified POLYGON objects back to type ‘MULTIPOLYGON’, and remove the rows with empty geometries

world_simple <- st_cast(world_simple, "MULTIPOLYGON")
world_simple <- world_simple[!st_is_empty(world_simple),] # subset rows

Extend the boundaries of vectors

world_buffered <- st_buffer(world_simple, dist = 1e5) # 100km

Remove overlaps among the polygons in world_simple and merge all polygons

world_diff <- st_difference(world_simple) # exclude overlaps
world_union <- st_union(world_simple)

Plot the manipulated objects derived from the country polygons

tmap_options(check.and.fix = TRUE)
sf_use_s2(FALSE)

tm_shape(world_pts) + tm_dots() + # polygon centroids
  tm_shape(world_simple) + tm_polygons(alpha = 0.5) + # simplified
  tm_shape(world_buffered) + tm_polygons(alpha = 0.5) + # buffered
  tm_shape(world_union) + tm_polygons(alpha = 0.5) + # merged (union)
  tm_view(set.view = 2) # zoom in

2.5 Spatial operations

Subset rows of an sf (‘target’) object based on its spatial relationship with another

uk <- st_filter(world_projected, # target: world polygons
                st_transform(cycle_hire_osm, st_crs(world_projected)), # bicycle rental points in London - crs must be the same
                .predicate = st_intersects) # default setting

tm_shape(uk) + tm_polygons() +
  tm_shape(cycle_hire_osm) + tm_dots()

Match rows (and append columns) to an sf (‘target’) object based on its spatial relationship with another

cycle_hire_osm <- st_join(cycle_hire_osm, # target
                          st_transform(world_projected, st_crs(cycle_hire_osm)))

head(cycle_hire_osm)
## Simple feature collection with 6 features and 15 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 699095.4 ymin: 5712291 xmax: 701799.3 ymax: 5713119
## Projected CRS: WGS 84 / UTM zone 30N
##      osm_id                       name capacity cyclestreets_id description
## 1    108539            Windsor Terrace       14            <NA>        <NA>
## 2 598093293 Pancras Road, King's Cross       NA            <NA>        <NA>
## 3 772536185 Clerkenwell, Ampton Street       11            <NA>        <NA>
## 4 772541878                       <NA>       NA            <NA>        <NA>
## 5 781506147                       <NA>       NA            <NA>        <NA>
## 6 783824668      Finsbury Library, EC1       NA            <NA>        <NA>
##   iso_a2      name_long continent region_un       subregion    type area_km2
## 1     GB United Kingdom    Europe    Europe Northern Europe Country 249986.4
## 2     GB United Kingdom    Europe    Europe Northern Europe Country 249986.4
## 3     GB United Kingdom    Europe    Europe Northern Europe Country 249986.4
## 4     GB United Kingdom    Europe    Europe Northern Europe Country 249986.4
## 5     GB United Kingdom    Europe    Europe Northern Europe Country 249986.4
## 6     GB United Kingdom    Europe    Europe Northern Europe Country 249986.4
##        pop  lifeExp gdpPercap                 geometry
## 1 64613160 81.30488  38251.79 POINT (701607.8 5712674)
## 2 64613160 81.30488  38251.79 POINT (699095.4 5713119)
## 3 64613160 81.30488  38251.79 POINT (699892.7 5712401)
## 4 64613160 81.30488  38251.79 POINT (701799.3 5712314)
## 5 64613160 81.30488  38251.79 POINT (699685.1 5712696)
## 6 64613160 81.30488  38251.79 POINT (700897.9 5712291)

Combine the official dataset cycle_hire with the OSM dataset cycle_hire_osm of bicycle rental points across London, but only those beyond a distance of 10 m (remove duplicate points)

3 Working with rasters

Load package

library(terra)

3.1 Continuous

Example raster: Elevation in Luxembourg

filepath <- system.file("ex/elev.tif", package = "terra")
elevation <- rast(filepath)
elevation_utm <- project(elevation, "EPSG:23032") # project to UTM zone 32N for Luxembourg

Example vector: Polygons of Luxembourg

filepath <- system.file("ex/lux.shp", package = "terra")
luxembourg <- vect(filepath)
luxembourg
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 12, 6  (geometries, attributes)
##  extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
##  source      : lux.shp
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP
##  type        : <num>    <chr> <num>    <chr> <num> <int>
##  values      :     1 Diekirch     1 Clervaux   312 18081
##                    1 Diekirch     2 Diekirch   218 32543
##                    1 Diekirch     3  Redange   259 18664

Summarise the mean raster values for each sub-region (polygon) in Luxembourg

elev_mean <- extract(elevation, luxembourg, 
                     fun = mean, na.rm = TRUE)
luxembourg$elevation_mean <- elev_mean$elevation # add results to new column

Plot the rasters and polygons

tm_shape(elevation) +
  tm_raster(palette = "RdYlGn", n = 10) + # use similar color palette
  tm_shape(st_as_sf(luxembourg)) + # convert to sf, as tmap currently does not support SpatVect
  tm_polygons(col = "elevation_mean", palette = "RdYlGn", alpha = 0.7)

Narrow down raster data to a specific area of interest

luxembourg_subset <- luxembourg[luxembourg$NAME_2 == "Mersch"]
elevation_subset <- crop(elevation, luxembourg_subset) # reduce spatial extent
elevation_subset <- mask(elevation_subset, luxembourg_subset) # convert values outside polygon to NA
plot(elevation_subset)

3.2 Discrete

Classify values into a discrete (categorical) raster

elevation_2class <- ifel(elevation > 400, 1, 0) # E.g., threshold for 'high' elevation at 400m
plot(elevation_2class) # low/high

elevation_3class <- 
  classify(elevation, 
           rbind(c(400, Inf, 2), # matrix of 3 columns, 'from', 'to' and 'becomes'
                 c(300, 400, 1),
                 c(-Inf, 300, 0)))  
plot(elevation_3class) # low/mid/high

4 Data conversions

Vectorize rasters into points, lines and polygons

elevation_pts <- as.points(elevation_subset) # points
elevation_contours <- as.contour(elevation_subset) # lines
elevation_polygons <- as.polygons(elevation_subset) # polygons

Visualise all output

tm_shape(elevation_subset) + tm_raster(palette = "RdYlGn", n = 10) + # original raster
  tm_shape(st_as_sf(elevation_pts)) + tm_dots() + 
  tm_shape(st_as_sf(elevation_contours)) + tm_lines() + 
  tm_shape(st_as_sf(elevation_polygons)) + tm_polygons(alpha = 0) + 
  tm_view(set.view = 11)

Rasterize the points dataset cycle_hire_osm based on the column ‘capacity’

# first create a raster template at a specified resolution and crs
rast_template <- rast(ext(cycle_hire_osm), 
                      resolution = 1000,
                      crs = st_crs(cycle_hire_osm)$wkt) # input format is WKT

# then, rasterize the points dataset based on the template
cyclehire_capacity <- 
  rasterize(cycle_hire_osm, rast_template,
            field = "capacity", 
            fun = sum, # add up values of points within each pixel
            na.rm = TRUE)

Overlay the original points over the per-pixel capacity of bicycle docking stations

tm_shape(cyclehire_capacity) + 
  tm_raster(title = "Capacity (sum)", palette = "YlGn") +
  tm_shape(st_as_sf(cycle_hire_osm)) + 
  tm_dots(size = 0.01)

© XP Song
Disclaimer: These notes exclusively reserved for use with the Geospatial Data Science with R online course. Please do not distribute or reproduce these materials, they are intended solely for personal use.