st_interpolate_aw returns error - replacement has x rows, data has y

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

Solution

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