1

Import the ESRI shapefile of German districts and the district attribute table. Join the two data frames, transform the CRS to “EPSG:3035” and check your changes.

You need to rename one of the id variables or adjust your join accordingly (“id = district_id”).

solution 

# load libraries
library(sf)
## Warning: package 'sf' was built under R version 4.0.3
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# Import data
german_districts <- st_read(dsn = "../data",
                           layer = "GER_DISTRICTS") %>% 
                     rename(., district_id = id) #
## Reading layer `GER_DISTRICTS' from data source `C:\Users\annes\Documents\gesis-workshop-geospatial-techniques-R\data' using driver `ESRI Shapefile'
## Simple feature collection with 401 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 5.86625 ymin: 47.27012 xmax: 15.04182 ymax: 55.05838
## geographic CRS: ETRS89
attributes_districts <-  read.csv("../data/attributes_districts.csv", 
                                  header = T, fill = T, sep = ",") 


# Join data and transform
german_districts_enhanced <- 
  german_districts %>% 
  left_join(., attributes_districts, by = "district_id") %>% 
  st_transform(., crs = 3035)

# Check
st_crs(german_districts_enhanced)
## Coordinate Reference System:
##   User input: EPSG:3035 
##   wkt:
## PROJCRS["ETRS89-extended / LAEA Europe",
##     BASEGEOGCRS["ETRS89",
##         DATUM["European Terrestrial Reference System 1989",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["Europe Equal Area 2001",
##         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["northing (Y)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (X)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Europe - LCC & LAEA"],
##         BBOX[24.6,-35.58,84.17,44.83]],
##     ID["EPSG",3035]]
head(german_districts_enhanced, 2)
## Simple feature collection with 2 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 4279627 ymin: 3460480 xmax: 4335232 ymax: 3524426
## projected CRS:  ETRS89-extended / LAEA Europe
##   district_id population cases deaths cases_per_100k cases_7days death_7days
## 1        1001      90164  1044     21       1157.890         108           1
## 2        1002     246794  3037     71       1230.581         130           2
##   afd_voteshare_2017                       geometry
## 1           7.449591 MULTIPOLYGON (((4283235 352...
## 2           6.877478 MULTIPOLYGON (((4331981 348...

2

We want a first descriptive visual of the distribution of Covid-19 cases in Cologne and the surrounding districts. Calculate the number of Covid-19 cases in the last seven days (cases_7days) by population (population) and multiply with 100,000 (in Germany usually called “7 Tages Inzidenzzahl”: number of persons who were infected with Covid-19 in the last seven days per 100k).

Select Cologne (district_id == 5315), find the surrounding districts, and plot Cologne and its surrounding districts.

You can use the dplyr function bind_rows to combine the two spatial objects, “Cologne” and “Cologne Surroundings”.

solution 

# calculate Covid-19 rate
german_districts_enhanced <-
  german_districts_enhanced %>% 
  mutate(covid7d_rate = (cases_7days / population) * 100000)

# filter Cologne
cologne <-
  german_districts_enhanced %>% 
  filter(. , district_id == 5315)

# filter surrounding districts, append with Cologne data and select the Covid column
cologne_sur <-
  german_districts_enhanced %>%
  filter(lengths(st_touches(., cologne)) > 0) %>% 
  bind_rows(., cologne) %>%   
  select(. , covid7d_rate)

# plot  
plot(cologne_sur)

3

Save your data set of Cologne and its surrounding districts as an ESRI Shapefile.

solution 

# Export as shapefile
st_write(cologne_sur, 
         dsn = "./data/own_material/cologne_covid19_epsg3035", 
         delete_layer = TRUE) #optional