I am working to re-aggregate population counts from one administrative level to another. One particular row keeps returning the error "replacement has 4 rows, data has 2".
I have tried validating the geometries, but it has not worked. Neither did a buffer with 0 distance.
I was preparing the dummy data for posting this question using dput(). After saving the dput() of the sf layers into a text file, the command line works after re-assigning the variables. This only works when using the dput output saved in a textfile, but it does not work if I directly do
st_interpolate_aw(dput(geography1),dput(geography2), extensive = TRUE)
.
This returns the same error.
Is there any known reason why the error might be happening, and how saving the dput() output into a .txt file might be removing the error?
note - Is there any way to share the features other than dput()? the high amount of vertices in one of the layers surpasses the character limit for the question
edit I think I managed to export a geojson with the error, in case you are interested in checking it: https://drive.google.com/file/d/1wF6AB1oHVEkcuI4iRAloLO56DL_8EZBr/view?usp=sharing
CodePudding user response:
I believe the issue lies with the nature of the intersection of the two objects; it seems that the intersection produces two polygon type objects (perfectly OK) and two line type objects; these have by definition zero area which kind of messes up the interpolation process.
You have two options:
- tweak a bit the code behind
sf::st_interpolate_aw()
to force a filter to polygon-geometries behind the scenes - switch the processing engine from shiny new
{s2}
to trusty old GEOS
I suggest doing the latter, as shown in the code bellow, but making a pull request to {sf}
would be the cowboy thing to do :)
library(sf)
library(dplyr)
big <- st_read("SA1_suburb_1221.geojson")
small <- st_read("suburb_1221.geojson")
st_intersection(st_geometry(big), st_geometry(small))
# Geometry set for 4 features
# Geometry type: GEOMETRY
# Dimension: XY
# Bounding box: xmin: 142.8129 ymin: -36.75028 xmax: 142.9655 ymax: -36.63092
# Geodetic CRS: WGS 84
# LINESTRING (142.9648 -36.75015, 142.9648 -36.75...
# POLYGON ((142.8673 -36.63148, 142.8371 -36.6315...
# MULTILINESTRING ((142.9648 -36.74569, 142.9648 ...
# MULTIPOLYGON (((142.9515 -36.74478, 142.9515 -3...
sf_use_s2(F)
st_intersection(st_geometry(big), st_geometry(small))
# although coordinates are longitude/latitude, st_intersection assumes that they are planar
# Geometry set for 2 features
# Geometry type: GEOMETRY
# Dimension: XY
# Bounding box: xmin: 142.8129 ymin: -36.75028 xmax: 142.9655 ymax: -36.63092
# Geodetic CRS: WGS 84
# POLYGON ((142.9648 -36.75015, 142.9645 -36.7499...
# GEOMETRYCOLLECTION (POLYGON ((142.9504 -36.7467...
result <- st_interpolate_aw(big["below_poverty"],
st_geometry(small),
extensive = F)
result
# // a perfectly reasoneable result