1 Load libraries

rm(list = ls())

library(dplyr) # data wrangling
library(tidyr) # data wrangling
library(ggplot2) # data visualisation
library(sf) # simple features - geospatial geometries
library(osmdata) # obtaining OpenStreetMap vector data
library(units) # working with units
library(mapview) # interactive geometry viewing
library(ggmap) # downloading raster maps from a variety of sources
library(ggspatial) # map backgrounds and annotations for ggplot
library(tmap) # static/interactive map library with ggplot-like syntax

2 Read the file

greensp_sf <- readRDS("greensp_sf.rds")
head(greensp_sf)
## Simple feature collection with 6 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -3.388398 ymin: 55.90319 xmax: -3.140106 ymax: 55.98513
## Geodetic CRS:  WGS 84
##                           name  osm_id     leisure
## 4271400            Dundas Park 4271400        park
## 4288119     Harrison Park East 4288119        park
## 4348244      Bruntsfield Links 4348244        park
## 4891768   Baberton Golf Course 4891768 golf_course
## 4891786 Kingsknowe Golf Course 4891786 golf_course
## 4892551                   <NA> 4892551        <NA>
##                               geometry
## 4271400 MULTIPOLYGON (((-3.387083 5...
## 4288119 MULTIPOLYGON (((-3.224193 5...
## 4348244 MULTIPOLYGON (((-3.203498 5...
## 4891768 MULTIPOLYGON (((-3.28719 55...
## 4891786 MULTIPOLYGON (((-3.265052 5...
## 4892551 MULTIPOLYGON (((-3.146338 5...
unique(greensp_sf$leisure)
## [1] "park"           "golf_course"    NA               "nature_reserve"

3 Coordinate reference systems (CRS)

st_crs(greensp_sf)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   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["unknown"],
##         AREA["World"],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

Let’s transform our data into the British National Grid CRS, using it’s ESPG number (27700):

greensp_sf <- st_transform(greensp_sf, 27700)
st_crs(greensp_sf)
## Coordinate Reference System:
##   User input: EPSG:27700 
##   wkt:
## PROJCRS["OSGB36 / British National Grid",
##     BASEGEOGCRS["OSGB36",
##         DATUM["Ordnance Survey of Great Britain 1936",
##             ELLIPSOID["Airy 1830",6377563.396,299.3249646,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4277]],
##     CONVERSION["British National Grid",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",49,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-2,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996012717,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",400000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",-100000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."],
##         BBOX[49.75,-9,61.01,2.01]],
##     ID["EPSG",27700]]

4 Spatial Operations

4.1 Calculate Each Area

4.1.1 Area calculation

Calculating the area of polygon or multipolygon geometries is done using the st_area() function:

# Create and calculate a new column for feature area
greensp_sf <- mutate(greensp_sf, area = st_area(greensp_sf))

head(greensp_sf$area)
## Units: [m^2]
## [1]  29515.82  29073.37 150754.53 477437.67 398419.52 313020.65

4.1.2 Convert m2 to hectares

# convert m2 to hectares 
greensp_sf <- greensp_sf %>% 
  mutate(area_ha = set_units(area, "ha")) %>% 
  select(-area) # drop "area" column

4.1.3 View sf object

mapview(greensp_sf)

4.1.4 Clean data

# Filter out unneeded shapes
greensp_sf <-
  greensp_sf %>%
  filter(is.na(leisure) == FALSE) %>%  
  # remove leisure NAs
  rename(greensp_type = leisure) %>%
  # rename leisure to greensp_type
  st_make_valid()
  # a good function to use after importing data to make sure shapes are valid
# Remove green spaces with <2 ha
greensp_sf <- filter(greensp_sf, as.numeric(area_ha) >= 2)

4.1.5 View sf object

mapview(greensp_sf)

4.1.6 Split sf by green space type

# Separate into a list of multiple sf grouped by type
greensp_sf_list <- greensp_sf %>% split(.$greensp_type)

# Each sf object in the list can be accessed using the $ operator.
# E.g. greensp_sf_list$nature_reserve to get the sf object
# containing only nature reserves.

4.2 Remove overlap

4.2.1 Remove overlap

Using st_difference(x, y)

# Remove the parts of parks where they are overlapped by nature reserves
greensp_sf_list$park <- st_difference(greensp_sf_list$park,
                                      st_union(greensp_sf_list$nature_reserve))
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
# Remove the parts of parks where they are overlapped by golf courses
greensp_sf_list$park <- st_difference(greensp_sf_list$park,
                                      st_union(greensp_sf_list$golf_course))
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

