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