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.6.0 (2026-04-24 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)
Matrix products: default
LAPACK version 3.12.1
locale:
[1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=German_Germany.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.2.0 ggplot2_4.0.3
[4] viridisLite_0.4.3 tmap_4.4 spatialreg_1.4-3
[7] Matrix_1.7-5 spdep_1.4-2 spData_2.3.5
[10] areal_0.1.8 mapview_2.11.4 sf_1.1-1
loaded via a namespace (and not attached):
[1] tidyselect_1.2.1 dplyr_1.2.1
[3] farver_2.1.2 S7_0.2.2
[5] fastmap_1.2.0 leaflegend_1.2.8
[7] leaflet_2.2.3 TH.data_1.1-5
[9] XML_3.99-0.23 digest_0.6.39
[11] lifecycle_1.0.5 LearnBayes_2.15.2
[13] survival_3.8-6 terra_1.9-27
[15] magrittr_2.0.5 compiler_4.6.0
[17] rlang_1.2.0 tools_4.6.0
[19] data.table_1.18.4 knitr_1.51
[21] htmlwidgets_1.6.4 sp_2.2-1
[23] classInt_0.4-11 RColorBrewer_1.1-3
[25] multcomp_1.4-30 abind_1.4-8
[27] KernSmooth_2.23-26 purrr_1.2.2
[29] withr_3.0.2 leafsync_0.1.0
[31] grid_4.6.0 stats4_4.6.0
[33] cols4all_0.10 e1071_1.7-17
[35] leafem_0.2.5 colorspace_2.1-2
[37] spacesXYZ_1.6-0 scales_1.4.0
[39] MASS_7.3-65 cli_3.6.6
[41] mvtnorm_1.3-7 rmarkdown_2.31
[43] generics_0.1.4 otel_0.2.0
[45] rstudioapi_0.18.0 tmaptools_3.3
[47] DBI_1.3.0 proxy_0.4-29
[49] stringr_1.6.0 splines_4.6.0
[51] stars_0.7-2 parallel_4.6.0
[53] s2_1.1.11 base64enc_0.1-6
[55] marginaleffects_0.32.0 vctrs_0.7.3
[57] boot_1.3-32 sandwich_3.1-1
[59] jsonlite_2.0.0 crosstalk_1.2.2
[61] maptiles_0.11.0 units_1.0-1
[63] glue_1.8.1 lwgeom_0.2-16
[65] codetools_0.2-20 stringi_1.8.7
[67] gtable_0.3.6 deldir_2.0-4
[69] raster_3.6-32 tibble_3.3.1
[71] logger_0.4.2 pillar_1.11.1
[73] htmltools_0.5.9 satellite_1.0.6
[75] R6_2.6.1 wk_0.9.5
[77] microbenchmark_1.5.0 evaluate_1.0.5
[79] lattice_0.22-9 png_0.1-9
[81] backports_1.5.1 class_7.3-23
[83] Rcpp_1.1.1-1.1 coda_0.19-4.1
[85] nlme_3.1-169 xfun_0.57
[87] zoo_1.8-15 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” from the _data folder) and make a simple map of the geographies / adminstrative 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? Please print the CRS.
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 the share of non-EU immigrants (the variable is `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))
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.
Lade nötiges Paket: 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]