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")

Jot eN · Accepted Answer · 2021-03-21 13:06:49Z

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.

answered Mar 21, 2021 at 13:06

Jot eN

8268 silver badges21 bronze badges

1

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).
– cschwab98
Mar 23, 2021 at 23:33

Add a comment |

cschwab98 · Accepted Answer · 2021-03-20 02:41:57Z

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)

Stack Exchange Network

current community

your communities

more stack exchange communities

Raster extent based on SpatialPolygonsDataFrame generating incorrect extent

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged raster r sf rasterization sp or ask your own question.

Hot Network Questions

Raster extent based on SpatialPolygonsDataFrame generating incorrect extent

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged raster r sf rasterization sp or ask your own question.

Related

Hot Network Questions