1

Import the corona_cologne dataset into your R session.

There are several corona_cologne files in the ../data/ folder, you would only need the one with the .shp extension.

solution 

corona_cologne <-
  sf::read_sf(
    "./data/corona_cologne.shp"
  )

2

Use the same command as before, but instead of defining the path to the file location plug in the following string:

"https://geoportal.stadt-koeln.de/arcgis/rest/services/Politik_und_Verwaltung/covid_stadtteile/MapServer/1/query?where=id+is+not+null&text=&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&distance=&units=esriSRUnit_Foot&relationParam=&outFields=*&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=4326&havingClause=&returnIdsOnly=false&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&historicMoment=&returnDistinctValues=false&resultOffset=&resultRecordCount=&returnExtentOnly=false&datumTransformation=&parameterValues=&rangeValues=&quantizationParameters=&featureEncoding=esriDefault&f=geojson"

Have a look at the results. What do you think has happened?

If you don’t like long strings occupying a lot of horizontal screen space, consider cutting it in pieces. I prefer using the paste()/paste0() function or functions from additional packages such as glue::glue() from the glue package.

solution 

library(dplyr)

corona_cologne <- 
  glue::glue(
    "https://geoportal.stadt-koeln.de/arcgis/rest/services/\\
    Politik_und_Verwaltung/covid_stadtteile/MapServer/1/query?\\
    where=id+is+not+null&text=&objectIds=&time=&geometry=&\\
    geometryType=esriGeometryEnvelope&inSR=&spatialRel=\\
    esriSpatialRelIntersects&distance=&units=esriSRUnit_Foot&relationParam=\\
    &outFields=*&returnGeometry=true&returnTrueCurves=false&\\
    maxAllowableOffset=&geometryPrecision=&outSR=4326&havingClause=\\
    &returnIdsOnly=false&returnCountOnly=false&orderByFields=\\
    &groupByFieldsForStatistics=&outStatistics=&returnZ=false\\
    &returnM=false&gdbVersion=&historicMoment=&returnDistinctValues=false\\
    &resultOffset=&resultRecordCount=&returnExtentOnly=false\\
    &datumTransformation=&parameterValues=&rangeValues=\\
    &quantizationParameters=&featureEncoding=esriDefault&f=geojson"
  ) %>% 
  sf::st_read()
## Reading layer `OGRGeoJSON' from data source 
##   `https://geoportal.stadt-koeln.de/arcgis/rest/services/Politik_und_Verwaltung/covid_stadtteile/MapServer/1/query?where=id+is+not+null&text=&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&distance=&units=esriSRUnit_Foot&relationParam=&outFields=*&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=4326&havingClause=&returnIdsOnly=false&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&historicMoment=&returnDistinctValues=false&resultOffset=&resultRecordCount=&returnExtentOnly=false&datumTransformation=&parameterValues=&rangeValues=&quantizationParameters=&featureEncoding=esriDefault&f=geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 86 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 6.77253 ymin: 50.83045 xmax: 7.162028 ymax: 51.08496
## Geodetic CRS:  WGS 84
# This exercise shows the beauty of the `sf::st_read` function. It not only
# detects several different geospatial data formats; it also can download data
# directly from the internet (in this case a .geojson data file). Downloading 
# data directly is something that comes quite handy in the geospatial data world 
# as a lot of data are distributed in the form of download links or even APIs.

3

Transform the CRS of the dataset to EPSG:3035 using the sf::st_transform() function.

You only need the four digit number for defining the CRS.

solution 

corona_cologne <-
  corona_cologne %>% 
  sf::st_transform(3035)

4

Load both the immigrants_cologne and the inhabitants_cologne raster datasets in your R session.

You would need the stars::read_stars() function.

solution 

immigrants_cologne <-
  stars::read_stars("./data/immigrants_cologne.tif")

inhabitants_cologne <-
  stars::read_stars("./data/inhabitants_cologne.tif")

5

Aggregate the number of immigrants in each district using sum as function. Does it work?

It’s the same function as on slide XXX.

solution 

immigrants_count <-
  aggregate(
    x = immigrants_cologne,
    by = corona_cologne,
    FUN = sum,
    na.rm = TRUE
  )
## Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, : st_crs(x) == st_crs(y) ist nicht TRUE
# It does not work as R is complaining about non-matching CRS.

Let’s see what is happening here by comparing the exact CRS strings:

sf::st_crs(corona_cologne)
## 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["Statistical analysis."],
##         AREA["Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Turkey; United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State."],
##         BBOX[24.6,-35.58,84.17,44.83]],
##     ID["EPSG",3035]]
sf::st_crs(immigrants_cologne)
## Coordinate Reference System:
##   User input: ETRS89-extended / LAEA Europe (with axis order normalized for visualization) 
##   wkt:
## PROJCRS["ETRS89-extended / LAEA Europe (with axis order normalized for visualization)",
##     BASEGEOGCRS["ETRS89",
##         DATUM["European Terrestrial Reference System 1989",
##             ELLIPSOID["GRS 1980",6378137,298.257222101004,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     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]]]

Indeed, they are a bit different. We should change that.

6

Transform the raster data (again) to EPSG:3035.

stars and sf objects talk fluently to each other. You could either use sf::st_transform or stars::st_transform_proj() for this task. Try it out!

solution 

immigrants_cologne <-
  immigrants_cologne %>% 
  stars::st_transform_proj(3035)

inhabitants_cologne <-
  inhabitants_cologne %>% 
  stars::st_transform_proj(3035)

7

Aggregate the number of immigrants and inhabitants in each district using sum as function. Does it work now?

It’s the same function as on slide XXX.

solution 

immigrants_count <-
  aggregate(
    x = immigrants_cologne,
    by = corona_cologne,
    FUN = sum,
    na.rm = TRUE
  )

inhabitants_count <-
  aggregate(
    x = inhabitants_cologne,
    by = corona_cologne,
    FUN = sum,
    na.rm = TRUE
  )

8

Create an immigrant share variable from the two aggregated layers and add them to the corona_cologne dataset. Plot it if you will.

You can either do it on-the-fly using dplyr::mutate() as on the slides (no. 33), or create first a new object and then add it using base-R (corona_cologne$immigrant_share <- ...). You could even combine the approaches.

Either way, what might be a bit cumbersome is that you have to extract the first list element in each of the raster layers (e.g., immigrants_count[[1]])

solution 

corona_cologne <-
  corona_cologne %>% 
  dplyr::mutate(
    immigrant_share = immigrants_count[[1]] * 100 / inhabitants_count[[1]]
  )

plot(corona_cologne["immigrant_share"])