1 Refresher
Required packages
Session info
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)
Matrix products: default
locale:
[1] LC_COLLATE=English_United Kingdom.utf8
[2] LC_CTYPE=English_United Kingdom.utf8
[3] LC_MONETARY=English_United Kingdom.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=English_United Kingdom.utf8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[1] texreg_1.39.3 tidyr_1.3.1 osmdata_0.2.5
[4] nomisr_0.4.7 dplyr_1.1.4 rnaturalearth_1.0.1
[7] nngeo_0.4.8 mapview_2.11.2 gstat_2.1-1
[10] sf_1.0-16
loaded via a namespace (and not attached):
[1] xfun_0.45 raster_3.6-26 htmlwidgets_1.6.4
[4] lattice_0.22-6 vctrs_0.6.5 tools_4.4.1
[7] crosstalk_1.2.1 generics_0.1.3 stats4_4.4.1
[10] tibble_3.2.1 proxy_0.4-27 spacetime_1.3-1
[13] fansi_1.0.6 xts_0.14.0 pkgconfig_2.0.3
[16] KernSmooth_2.23-24 satellite_1.0.5 data.table_1.15.4
[19] leaflet_2.2.2 lifecycle_1.0.4 compiler_4.4.1
[22] FNN_1.1.4 rsdmx_0.6-3 munsell_0.5.1
[25] terra_1.7-78 codetools_0.2-20 snakecase_0.11.1
[28] htmltools_0.5.8.1 class_7.3-22 pillar_1.9.0
[31] classInt_0.4-10 tidyselect_1.2.1 digest_0.6.35
[34] purrr_1.0.2 fastmap_1.2.0 grid_4.4.1
[37] colorspace_2.1-0 cli_3.6.2 magrittr_2.0.3
[40] base64enc_0.1-3 XML_3.99-0.16.1 utf8_1.2.4
[43] leafem_0.2.3 e1071_1.7-14 scales_1.3.0
[46] sp_2.1-4 rmarkdown_2.27 httr_1.4.7
[49] zoo_1.8-12 png_0.1-8 evaluate_0.24.0
[52] knitr_1.47 rlang_1.1.4 Rcpp_1.0.12
[55] glue_1.7.0 DBI_1.2.3 rstudioapi_0.16.0
[58] jsonlite_1.8.8 R6_2.5.1 plyr_1.8.9
[61] intervals_0.15.4 units_0.8-5 1.1 Packages
Please make sure that you have installed the following packages:
pks <- c("dplyr",
"gstat",
"mapview",
"nngeo",
"nomisr",
"osmdata",
"rnaturalearth",
"sf",
"spatialreg",
"spdep",
"texreg",
"tidyr",
"tmap",
"viridisLite")The most important package is sf: Simple Features for R. users are strongly encouraged to install the sf binary packages from CRAN. If that does not work, please have a look at the installation instructions. It requires software packages GEOS, GDAL and PROJ.
1.2 Coordinates
In general, spatial data is structured like conventional/tidy data (e.g. data.frames, matrices), but has one additional dimension: every observation is linked to some sort of geo-spatial information. Most common types of spatial information are:
Points (one coordinate pair)
Lines (two coordinate pairs)
Polygons (at least three coordinate pairs)
Regular grids (one coordinate pair for centroid + raster / grid size)
1.2.1 Coordinate reference system (CRS)
In its raw form, a pair of coordinates consists of two numerical values. For instance, the pair c(51.752595, -1.262801) describes the location of Nuffield College in Oxford (one point). The fist number represents the latitude (north-south direction), the second number is the longitude (west-east direction), both are in decimal degrees.
However, we need to specify a reference point for latitudes and longitudes (in the Figure above: equator and Greenwich). For instance, the pair of coordinates above comes from Google Maps which returns GPS coordinates in ‘WGS 84’ (EPSG:4326).
# Coordinate pairs of two locations
coords1 <- c(51.752595, -1.262801)
coords2 <- c(51.753237, -1.253904)
coords <- rbind(coords1, coords2)
# Conventional data frame
nuffield.df <- data.frame(name = c("Nuffield College", "Radcliffe Camera"),
address = c("New Road", "Radcliffe Sq"),
lat = coords[,1], lon = coords[,2])
head(nuffield.df) name address lat lon
coords1 Nuffield College New Road 51.75259 -1.262801
coords2 Radcliffe Camera Radcliffe Sq 51.75324 -1.253904# Combine to spatial data frame
nuffield.spdf <- st_as_sf(nuffield.df,
coords = c("lon", "lat"), # Order is important
crs = 4326) # EPSG number of CRS
# Map
mapview(nuffield.spdf, zcol = "name")1.2.2 Projected CRS
However, different data providers use different CRS. For instance, spatial data in the UK usually uses ‘OSGB 1936 / British National Grid’ (EPSG:27700). Here, coordinates are in meters, and projected onto a planar 2D space.
There are a lot of different CRS projections, and different national statistics offices provide data in different projections. Data providers usually specify which reference system they use. This is important as using the correct reference system and projection is crucial for plotting and manipulating spatial data.
If you do not know the correct CRS, try starting with a standards CRS like EPSG:4326 if you have decimal degree like coordinates. If it looks like projected coordinates, try searching for the country or region in CRS libraries like https://epsg.io/. However, you must check if the projected coordinates match their real location, e.g. using mapview().
1.2.3 Why different projections?
By now, (most) people agree that the earth is not flat. So, to plot data on a 2D planar surface and to perform certain operations on a planar world, we need to make some re-projections. Depending on where we are, different re-projections of our data (globe in this case) might work better than others.
world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)[1] "sf" "data.frame"st_crs(world)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["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]# Extract a country and plot in current CRS (WGS84)
ger.spdf <- world[world$name == "Germany", ]
plot(st_geometry(ger.spdf))
# Now, let's transform Germany into a CRS optimized for Iceland
ger_rep.spdf <- st_transform(ger.spdf, crs = 5325)
plot(st_geometry(ger_rep.spdf))
Depending on the angle, a 2D projection of the earth looks different. It is important to choose a suitable projection for the available spatial data. For more information on CRS and re-projection, see e.g. Lovelace, Nowosad, and Muenchow (2019) or Stefan Jünger & Anne-Kathrin Stroppe’s GESIS workshop materials.
1.3 Importing some real world data
sf imports many of the most common spatial data files, like geojson, gpkg, or shp.
1.3.1 London shapefile (polygon)
Let’s get some administrative boundaries for London from the London Datastore. We use the sf package and its funtion st_read() to import the data.
# Create subdir (all data withh be stored in "_data")
dn <- "_data"
ifelse(dir.exists(dn), "Exists", dir.create(dn))
# Download zip file and unzip
tmpf <- tempfile()
boundary.link <- "https://data.london.gov.uk/download/statistical-gis-boundary-files-london/9ba8c833-6370-4b11-abdc-314aa020d5e0/statistical-gis-boundaries-london.zip"
download.file(boundary.link, tmpf)
unzip(zipfile = tmpf, exdir = paste0(dn))
unlink(tmpf)dn <- "_data"
# This is a shapefile
# We only need the MSOA layer for now
msoa.spdf <- st_read(dsn = paste0(dn, "/statistical-gis-boundaries-london/ESRI"),
layer = "MSOA_2011_London_gen_MHW" # Note: no file ending
)Reading layer `MSOA_2011_London_gen_MHW' from data source
`C:\work\Lehre\Geodata_Spatial_Regression\_data\statistical-gis-boundaries-london\ESRI'
using driver `ESRI Shapefile'
Simple feature collection with 983 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 503574.2 ymin: 155850.8 xmax: 561956.7 ymax: 200933.6
Projected CRS: OSGB36 / British National GridThe object msoa.spdf is our spatial data.frame. It looks essentially like a conventional data.frame, but has some additional attributes and geo-graphical information stored with it. Most importantly, notice the column geometry, which contains a list of polygons. In most cases, we have one polygon for each line / observation.
head(msoa.spdf) MSOA11CD MSOA11NM LAD11CD LAD11NM
1 E02000001 City of London 001 E09000001 City of London
2 E02000002 Barking and Dagenham 001 E09000002 Barking and Dagenham
3 E02000003 Barking and Dagenham 002 E09000002 Barking and Dagenham
4 E02000004 Barking and Dagenham 003 E09000002 Barking and Dagenham
5 E02000005 Barking and Dagenham 004 E09000002 Barking and Dagenham
6 E02000007 Barking and Dagenham 006 E09000002 Barking and Dagenham
RGN11CD RGN11NM USUALRES HHOLDRES COMESTRES POPDEN HHOLDS
1 E12000007 London 7375 7187 188 25.5 4385
2 E12000007 London 6775 6724 51 31.3 2713
3 E12000007 London 10045 10033 12 46.9 3834
4 E12000007 London 6182 5937 245 24.8 2318
5 E12000007 London 8562 8562 0 72.1 3183
6 E12000007 London 8791 8672 119 50.6 3441
AVHHOLDSZ
1 1.6
2 2.5
3 2.6
4 2.6
5 2.7
6 2.5
geometry
1 531667.6, 531647.2, 531626.7, 531667.6, 180535.0, 180532.9, 180539.0, 180535.0, 532135.1, 532158.2, 532252.9, 532293.1, 532310.1, 532419.6, 532435.8, 532445.8, 532572.3, 532634.5, 532619.1, 532746.8, 532883.5, 532946.1, 533078.8, 533195.3, 533204.7, 533230.5, 533410.7, 533352.0, 533443.0, 533460.4, 533499.1, 533604.2, 533611.0, 533743.7, 533782.6, 533803.2, 533837.9, 533842.3, 533807.9, 533649.1, 533602.4, 533551.2, 533556.3, 533483.9, 533420.1, 533403.4, 533434.1, 533395.0, 533385.4, 533382.0, 533236.7, 532916.8, 532643.1, 532474.2, 532274.6, 532279.1, 532217.0, 532016.0, 531625.6, 531145.1, 531147.4, 531074.7, 531127.9, 531112.3, 531172.2, 531133.7, 530966.7, 531158.5, 531186.6, 531275.0, 531410.7, 531521.6, 531554.9, 531627.8, 531796.4, 532022.4, 532082.9, 532105.1, 532089.3, 532071.3, 532135.1, 182198.1, 182151.1, 182179.3, 182068.4, 182003.6, 181998.3, 182003.3, 181957.9, 181938.9, 181926.0, 181847.2, 181786.9, 181715.4, 181894.8, 181840.5, 181970.9, 182000.0, 182082.5, 182037.9, 181739.7, 181755.5, 181669.5, 181582.4, 181418.1, 181411.0, 181261.3, 181161.6, 181053.4, 180857.9, 180789.3, 180767.8, 180733.9, 180815.2, 180812.8, 180778.0, 180758.2, 180691.8, 180649.9, 180623.2, 180515.2, 180516.5, 180510.7, 180573.3, 180594.6, 180624.6, 180677.0, 180742.3, 180789.8, 180749.4, 180800.8, 180798.0, 180782.1, 180882.2, 180990.3, 181049.5, 181104.3, 181124.6, 181279.0, 181596.7, 181603.5, 181599.9, 181581.1, 181576.1, 181639.8, 181659.3, 181695.7, 181784.5, 181893.5, 181911.8, 182011.2, 182077.0, 182159.6, 182198.1
2 548881.6, 548881.1, 549000.7, 549089.4, 549084.2, 549071.8, 549116.5, 549121.0, 548464.4, 548481.9, 548466.0, 548462.0, 548452.0, 548443.9, 548344.0, 548357.8, 548397.1, 548394.1, 548234.9, 548114.1, 548079.0, 548044.6, 548044.0, 548058.5, 547960.3, 547804.4, 547778.3, 547735.4, 547770.7, 547921.6, 547882.6, 547856.8, 547813.7, 547514.6, 547515.1, 547533.1, 547536.9, 547510.1, 547482.1, 547433.0, 547412.2, 547519.7, 547729.3, 547618.2, 547865.1, 547973.4, 548453.8, 548767.4, 548877.7, 548845.7, 548881.6, 190845.3, 190802.5, 190873.5, 190714.1, 190584.4, 190315.8, 190074.8, 189772.3, 189672.7, 189556.8, 189263.0, 189132.0, 189024.7, 188953.8, 188916.1, 188845.1, 188848.7, 188776.6, 188749.5, 188728.0, 188930.0, 188924.2, 188793.3, 188751.7, 188743.2, 188626.8, 188711.1, 188840.8, 188886.7, 188943.2, 189049.4, 189093.7, 189096.5, 189119.7, 189127.6, 189227.1, 189301.0, 189373.9, 189374.1, 189612.1, 189680.7, 189711.4, 189779.2, 190059.9, 190258.1, 190348.2, 190810.1, 191139.0, 191097.3, 190846.8, 190845.3
3 549102.4, 548954.5, 548889.7, 548874.2, 548952.5, 548954.4, 548946.4, 549019.1, 549051.7, 549085.6, 548972.7, 548871.1, 548880.3, 548881.2, 548919.3, 549019.5, 549057.9, 549099.6, 549057.4, 548997.8, 549009.9, 548956.7, 548911.0, 548895.8, 548836.2, 548849.4, 548884.3, 548798.8, 548754.9, 548587.0, 548651.2, 548410.4, 548307.5, 548205.4, 548113.3, 547828.7, 547706.5, 547682.9, 547704.3, 547699.6, 547700.9, 547721.5, 547746.5, 547715.1, 547700.3, 547764.8, 547817.9, 547767.6, 547701.8, 547662.6, 547656.6, 547655.1, 547642.4, 547697.0, 547804.4, 547960.3, 548058.5, 548044.0, 548044.6, 548079.0, 548114.1, 548234.9, 548394.1, 548397.1, 548357.8, 548344.0, 548443.9, 548452.0, 548462.0, 548466.0, 548481.9, 548464.4, 549121.0, 549120.8, 549102.4, 189324.6, 189063.2, 189065.7, 188496.2, 188492.3, 188246.3, 188214.2, 188228.3, 187899.8, 187836.7, 187816.6, 187799.0, 187766.3, 187727.4, 187717.8, 187715.0, 187664.8, 187656.1, 187578.5, 187577.1, 187528.1, 187476.1, 187463.3, 187495.6, 187479.0, 187412.7, 187278.5, 187271.9, 187368.4, 187432.0, 187581.6, 187606.0, 187596.6, 187517.0, 187469.0, 187450.8, 187444.4, 187443.2, 187598.1, 187618.6, 187707.7, 187759.0, 187885.7, 187914.5, 188007.7, 188024.1, 188037.6, 188214.7, 188170.1, 188271.9, 188376.2, 188431.4, 188528.7, 188556.9, 188626.8, 188743.2, 188751.7, 188793.3, 188924.2, 188930.0, 188728.0, 188749.5, 188776.6, 188848.7, 188845.1, 188916.1, 188953.8, 189024.7, 189132.0, 189263.0, 189556.8, 189672.7, 189772.3, 189640.2, 189324.6
4 551550.0, 551478.0, 551441.7, 551397.2, 551528.4, 551543.0, 551526.7, 551572.3, 551554.4, 551619.4, 551740.0, 551779.9, 551943.8, 551905.8, 551765.0, 551670.8, 551204.5, 551197.9, 550485.1, 550425.2, 550335.8, 550247.1, 550200.7, 550239.0, 550531.5, 550406.2, 550272.8, 550244.1, 550171.0, 550134.7, 550106.5, 550121.0, 550076.0, 550100.0, 550064.6, 550083.7, 550141.2, 550140.0, 550167.1, 550299.4, 550295.0, 550324.1, 550350.4, 550380.5, 550414.1, 550429.2, 550413.7, 550411.0, 550392.6, 550500.4, 550519.0, 550739.3, 550978.8, 551129.7, 551430.6, 551551.6, 551550.0, 187364.6, 187263.7, 187137.0, 186966.9, 186852.2, 186762.3, 186570.8, 186477.1, 186432.7, 186429.7, 186318.5, 186178.4, 186027.7, 185894.1, 185746.5, 185446.8, 185301.9, 185502.7, 185441.1, 185385.0, 185516.6, 185658.5, 185624.8, 185744.4, 185756.4, 185898.1, 185953.1, 185994.0, 185971.0, 185957.3, 185984.5, 186019.8, 186039.5, 186065.0, 186096.4, 186125.6, 186101.6, 186225.4, 186262.4, 186255.3, 186375.0, 186374.3, 186269.7, 186430.0, 186481.0, 186734.1, 186968.6, 187021.3, 187189.0, 187258.3, 187270.0, 187373.8, 187347.3, 187363.3, 187345.3, 187445.2, 187364.6
5 549099.6, 549161.4, 549210.3, 549418.6, 549145.8, 549138.7, 549135.0, 549128.7, 548983.4, 548970.8, 548933.1, 548958.3, 548992.2, 549027.3, 549074.5, 549156.0, 549153.0, 549153.9, 549077.1, 549026.8, 548950.0, 548933.9, 548829.6, 548787.2, 548741.3, 548693.7, 548675.6, 548724.9, 548708.2, 548697.8, 548802.6, 548816.9, 548637.0, 548619.9, 548479.9, 548477.4, 548423.0, 548219.5, 548136.3, 548171.4, 548152.6, 548058.0, 548080.9, 548046.4, 548065.5, 548166.7, 548172.6, 548183.8, 548179.2, 548201.9, 548262.1, 548205.4, 548307.5, 548410.4, 548651.2, 548587.0, 548754.9, 548798.8, 548884.3, 548849.4, 548836.2, 548895.8, 548911.0, 548956.7, 549009.9, 548997.8, 549057.4, 549099.6, 187656.1, 187625.6, 187648.6, 187442.5, 187382.9, 187195.2, 186976.1, 186912.5, 186875.1, 186589.7, 186439.2, 186426.8, 186459.4, 186392.2, 186436.5, 186429.5, 186380.6, 186255.1, 186262.7, 186159.8, 186174.0, 186160.1, 186002.1, 186018.2, 185971.8, 185979.0, 186019.8, 186041.2, 186124.9, 186138.5, 186261.4, 186283.4, 186272.0, 186270.9, 186263.0, 186320.0, 186385.5, 186424.0, 186451.3, 186545.7, 186621.7, 186664.5, 186752.6, 186780.8, 186815.5, 186868.7, 186879.0, 187040.4, 187081.8, 187084.3, 187090.8, 187517.0, 187596.6, 187606.0, 187581.6, 187432.0, 187368.4, 187271.9, 187278.5, 187412.7, 187479.0, 187495.6, 187463.3, 187476.1, 187528.1, 187577.1, 187578.5, 187656.1
6 549819.9, 549950.8, 550250.9, 550392.6, 550411.0, 550413.7, 550429.2, 550414.1, 550380.5, 550350.4, 550324.1, 550295.0, 550299.4, 550167.1, 550140.0, 550141.2, 550083.7, 550064.6, 550100.0, 550076.0, 550121.0, 550106.5, 549955.3, 549892.7, 549899.1, 549912.6, 549604.0, 549341.7, 549285.5, 549140.9, 549100.9, 549023.0, 549001.5, 549020.4, 549089.7, 549143.7, 549039.9, 549052.3, 549026.8, 549077.1, 549153.9, 549153.0, 549156.0, 549074.5, 549027.3, 548992.2, 548958.3, 548933.1, 548970.8, 548983.4, 549128.7, 549135.0, 549138.7, 549145.8, 549418.6, 549547.2, 549819.9, 187080.7, 186961.1, 187148.6, 187189.0, 187021.3, 186968.6, 186734.1, 186481.0, 186430.0, 186269.7, 186374.3, 186375.0, 186255.3, 186262.4, 186225.4, 186101.7, 186125.6, 186096.4, 186065.0, 186039.5, 186019.8, 185984.5, 185936.2, 185973.8, 185885.6, 185841.2, 185773.2, 185800.9, 185796.6, 185788.5, 185827.9, 185817.6, 185872.2, 185931.3, 185924.6, 185959.8, 186065.6, 186131.9, 186159.8, 186262.7, 186255.1, 186380.6, 186429.5, 186436.5, 186392.2, 186459.4, 186426.8, 186439.2, 186589.7, 186875.1, 186912.5, 186976.1, 187195.2, 187382.9, 187442.5, 187336.8, 187080.7Shapefiles are still among the most common formats to store and transmit spatial data, despite them being inefficient (file size and file number).
However, sf reads everything spatial, such as geo.json, which usually is more efficient, but less common (but we’re getting there).
# Download file
ulez.link <- "https://data.london.gov.uk/download/ultra_low_emissions_zone/936d71d8-c5fc-40ad-a392-6bec86413b48/CentralUltraLowEmissionZone.geojson"
download.file(ulez.link, paste0(dn, "/ulez.json"))# Read geo.json
st_layers(paste0(dn, "/ulez.json"))Driver: GeoJSON
Available layers:
layer_name geometry_type features fields
1 CentralUltraLowEmissionZone Multi Polygon 1 4
crs_name
1 OSGB36 / British National Gridulez.spdf <- st_read(dsn = paste0(dn, "/ulez.json")) # here dsn is simply the fileReading layer `CentralUltraLowEmissionZone' from data source
`C:\work\Lehre\Geodata_Spatial_Regression\_data\ulez.json'
using driver `GeoJSON'
Simple feature collection with 1 feature and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 527271.5 ymin: 178041.5 xmax: 533866.3 ymax: 183133.4
Projected CRS: OSGB36 / British National Gridhead(ulez.spdf)Simple feature collection with 1 feature and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 527271.5 ymin: 178041.5 xmax: 533866.3 ymax: 183133.4
Projected CRS: OSGB36 / British National Grid
fid OBJECTID BOUNDARY Shape_Area geometry
1 1 1 CSS Area 21.37557 MULTIPOLYGON (((531562.7 18...Again, this looks like a conventional data.frame but has the additional column geometry containing the coordinates of each observation. st_geometry() returns only the geographic object and st_drop_geometry() only the data.frame without the coordinates. We can plot the object using mapview().
mapview(msoa.spdf[, "POPDEN"])1.3.2 Census API (admin units)
Now that we have some boundaries and shapes of spatial units in London, we can start looking for different data sources to populate the geometries.
A good source for demographic data is for instance the 2011 census. Below we use the nomis API to retrieve population data for London, See the Vignette for more information (Guest users are limited to 25,000 rows per query). Below is a wrapper to avoid some errors with sex and urban-rural cross-tabulation in some of the data.
### For larger request, register and set key
# Sys.setenv(NOMIS_API_KEY = "XXX")
# nomis_api_key(check_env = TRUE)
x <- nomis_data_info()
# Get London ids
london_ids <- msoa.spdf$MSOA11CD
### Get key statistics ids
# select requires tables (https://www.nomisweb.co.uk/sources/census_2011_ks)
# Let's get KS201EW (ethnic group), KS205EW (passport held), and KS402EW (housing tenure)
# Get internal ids
stats <- c("KS201EW", "KS402EW", "KS205EW")
oo <- which(grepl(paste(stats, collapse = "|"), x$name.value))
ksids <- x$id[oo]
ksids # This are the internal ids
### look at meta information
q <- nomis_overview(ksids[1])
head(q)
a <- nomis_get_metadata(id = ksids[1], concept = "GEOGRAPHY", type = "type")
a # TYPE297 is MSOA level
b <- nomis_get_metadata(id = ksids[1], concept = "MEASURES", type = "TYPE297")
b # 20100 is the measure of absolute numbers
### Query data in loop over the required statistics
for(i in ksids){
# Determin if data is divided by sex or urban-rural
nd <- nomis_get_metadata(id = i)
if("RURAL_URBAN" %in% nd$conceptref){
UR <- TRUE
}else{
UR <- FALSE
}
if("C_SEX" %in% nd$conceptref){
SEX <- TRUE
}else{
SEX <- FALSE
}
# make data request
if(UR == TRUE){
if(SEX == TRUE){
tmp_en <- nomis_get_data(id = i, time = "2011",
geography = london_ids, # replace with "TYPE297" for all MSOAs
measures = 20100, RURAL_URBAN = 0, C_SEX = 0)
}else{
tmp_en <- nomis_get_data(id = i, time = "2011",
geography = london_ids, # replace with "TYPE297" for all MSOAs
measures = 20100, RURAL_URBAN = 0)
}
}else{
if(SEX == TRUE){
tmp_en <- nomis_get_data(id = i, time = "2011",
geography = london_ids, # replace with "TYPE297" for all MSOAs
measures = 20100, C_SEX = 0)
}else{
tmp_en <- nomis_get_data(id = i, time = "2011",
geography = london_ids, # replace with "TYPE297" for all MSOAs
measures = 20100)
}
}
# Append (in case of different regions)
ks_tmp <- tmp_en
# Make lower case names
names(ks_tmp) <- tolower(names(ks_tmp))
names(ks_tmp)[names(ks_tmp) == "geography_code"] <- "msoa11"
names(ks_tmp)[names(ks_tmp) == "geography_name"] <- "name"
# replace weird cell codes
onlynum <- which(grepl("^[[:digit:]]+$", ks_tmp$cell_code))
if(length(onlynum) != 0){
code <- substr(ks_tmp$cell_code[-onlynum][1], 1, 7)
if(is.na(code)){
code <- i
}
ks_tmp$cell_code[onlynum] <- paste0(code, "_", ks_tmp$cell_code[onlynum])
}
# save codebook
ks_cb <- unique(ks_tmp[, c("date", "cell_type", "cell", "cell_code", "cell_name")])
### Reshape
ks_res <- tidyr::pivot_wider(ks_tmp, id_cols = c("msoa11", "name"),
names_from = "cell_code",
values_from = "obs_value")
### Merge
if(i == ksids[1]){
census_keystat.df <- ks_res
census_keystat_cb.df <- ks_cb
}else{
census_keystat.df <- merge(census_keystat.df, ks_res, by = c("msoa11", "name"), all = TRUE)
census_keystat_cb.df <- rbind(census_keystat_cb.df, ks_cb)
}
}
# Descriptions are saved in the codebook
save(census_keystat.df, file = "_data/Census_ckeystat.RData")
save(census_keystat_cb.df, file = "_data/Census_codebook.RData")Now, we have one file containing the geometries of MSOAs and one file with the census information on ethnic groups. Obviously, we can easily merge them together using the MSOA identifiers.
And we can, for instance, plot the spatial distribution of ethnic groups.
msoa.spdf$per_white <- msoa.spdf$KS201EW_100 / msoa.spdf$KS201EW0001 * 100
msoa.spdf$per_mixed <- msoa.spdf$KS201EW_200 / msoa.spdf$KS201EW0001 * 100
msoa.spdf$per_asian <- msoa.spdf$KS201EW_300 / msoa.spdf$KS201EW0001 * 100
msoa.spdf$per_black <- msoa.spdf$KS201EW_400 / msoa.spdf$KS201EW0001 * 100
msoa.spdf$per_other <- msoa.spdf$KS201EW_500 / msoa.spdf$KS201EW0001 * 100
mapview(msoa.spdf[, "per_white"])If you’re interested in more data sources, see for instance APIs for social scientists: A collaborative review by Paul C. Bauer, Camille Landesvatter, Lion Behrens. It’s a collection of several APIs for social sciences.
1.3.3 Gridded data
So far, we have queried data on administrative units. However, often data comes on other spatial scales. For instance, we might be interested in the amount of air pollution, which is provided on a regular grid across the UK from Defra.
# Download
pol.link <- "https://uk-air.defra.gov.uk/datastore/pcm/mapno22011.csv"
download.file(pol.link, paste0(dn, "/mapno22011.csv"))pol.df <- read.csv(paste0(dn, "/mapno22011.csv"), skip = 5, header = T, sep = ",",
stringsAsFactors = F, na.strings = "MISSING")
head(pol.df) ukgridcode x y no22011
1 54291 460500 1221500 NA
2 54292 461500 1221500 NA
3 54294 463500 1221500 NA
4 54979 458500 1220500 NA
5 54980 459500 1220500 NA
6 54981 460500 1220500 NAThe data comes as point data with x and y as coordinates. We have to transform this into spatial data first. We first setup a spatial points object with st_as_sf. Subsequently, we transform the point coordinates into a regular grid. We use a buffer method st_buffer with “diameter”, and only one segment per quadrant (nQuadSegs). This gives us a 1x1km regular grid.
1.3.4 OpenStreetMap (points)
Another interesting data source is the OpenStreetMap API, which provides information about the geographical location of a serious of different indicators. Robin Lovelace provides a nice introduction to the osmdata API. Available features can be found on OSM wiki.
First we create a bounding box of where we want to query data. st_bbox() can be used to get bounding boxes of an existing spatial object (needs CRS = 4326). An alternative would be to use opq(bbox = 'greater london uk').
# bounding box of where we want to query data
q <- opq(bbox = st_bbox(st_transform(msoa.spdf, 4326)))And we want to get data for all pubs and bars which are within this bounding box.
# First build the query of location of pubs in London
osmq <- add_osm_feature(q, key = "amenity", value = "pub")
# And then query the data
pubs.osm <- osmdata_sf(osmq)Right now there are some results in polygons, some in points, and they overlap. Often, data from OSM needs some manual cleaning. Sometimes the same features are represented by different spatial objects (e.g. points + polygons).
# Make unique points / polygons
pubs.osm <- unique_osmdata(pubs.osm)
# Get points and polygons (there are barley any pubs as polygons, so we ignore them)
pubs.points <- pubs.osm$osm_points
pubs.polys <- pubs.osm$osm_multipolygons
# # Drop OSM file
# rm(pubs.osm); gc()
# Reduce to point object only
pubs.spdf <- pubs.points
# Reduce to a few variables
pubs.spdf <- pubs.spdf[, c("osm_id", "name", "addr:postcode", "diet:vegan")]Again, we can inspect the results with mapview.
# Reload (as I don't run the above here)
load("_data/osm_d.RData")mapview(st_geometry(pubs.spdf))Note that OSM is solely based on contribution by users, and the quality of OSM data varies. Usually data quality is better in larger cities, and better for more stable features (such as hospitals, train stations, highways) rahter than pubs or restaurants which regularly appear and disappear. However, data from London Datastore would indicate more pubs than what we find with OSM.
1.3.5 Save
We will store the created data to use them again in the next session.