This script provides an example of downloading and importing administrative boundaries from the Office for National Statistics Open Geography portal into RStudio and plotting on a map. It is based upon a useful example by Trafford Data Lab.
In this example we want to load Local Authority District geography. The Open Geography portal helpfully provides an API explorer to help us structure the query.
In the code chunk below, we break up the query into parts that we can change depending on the type of geography we want, the areas we want to include, the fields we want the query to return etc.
# Elements for the query
geo_endpoint <- "https://ons-inspire.esriuk.com/arcgis/rest/services/"
# The geo boundary layer will change depending on which
geo_boundarylayer <- "Administrative_Boundaries/Local_Authority_Districts_December_2020_UK_BGC/"
geo_server <- "FeatureServer/0/"
geo_search <- "LAD20NM IN "
# Construct a vector of local authorities to load
# the following local authorities are the 'Solent' region
las_to_load <- c("Southampton","Portsmouth","Winchester",
"Eastleigh","Isle of Wight","Fareham",
"Gosport","Test Valley","East Hampshire",
"Havant","New Forest","Hart","Basingstoke and Deane")
geo_where <- las_to_load # sometimes we don't want all boundaries
geo_outfields <- "*" # returns all fields
#geo_outfields <- c("LAD20CD","LAD20NM","LONG","LAT") # use in place of line above to return selected fields only
geo_outSR <- "4326"
geo_format <- "json"We then paste the elements together to construct the API query URL …
# Assemble the full URL for the query from elements above
geo_query_string <- paste0(geo_endpoint,geo_boundarylayer,geo_server,
"query?where=",geo_search,"(",paste(paste0("'",geo_where,"'"), collapse = ","),
")&outFields=",(paste(geo_outfields, collapse = ",")),"&outSR=",geo_outSR,"&f=",geo_format)
# Format the URL to remove spaces
geo_query <- utils::URLencode(geo_query_string)message("Loading LA geometry from ONS Open Geography API")## Loading LA geometry from ONS Open Geography API# API query
sf_data <- st_read(geo_query)## Reading layer `ESRIJSON' from data source `https://ons-inspire.esriuk.com/arcgis/rest/services/Administrative_Boundaries/Local_Authority_Districts_December_2020_UK_BGC/FeatureServer/0/query?where=LAD20NM%20IN%20('Southampton','Portsmouth','Winchester','Eastleigh','Isle%20of%20Wight','Fareham','Gosport','Test%20Valley','East%20Hampshire','Havant','New%20Forest','Hart','Basingstoke%20and%20Deane')&outFields=*&outSR=4326&f=json' using driver `ESRIJSON'
## Simple feature collection with 13 features and 8 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -1.957224 ymin: 50.57493 xmax: -0.7446835 ymax: 51.38392
## geographic CRS: WGS 84So what did the API query return? Table 1.1 shows selected columns from the data. The data.frame contains a row for each local authority district, with the geometry column containing the geography; in this case the polygon data relating to the boundary of each local authority.
kable(sf_data[c(2:3,7:9)], caption = "Selected columns of data by API query")| lad20cd | lad20nm | long | lat | geometry |
|---|---|---|---|---|
| E06000044 | Portsmouth | -1.07006 | 50.80800 | MULTIPOLYGON (((-1.058103 5… |
| E06000045 | Southampton | -1.39952 | 50.92120 | MULTIPOLYGON (((-1.401357 5… |
| E06000046 | Isle of Wight | -1.33366 | 50.67129 | MULTIPOLYGON (((-1.289054 5… |
| E07000084 | Basingstoke and Deane | -1.22021 | 51.25937 | MULTIPOLYGON (((-1.079991 5… |
| E07000085 | East Hampshire | -0.93971 | 51.05765 | MULTIPOLYGON (((-0.8608213 … |
| E07000086 | Eastleigh | -1.32943 | 50.94872 | MULTIPOLYGON (((-1.392975 5… |
| E07000087 | Fareham | -1.23742 | 50.85388 | MULTIPOLYGON (((-1.269276 5… |
| E07000088 | Gosport | -1.16731 | 50.80636 | MULTIPOLYGON (((-1.171547 5… |
| E07000089 | Hart | -0.88321 | 51.28197 | MULTIPOLYGON (((-0.9319052 … |
| E07000090 | Havant | -1.00398 | 50.87231 | MULTIPOLYGON (((-1.019384 5… |
| E07000091 | New Forest | -1.59293 | 50.85748 | MULTIPOLYGON (((-1.819856 5… |
| E07000093 | Test Valley | -1.50214 | 51.13417 | MULTIPOLYGON (((-1.516282 5… |
| E07000094 | Winchester | -1.24393 | 51.03002 | MULTIPOLYGON (((-1.295087 5… |
Useful lookup spatial reference for CRS https://spatialreference.org/.
Sometimes transformation is required using the st_transform() function.
st_coord_sys <- st_crs(sf_data) # check coord system
st_coord_sys # current coord system EPSG: 4326 (is what leaflet wants - good)## 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["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]# transform the coord system if required
if(st_coord_sys$epsg != 4326){
sf_data <- st_transform(sf_data, "+proj=longlat +datum=WGS84")
}Now we can create a map, here using the Leaflet package.
Find a cheatsheet here.
Optional: first we can create popups by adding a column to sf_data (uses htmltools).
sf_data$popup_text <-
paste("Locial authority area code: ","<b>", sf_data$lad20cd, "</b>",
'<br/>', 'Local authority: ', '<b>', sf_data$lad20nm, '</b>', ' ') %>%
lapply(htmltools::HTML)Finally we can create the map …
leaflet(sf_data) %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addPolygons(color = "blue", fillColor = "blue", fillOpacity = 0.2, weight = 1.5, popup = ~(lad20nm), # popups clicked
label = ~(popup_text), # define labels
labelOptions = labelOptions( # label options
style = list("font-weight" = "normal", padding = "2px 2px"),
direction = "auto"),
highlight = highlightOptions(
weight = 5,
color = "#666",
fillOpacity = 0.7,
bringToFront = TRUE))