1

Pick a address in Germany and geocode it through OSM. Make sure that the CRS is transformed to EPSG:3035.

If you’re not familiar with German places, you could, for example, choose between these three: - “B2,1 68159 Mannheim” - “Waldspirale, 64289 Darmstadt” - “Gustav-Stresemann-Ring 11, 65189 Wiesbaden”

solution 

gesis_cologne <-
  tmaptools::geocode_OSM(
    "Unter Sachsenhausen 6, 50667 Köln", 
    as.sf = TRUE,
    geometry = "point"
    ) %>% 
  sf::st_transform(3035) %>% 
  dplyr::select(query, geometry = point)

gesis_cologne
## Simple feature collection with 1 feature and 1 field
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 4106871 ymin: 3096821 xmax: 4106871 ymax: 3096821
## projected CRS:  ETRS89-extended / LAEA Europe
##                               query                geometry
## 1 Unter Sachsenhausen 6, 50667 Köln POINT (4106871 3096821)
bboxing <- function(input_point) {
  sf::st_bbox(
    c(
      xmin = sf::st_coordinates(input_point)[1] - 1000,
      xmax = sf::st_coordinates(input_point)[1] + 1000,
      ymax = sf::st_coordinates(input_point)[2] + 1000,
      ymin = sf::st_coordinates(input_point)[2] - 1000
    ),
    crs = sf::st_crs(3035)
  )
}

2

Apply the bboxing function and store the results in a new object. What do you think is the function doing?

Note that we use a cut-out map of a specific place’s surroundings.

solution 

gesis_bbox <- bboxing(gesis_cologne)

gesis_bbox
##    xmin    ymin    xmax    ymax 
## 4105871 3095821 4107871 3097821
# It's a bounding box for the areal surrounding of our geocoded address -- 1000 
# meters in each direction. We will use it to cut out a slice of some data we
# will load next.

3

Please download data on roads using OpenStreetMap. You can use the “highways” plus “trunk”, “primary”, “secondary”, and “tertiary” as key/value pairs. Again, make sure that the CRS is transformed to EPSG:3035

This task was part of the session this morning. It’s part of this crazy pipe where you first define a bounding box using osmdata::getbb() and so on.

solution 

cologne_roads <-
  osmdata::getbb("Köln") %>% 
  osmdata::opq(timeout = 25*100) %>%
  osmdata::add_osm_feature(
    "highway", 
    c("trunk", "primary", "secondary", "tertiary")
  )%>% 
  osmdata::osmdata_sf() %>% 
  .$osm_lines %>% 
  sf::st_transform(3035)

4

Crop the downloaded road data to the extent of the object you created in exercise “2”. We haven’t used the function before, but it’s the same as in the use of cropping raster data – it’s the sf::st_crop() function.

The first argument of the function is your downloaded road data, the second one the – psst! – bounding box from exercise “2”.

solution 

gesis_roads <-
  cologne_roads %>% 
  sf::st_crop(gesis_bbox)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

Use the following function for shearing and rotating your layers. Run the code first before applying it.

rotate_data <- function(data, x_add = 0, y_add = 0) {
  
  shear_matrix <- function () { 
    matrix(c(2, 1.2, 0, 1), 2, 2) 
  }
  
  rotate_matrix <- function(x) { 
    matrix(c(cos(x), sin(x), -sin(x), cos(x)), 2, 2) 
  }
  
  data %>% 
    dplyr::mutate(
      geometry = 
        .$geometry * shear_matrix() * rotate_matrix(pi/20) + c(x_add, y_add)
    )
}

5

Finally, plot the data as is, and then after or while rotating it with the rotate_data() function.

You can rotate the data either before calling gpplot() and store them in a separate object. Or, you can do that on the fly within the geom_sf() call.

solution 

# unrotated
ggplot() +
  geom_sf(data = gesis_roads)

# rotated
ggplot() +
  geom_sf(data = rotate_data(gesis_roads))