3 Exercises I
\[ \newcommand{\Exp}{\mathrm{E}} \newcommand\given[1][]{\:#1\vert\:} \newcommand{\Cov}{\mathrm{Cov}} \newcommand{\Var}{\mathrm{Var}} \newcommand{\rank}{\mathrm{rank}} \newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}} \]
Required packages
Session info
R version 4.5.1 (2025-06-13 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)
Matrix products: default
LAPACK version 3.12.1
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] gridExtra_2.3 ggthemes_5.1.0 ggplot2_3.5.2
[4] viridisLite_0.4.2 tmap_4.1 spatialreg_1.3-6
[7] Matrix_1.7-3 spdep_1.3-13 spData_2.3.4
[10] areal_0.1.8 mapview_2.11.2 sf_1.0-21
loaded via a namespace (and not attached):
[1] tidyselect_1.2.1 dplyr_1.1.4
[3] farver_2.1.2 fastmap_1.2.0
[5] leaflegend_1.2.1 leaflet_2.2.2
[7] TH.data_1.1-3 XML_3.99-0.18
[9] digest_0.6.37 lifecycle_1.0.4
[11] LearnBayes_2.15.1 survival_3.8-3
[13] terra_1.8-54 magrittr_2.0.3
[15] compiler_4.5.1 rlang_1.1.6
[17] tools_4.5.1 data.table_1.17.6
[19] knitr_1.50 htmlwidgets_1.6.4
[21] sp_2.2-0 classInt_0.4-11
[23] RColorBrewer_1.1-3 multcomp_1.4-28
[25] abind_1.4-8 KernSmooth_2.23-26
[27] purrr_1.0.4 withr_3.0.2
[29] leafsync_0.1.0 grid_4.5.1
[31] stats4_4.5.1 cols4all_0.8
[33] e1071_1.7-16 leafem_0.2.4
[35] colorspace_2.1-1 spacesXYZ_1.6-0
[37] scales_1.4.0 MASS_7.3-65
[39] dichromat_2.0-0.1 cli_3.6.5
[41] mvtnorm_1.3-3 rmarkdown_2.29
[43] generics_0.1.4 rstudioapi_0.17.1
[45] tmaptools_3.2 DBI_1.2.3
[47] proxy_0.4-27 stringr_1.5.1
[49] splines_4.5.1 stars_0.6-8
[51] parallel_4.5.1 s2_1.1.9
[53] base64enc_0.1-3 vctrs_0.6.5
[55] boot_1.3-31 sandwich_3.1-1
[57] jsonlite_2.0.0 crosstalk_1.2.1
[59] units_0.8-7 maptiles_0.10.0
[61] glue_1.8.0 lwgeom_0.2-14
[63] codetools_0.2-20 leaflet.providers_2.0.0
[65] stringi_1.8.7 gtable_0.3.6
[67] deldir_2.0-4 raster_3.6-32
[69] tibble_3.3.0 logger_0.4.0
[71] pillar_1.10.2 htmltools_0.5.8.1
[73] satellite_1.0.5 R6_2.6.1
[75] wk_0.9.4 microbenchmark_1.5.0
[77] evaluate_1.0.4 lattice_0.22-7
[79] png_0.1-8 class_7.3-23
[81] Rcpp_1.0.14 coda_0.19-4.1
[83] nlme_3.1-168 xfun_0.52
[85] zoo_1.8-14 pkgconfig_2.0.3 Reload data from pervious session
load("_data/msoa2_spatial.RData")3.1 General Exercises
3.1.1 1) Can you import the spatial administrative units of Germany (“Kreisgrenzen_2020_mit_Einwohnerzahl” in _data folder) and make a simple plot of the boundaries? {.unnumbered}
# Import shape file layer
ger.sdpf <- st_read(dsn = "_data/Kreisgrenzen_2020_mit_Einwohnerzahl",
layer = "KRS_ew_20")Reading layer `KRS_ew_20' from data source
`C:\work\Lehre\Geodata_Spatial_Regression\_data\Kreisgrenzen_2020_mit_Einwohnerzahl'
using driver `ESRI Shapefile'
Simple feature collection with 401 features and 19 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 5.86625 ymin: 47.27012 xmax: 15.04182 ymax: 55.05878
Geodetic CRS: WGS 84# Plot via ggplot
gp <- ggplot(ger.sdpf)+
geom_sf( color = "magenta", fill = NA)+
coord_sf(datum = NA)+
theme_map()
gp
2) What is the Coordinate reference system of this German shape file?
st_crs(ger.sdpf)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]]3) Using the MSOA data in London from above: Please draw a map showing the distribution of share of non-EU immigrants (per_nonEU).
gp1 <- ggplot(msoa.spdf)+
geom_sf(aes(fill = per_nonEU))+
scale_fill_viridis_c(option = "A")+
coord_sf(datum = NA)+
theme_map()+
theme(legend.position = c(.9, .4))Warning: A numeric `legend.position` argument in `theme()` was deprecated in
ggplot2 3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()`
instead.gp1
4) Below we will load some additional data on heat islands.
We will add some Data about London’s Urban Heat Island. It contains information about the mean temperature at midnight during the summer of 2011.
This is a tif file that we need to read in with stars and then transform into sf.
Loading required package: abind# Read geo_tif with stars
urbclim <- read_stars("https://data.london.gov.uk/download/ae16d5af-5dce-49bc-b1e2-88bb41e8bfd0/f4e3a05d-fad7-4b56-8c42-ba2274f3bb3a/London_Tmin_midnight_2011.tif")
# urbclim <- read_stars("_data/London_Tmin_midnight_2011.tif")
# Transfer to sf
urbclim.spdf <- st_as_sf(urbclim)
names(urbclim.spdf)[1] <- "Tmin_midnight"On which projection is the urbclim temperature data?
Can you please calculate the average night-time temperature for each MSOA using area weighted interpolation. Make sure that the objects are on the same projections / crs.
Create a map showing the temperature for each MSOA, and plot it next to the maps of non-EU immigrant residents (e.g. using
grid.arrange()) .
# Check projection
st_crs(urbclim.spdf)Coordinate Reference System:
User input: unnamed
wkt:
PROJCRS["unnamed",
BASEGEOGCRS["Unknown datum based upon the GRS 1980 ellipsoid",
DATUM["Not specified (based on GRS 1980 ellipsoid)",
ELLIPSOID["GRS 1980",6378137,298.257222101004,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4019]],
CONVERSION["Lambert Azimuthal Equal Area",
METHOD["Lambert Azimuthal Equal Area",
ID["EPSG",9820]],
PARAMETER["Latitude of natural origin",52,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",10,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",4321000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",3210000,
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]]]# Add id
urbclim.spdf$id <- rownames(urbclim.spdf)
# Bring on common crs
urbclim.spdf <- st_transform(urbclim.spdf, crs = st_crs(msoa.spdf))
# Use area weights interpolation to merge
msoa.spdf <- aw_interpolate(
msoa.spdf,
tid = "MSOA11CD",
source = urbclim.spdf,
sid = "id",
weight = "sum",
output = "sf",
intensive = "Tmin_midnight"
)# Make map of Temperature
gp2 <- ggplot(msoa.spdf)+
geom_sf(aes(fill = Tmin_midnight))+
scale_fill_viridis_c(option = "C")+
coord_sf(datum = NA)+
theme_map()+
theme(legend.position = c(.9, .4))
# Plot two ggplot maps next to each other
gp <- grid.arrange(gp1, gp2, nrow = 1)
gpTableGrob (1 x 2) "arrange": 2 grobs
z cells name grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (1-1,2-2) arrange gtable[layout]