17 GIS with R

As a rule of thumb, use GIS with R if you’re familiar with R but not ArcGIS and Python, and your project has just some spatial content. If your project is predominantly spatial, you should consider upskilling to ArcGIS. The GIS capabilities of R packages are increasingly competent. However, you’ll be swimming against the tide if you don’t take advantage of the Council’s investment in the ArcGIS platform, tools and skills.

This chapter focusses on how you can include spatial data in your R project from public data sources and the Council’s ArcGIS platform.

17.1 Features

The term “feature” is useful to distinguish a record as “spatial”. By “spatial” we mean it has a point, a line or a polygon. A point is a pair of coordinates, a line a minimum of two pairs of coordinates and polygons at least three pairs.

Something to keep an eye out for is MULTIPOINT, multi-part polygons etc. This is a single record, a single feature, but it has more than one point, line or polygon. For example, Sheffield secondary school features might be MULTIPOINT so that the feature for King Edward VII school can have one point for the Lower School near Crosspool, and another point for the Upper School in Broomhill.

17.1.1 Simple features

Simple features is an open standard developed by the Open Geospatial Consortium. It is a hierarchical data model that represents a wide range of geometry types, all of which are supported by the sf R package. The sf package also plays nice with the tidyverse.

# Create a data frame with some example records 
df <- tibble::tribble(
  ~name,    ~postcode, ~longitude,  ~latitude,
  "SCC",    "S1 2HH",    -1.470006,   53.38038,
  "Blades", "S2 4SU",  -1.470512,   53.36986,
  "Owls",   "S6 1SW",    -1.500859,   53.41084
)

# Create simple features
sf_shef <- sf::st_as_sf(df, 
                        coords = c("longitude", "latitude"), #from postcodes.io
                        crs = 4326) #WSG84

# Plot
mapview::mapview(sf_shef, layer.name = "Sheffield institutions")

## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)

17.2 Coordinates & projections

The minimum you need to know is that any spatial data you use is likely to involve one or more of three coordinate systems:

WSG84 with an ESPG of 4326
Web Mercator with an EPSG of 3857
OSGB1936 (British National Grid) with an ESPG of 27700

Most web maps are based on features stored as WSG84 and displayed (projected) as Web Mercator. Postcodes.io uses WSG84, but the Council’s Portal and AGOL, and the ONS Open Geography Portal use OSGB1936 (BNG).

# Check the coordinate system of simple features
sf::st_crs(sf_shef)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     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]]

17.2.1 Geographic & projected coordinate systems

A reference ellipsoid as an approximation of the surface of the earth and provides a shape on which a coordinate system can be placed.

Geographic coordinate systems use the ellipsoid to define specific locations on the surface to create a grid. Datums are geographic coordinate systems based on a specific ellipsoid (so a more specific geographic coordinate system), with an origin at a specific location, and the ellipsoid at a specific orientation. These are also called “spatial reference systems” or “coordinate reference systems”.

Projected coordinate systems are like geographic coordinate systems. A projected coordinate system is also a grid used as a reference for locations on the planet, but it’s a translation of the three-dimensional grid onto a two-dimensional plane such as a paper map or screen.

A GCS defines where the data is located on the earth’s surface. A PCS tells the data how to draw on a flat surface such as a screen.

Coordinate systems, projected and geographic, are often identified by an EPSG code.

17.3 ArcGIS R & Python packages

Listed below are R and Python packages that allow us to leverage ArcGIS functionality. Generally speaking, they are in descending order of ease of use, and the scope of functionality is in ascending order.

Package	Language	ArcGIS Desktop Install Required	Scope
`esri2sf`	R	✕	Read from AGOL & ArcGIS Enterprise
`arcgisbinding`	R	✓	Read & write from AGOL & ArcGIS Enterprise
`arcgis`	Python	✓	Access AGOL & ArcGIS Enterprise functionality
`arcpy`	Python	✓	Access ArcGIS Desktop or Pro functionality

17.3.1 esri2sf

