Introduction

You can map data from the Improvement Service Spatial Hub using R. Below are examples showing how this can be achieved using WFS pathways and shapefiles.

Web Feature Services (WFS) are a means of accessing geographic data via a file pathway, drawing directly from the Spatial Hub. The advantage of utilising WFS files lies in their ability to fetch data in real-time from the Spatial Hub each time you refresh your map. This approach is more dynamic compared to using static file types that are downloaded and stored on your local device (such as those in your downloads folder). Static files risk becoming outdated unless you frequently check the source data and manually update your files with the latest data versions.

Reading Shapefiles

You can read shapefiles from the Spatial Hub using the sf library in R: https://r-spatial.github.io/sf/

Shapefiles in the Spatial Hub are stored with the British National Grid projection. For web mapping, re-project these to WGS84, as shown in the code snippet below.

# add packages to R
if (!require(dplyr)) install.packages(dplyr)
if (!require(sf)) install.packages(sf)
if (!require(leaflet)) install.packages(leaflet)

#load packages
library(dplyr) # a package for data wrangling
library(sf) # a package for handling GIS data
library(leaflet) # a package for mapping GIS data

# once I have downloaded my shapefile from the spatial hub, I can read it into R using st_read.
town_centre <- st_read("C:/Users/kirkland.aline/OneDrive - IS/Documents/pub_townc.shp")

## Reading layer `pub_townc' from data source 
##   `C:\Users\kirkland.aline\OneDrive - IS\Documents\pub_townc.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 397 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 124603.3 ymin: 540130.3 xmax: 413642.6 ymax: 1011389
## Projected CRS: OSGB36 / British National Grid

# I now want to change the projection so I can map the data. Transform to WGS84 (EPSG:4326)
town_centre <- st_transform(town_centre, crs = 4326)

# I want to check the geometery type. For mapping, I need a simple geometry type like point or polygon
st_geometry_type(town_centre, by_geometry = FALSE)

## [1] MULTIPOLYGON
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

# The shapefile has multipolygon geometry, so I need to change that to polygon so I can map the data. 
town_centre <- st_cast(town_centre, "POLYGON")

# Leaflet library allows you to create maps from the shapefile.  
leaflet()%>%
  addProviderTiles('Esri.WorldImagery', options = providerTileOptions(), group = "Imagery Basemap") %>% # add an ESRI imagery basemap. 
  addPolygons(data = town_centre, # add a polygon of the shapefile to the map
              color = "pink")

Reading WFS pathways

# add packages to R
if (!require(dplyr)) install.packages(dplyr)
if (!require(sf)) install.packages(sf)
if (!require(leaflet)) install.packages(leaflet)

#load packages
library(dplyr) # a package for data wrangling
library(sf) # a package for handling GIS data
library(leaflet) # a package for mapping GIS data

# read in the WFS from the Spatial Hub
tc_wfs <- st_read("WFS:https://geo.spatialhub.scot/geoserver/sh_townc/wfs?service=WFS&authkey=your_auth_key_here")

## Reading layer `sh_townc:pub_townc' from data source 
##   `WFS:https://geo.spatialhub.scot/geoserver/sh_townc/wfs?service=WFS&authkey=your_auth_key_here' 
##   using driver `WFS'
## Simple feature collection with 397 features and 10 fields
## Geometry type: MULTISURFACE
## Dimension:     XY
## Bounding box:  xmin: 124603.3 ymin: 540130.3 xmax: 413642.6 ymax: 1011389
## Projected CRS: OSGB36 / British National Grid

# I want to check the geometery type. For mapping, I need a simple geometry type like point or polygon
st_geometry_type(tc_wfs, by_geometry = FALSE)

## [1] MULTISURFACE
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

# the geometry type was multi-surface. I need to transform this first to a geometery collection and then to a simple polygon to make it mappable.
tc_wfs <- sf::st_cast(tc_wfs, "GEOMETRYCOLLECTION")

# extract polygon data from geometry collection
tc_wfs <- sf::st_collection_extract(tc_wfs, type = c("POLYGON"))

# re-project the data from British National Grid to WGS84 for web mapping
tc_wfs <- sf::st_transform(tc_wfs, crs = 4326)

# cast the data to a simple polygon for mapping
tc_wfs <- st_cast(tc_wfs, "POLYGON")

# Leaflet library allows you to create maps from the shapefile.  
leaflet()%>%
  addProviderTiles('Esri.WorldImagery', options = providerTileOptions(), group = "Imagery Basemap") %>% # add an ESRI imagery basemap. 
  addPolygons(data = tc_wfs, # add a polygon of the wfs to the map
              color = "red")