Raster extent based on SpatialPolygonsDataFrame generating incorrect extent

Question

I'm trying to rasterize data from Berlin's average price per parcel data. I've downloaded the shapefile from their data portal using an API with a function that generates an sf object, and then selected only the ID, the value (column BRW), and the geometry. here is the description of the object:

brw2003
Simple feature collection with 802 features and 2 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 13.10962 ymin: 52.33956 xmax: 13.74145 ymax: 52.66066
geographic CRS: WGS 84
First 10 features:
 gml_id BRW geometry
1 s_brw_2003.1001 160 MULTIPOLYGON (((13.53115 52...
2 s_brw_2003.1002 200 MULTIPOLYGON (((13.2037 52....
3 s_brw_2003.1003 80 MULTIPOLYGON (((13.56818 52...
4 s_brw_2003.1004 1500 MULTIPOLYGON (((13.42251 52...
5 s_brw_2003.1005 460 MULTIPOLYGON (((13.34401 52...
6 s_brw_2003.1007 400 MULTIPOLYGON (((13.58039 52...
7 s_brw_2003.1008 1800 MULTIPOLYGON (((13.43833 52...
8 s_brw_2003.1009 550 MULTIPOLYGON (((13.39865 52...
9 s_brw_2003.1010 140 MULTIPOLYGON (((13.50077 52...
10 s_brw_2003.1012 900 MULTIPOLYGON (((13.40922 52...

Then, in order to rasterize this data, I've tried to make a template raster based on the average block size of Berlin at 4.5 hectares = 45000 m^2. However, the raster doesn't crop correctly based on this code:

# changing to projected coordinate system
brw2003_ext <- st_transform(brw2003, 27700)
template03 <- raster(extent(brw2003_ext), resolution = 4500, crs = st_crs(brw2003_ext)$proj4string)
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
 Discarded datum Unknown based on Airy 1830 ellipsoid in Proj4 definition
> template03
class : RasterLayer 
dimensions : 8, 10, 80 (nrow, ncol, ncell)
resolution : 4500, 4500 (x, y)
extent : 1422360, 1467360, 383443.4, 419443.4 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs 

I don't understand why the xmin/xmax/ymin/ymax are all incorrectly placed, and why, at a resolution of only 4500 m^2, only 80 cells are generated.

I've wondered whether this was an issue because I was generating the extent based on a multipolygon or an sf object rather than an sp object, but neither option actually changes the outcome:

brw2003_ext <- st_union(brw2003_ext)
> brw2003_ext
Geometry set for 1 feature 
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 1422360 ymin: 385091.6 xmax: 1466630 ymax: 419443.4
projected CRS: OSGB 1936 / British National Grid
MULTIPOLYGON (((1445891 396577.3, 1445935 39656...
brw2003_ext <- as_Spatial(brw2003, cast = TRUE)
class : SpatialPolygonsDataFrame 
features : 802 
extent : 13.10962, 13.74145, 52.33956, 52.66066 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs 
variables : 2
names : gml_id, BRW 
min values : s_brw_2003.1001, 70 
max values : s_brw_2003.2221, 18000 
template03 <- raster(extent(brw2003_ext), resolution = 4500, crs = st_crs(brw2003_ext)$proj4string)
> template03
class : RasterLayer 
dimensions : 1, 1, 1 (nrow, ncol, ncell)
resolution : 4500, 4500 (x, y)
extent : 13.10962, 4513.11, -4447.339, 52.66066 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs 

This still gives me an incorrectly cropped raster. Any ideas on how to fix this?

Edit: incase anyone wanted to play with the data themselves, here is my code for downloading this dataset:

########################################
# creating functions to download data from FIS Broker 
########################################
get_X_Y_coordinates <- function(x) {
 sftype <- as.character(sf::st_geometry_type(x, by_geometry = FALSE))
 if(sftype == "POINT") {
 xy <- as.data.frame(sf::st_coordinates(x))
 dplyr::bind_cols(x, xy)
 } else {
 x
 }
}
sf_fisbroker <- function(url) {
 typenames <- basename(url)
 url <- httr::parse_url(url)
 url$query <- list(service = "wfs",
 version = "2.0.0",
 request = "GetFeature",
 srsName = "EPSG:25833",
 TYPENAMES = typenames)
 request <- httr::build_url(url)
 print(request)
 out <- sf::read_sf(request)
 out <- sf::st_transform(out, 4326)
 out <- get_X_Y_coordinates(out)
 out <- st_as_sf(as.data.frame(out))
 return(out)
}
export_format <- c(
 "geojson", 
 "sqlite"
)
sf_save <- function(z, fname) {
 ifelse(!dir.exists(fname), dir.create(fname), "Folder exists already")
 ff <- paste(file.path(fname, fname), export_format, sep = ".")
 purrr::walk(ff, ~{ sf::st_write(z, .x, delete_dsn = TRUE)})
 saveRDS(z, paste0(file.path(fname, fname), ".rds"))
}
brw2003 <- sf_fisbroker("https://fbinter.stadt-berlin.de/fb/wfs/data/senstadt/s_brw_2003")

Answer 1

You could use the stars package for rasterization:

library(sf)
library(stars)
brw2003 <- sf_fisbroker("https://fbinter.stadt-berlin.de/fb/wfs/data/senstadt/s_brw_2003")
brw2003_ext <- st_transform(brw2003, 27700)
brw2003_grid <- st_rasterize(brw2003_ext["BRW"], dx = 4500, dy = 4500)

BTW: 4500 does not create cell of size 4,500 m2, but cells of size 4,500 x 4,500 m = 20250000 m2. You should have cells of size ~67m to create cell of size 4,500 m2.

• cschwab98: Thanks for pointing me towards this — I found it pretty late the night I was working on this project after I managed to get things working, but I definitely plan to use it in the future. I also realized that I was using the resolution incorrectly haha, ended up using sqrt(45000).

Answer 2

Not totally an answer necessarily as I still don't know why this isn't working, but I managed a workaround by using a single polygon representing Berlin and including filter(!st_is_empty(.)) on the shapefile containing parcel values:

brw2003 <- sf_fisbroker("https://fbinter.stadt-berlin.de/fb/wfs/data/senstadt/s_brw_2003")
brw2003 <- brw2003 %>%
 dplyr::filter(NUTZUNG != "G - Gewerbe") %>%
 dplyr::select(gml_id, BRW, geometry) %>%
 filter(!st_is_empty(.))
load(url("https://userpage.fu-berlin.de/soga/300/30100_data_sets/berlin_district.RData"))
berlin.sf <- st_transform(berlin.sf, 3035)
berlin_template <- raster(extent(berlin.sf), 
 resolution = 216, 
 crs = st_crs(berlin.sf)$proj4string)
brw2003rast <- rasterize(brw2003, berlin_template, field = "BRW")
> brw2003rast
class : RasterLayer 
dimensions : 171, 211, 36081 (nrow, ncol, ncell)
resolution : 216, 216 (x, y)
extent : 4531042, 4576618, 3253780, 3290716 (xmin, xmax, ymin, ymax)
crs : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
source : memory
names : layer 
values : 70, 18000 (min, max)