4.2.2 Calculate Area

# Let's turn the list back into one sf object. We will also need to re-calculate area!
greensp_sf <- bind_rows(greensp_sf_list) %>%
  # bind the list into one sf object
  mutate(area_ha = set_units(st_area(.), "ha")) %>%
  # calculate area again
  filter(as.numeric(area_ha) >= 2) # remove area < 2 ha again

5 Draw maps

5.1 Data preparation

# Reorder greenspace types, capitalise, remove underscores
greensp_sf_forplot <-
  mutate(greensp_sf,
         greensp_type = factor(greensp_type,
                               levels = c("park", "nature_reserve", "golf_course"),
                               labels = c("Park", "Nature reserve", "Golf course")))

5.2 Static map

5.2.1 Download raster background for map

# Get the rectangular bounding box
city_rect <- getbb("City of Edinburgh", featuretype = "settlement")

stamen_raster <- get_stamenmap(city_rect, zoom = 12)

5.2.2 Plot with ggplot2

# Plot map with ggplot
(edi_greenspaces_map <-
  ggplot(data = greensp_sf_forplot) +
    inset_ggmap(stamen_raster) + # add ggmap background
    geom_sf(aes(fill = greensp_type)) + # add sf shapes, coloured by greensp_type
    coord_sf(crs = st_crs(4326), expand = FALSE) +
    # change the CRS of the sf back to WGS84 to match the ggmap raster
    scale_fill_manual(values = c("#44AA99", "#117733", "#AA4499")) +
    # add custom colours from Tol palette (colourblind-friendly)
    labs(title = "Green spaces in Edinburgh, Scotland",
         subtitle = "Parks, nature reserves and golf courses > 2 ha\n",
         caption = paste0("Map tiles by Stamen Design (stamen.com), CC BY 3.0. ",
                         "http://creativecommons.org/licenses/by/3.0\n",
                         "Map data © OpenStreetMap contributors, ODbL. ",
                         "http://www.openstreetmap.org/copyright")) +
    # add various labels
    annotation_scale(location = "bl") + # ggspatial scale on bottom left
    annotation_north_arrow(location = "tr") + # ggspatial arrow on top right
    theme_void() + # get rid of axis ticks, titles
    theme(legend.title = element_blank(),
          legend.position = c(.98, .02),
          legend.justification = c("right", "bottom"),
          legend.box.just = "right",
          legend.box.background = element_rect(fill = "white", colour = "gray"),
          legend.margin = margin(6, 6, 6, 6),
          # move legend to bottom right and customise
          plot.margin = margin(12,12,12,12))
          # add margin around plot
)

5.2.3 Image Save

ggsave("/Users/admin/Desktop/Linh Data Studio/Spatial Data/Spatial Modelling/edi_greenspaces_map.png", edi_greenspaces_map, width = 8, height = 6.5)

5.3 Interactive map with tmap

# Plot interactively with tmap
tmap_mode("view") # interactive mode

(edi_greenspace_tmap <-
  tm_basemap("Stamen.Terrain") + # add Stamen Terrain basemap
  tm_shape(greensp_sf_forplot) + # add the sf
  tm_sf(col = "greensp_type", # colour by green space type
        title = "", # no legend title
        palette = c("#44AA99", "#117733", "#AA4499"), # custom fill colours
        popup.vars = c("Area  " = "area_ha"), # customise popup to show area
        popup.format = list(digits=1)) + # limit area to 1 decimal digit
  tm_scale_bar() # add scale bar
)
tmap_save(tm = edi_greenspace_tmap, filename = "/Users/admin/Desktop/Linh Data Studio/Spatial Data/Spatial Modelling/edi_greenspace_tmap.html")

6 Calculate total area

Calculate the total area of each type of green space. Please note that we can’t sum the areas we calculated as there are overlapping polygons within the green space types.

(greensp_type_area <-
  lapply(greensp_sf_list,
      function(x) set_units(st_area(st_union(x)), "ha")) %>% # apply the function to each element of the list using lapply()
  as.data.frame() %>%
  pivot_longer(cols = everything(), names_to = "greensp_type",
               values_to = "area_ha"))
## # A tibble: 3 × 2
##   greensp_type   area_ha
##   <chr>             [ha]
## 1 golf_course       907.
## 2 nature_reserve    290.
## 3 park             1387.