The esri2sf R package is the only one of the four packages listed above that is not maintained by ESRI. It is the simplest to start using because it doesn’t require an installation and license for ArcGIS Desktop or Pro. The flip side, is that it has the least amount of functionality and the code isn’t as robust. The scope of the package is limited to reading simple features or data frames from AGOL & ArcGIS Enterprise (ours or other organisation’s).

The ONS Open Geography Portal is based on the same ArcGIS web tools as our Portal and AGOL. To get, for example, the local authority boundary for Sheffield you first need to navigate to the dataset on the website and get the GeoService API URL:

Then, using:

the copied GeoService URL,

https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/Local_Authority_Districts_December_2020_UK_BGC/FeatureServer/0/query?outFields=\\\\\\\*&where=1%3D1
some resources to understand how to query the API,

ArcGIS REST API - Query (Feature Service/Layer)

SQL 92 is used by ESRI based APIs

HTML encoding is used for building URL queries
and esri2sf,

we can get the Sheffield boundary as a simple feature:

# Base URL for ONS Open Geography Portal
ons_geog_base_url <- "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"

# URL part for how detailed the local authority boundary should be
generalised <- "Local_Authority_Districts_December_2020_UK_BGC"

# Get Sheffield MSOA boundaries as simple features  
shef_msoa <- esri2sf::esri2sf(url = stringr::str_c(ons_geog_base_url,
                                                   generalised,
                                                   "/FeatureServer/0"),
                              where = "LAD20NM LIKE 'Sheffield%'",
                              geomType = "esriGeometryPolygon")

# Plot
mapview::mapview(shef_msoa, layer.name = "Sheffield")

The example above reads open data. However, most of the Council’s spatial data on Portal and AGOL requires login credentials to access. Nevertheless, Portal and AGOL login credentials are a smaller requirement than a ArcGIS Desktop or Pro installation and license.

TODO - test and note authentication with token

Our public github.com/scc-pi/functions repository includes functions in ShefAreas.R that use esri2sf. The sub-headings immediately below outline how other ArcGIS packages are more powerful, but also more onerous in terms of licensing and setup. So where possible, the intention is to utilise esri2sf in preference to the other ArcGIS packages.

Further information on the esri2sf package, how to install and use it, is available from github.com/yonghah/esri2sf.

17.3.2 R-ArcGIS Bridge

The R-ArcGIS Bridge is the name for the arcgisbinding R package. It’s scope is limited to:

Transferring data in both directions between ArcGIS and R
Wrapping R tools for use in ArcGIS

For more information on using the arcgisbinding package see the vignette on Using the R-ArcGIS Bridge.

To use the R-ArcGIS Bridge you need an installation and license for ArcGIS Desktop or Pro. After loading the arcgisbinding package you have to call arc.check_product() to define a desktop license:

library(arcgisbinding)
arc.check_product()

Using arcgisbinding with ArcGIS Desktop means that you need to use 32-bit R, because ArcGIS Desktop is 32-bit. An ArcGIS Desktop 64-bit Background Geoprocessing option is a red herring because it only allows using the bridge from ArcGIS, not from within R itself. ArcGIS Pro is 64-bit and allows you to use use 64-bit R with arcgisbinding.

The README.md in the r-bridge GitHub repository is the best starting point for more detail on R-ArcGIS Bridge requirements and installation.

The R-ArcGIS Bridge cannot write spatial data directly to the Council’s Portal and AGOL. It can write to a local geodatabase, which can then be published to Portal or AGOL via ArcGIS Pro. The R-ArcGIS Bridge can read spatial data directly from the Council’s Portal and AGOL:

The arc.check_product() to define a desktop license means that the results of the example above can’t be rendered via GitHub Actions in these notes. However, the same example has been knitted locally and included in some notes on isochrones.

17.4 Leaflet

TODO

17.5 Rasters

TODO

17.6 Further resources

Introduction to GIS in R, 2-3hr course by ONS Geo (geography teams in ONS)

Geocomputation with R, by Robin Lovelace, Jakub Nowosad, & Jannes Muenchow

Spatial Data Science with applications in R, by Edzer Pebesma & Roger Bivand

mapview basics vignette

Geographic vs Projected Coordinate Systems, Heather Smith, ArcGIS Blog (27 February 2020)

16 GIS with Python

18 Postcodes & addresses