[R-sig-Geo] How can I manipulate the parameters of st_interpolate_aw so it doesn't fail as a result of point and line intersections?

[Edzer Pebesma]

Thanks for the clear diagnosis; a small reprex would be

library(sf)
a <- st_polygon(list(cbind(c(0,0,7.5,7.5,0),c(0,-1,-1,0,0))))
b <- st_polygon(list(cbind(c(0,1,2,3,4,5,6,7,7,0),c(1,0,.5,0,0,0.5,
        -0.5,-0.5,1,1))))
c<-st_sfc(a,b)

aa = st_sf(value = 1, geom = c[1])
bb = st_sf(value = 2, geom = c[2])
cc = st_sf(value = 3, geom = c[1] + c(.5,.5))
st_interpolate_aw(aa, cc, extensive = TRUE)

cc = st_sf(value = 3, geom = c[1] + c(1,1))
st_interpolate_aw(aa, cc, extensive = TRUE)

st_interpolate_aw(aa, bb, extensive = TRUE)


revealing two problems: 0 or 1 dimensionsal intersections, and
intersections giving GEOMETRYCOLLECTION.

This now should work when you install sf from github.

On 9/4/19 8:32 AM, Chris Fowler wrote:
>  I have also asked this question here:
> https://stackoverflow.com/questions/57767022/how-do-you-use-st-interpolate-aw-with-polygon-layers-that-legitimately-include-p?noredirect=1#comment101987833_57767022
> 
> but expect that the specialized knowledge on this list serve may achieve a
> better result and that members of this list who are, like me, still making
> the transition from sp to sf could benefit from the responses it generates.
> 
> I am trying to use st_interpolate_aw() to assign values from one polygon
> layer to another polygon layer (voting district vote totals assigned to
> census tracts). The operation fails because st_interpolate_aw() calls
> st_intersection() which finds point and line intersections that then get
> treated as points, lines, or geometrycollections and break the subsequent
> code casting geometries to multipolygons (I get the error "Error in
> st_cast.POINT(x[1], to, ...) :cannot create MULTIPOLYGON from POINT")
> 
> Within st_interpolate_aw the relevant code looks like this where 'x' is my
> first polygon layer and 'to' is my second:
> 
> i = st_intersection(st_geometry(x), st_geometry(to))
> idx = attr(i, "idx")
> i = st_cast(i, "MULTIPOLYGON")
> 
> When I use debug to step through the above lines of code in
> st_interpolate_aw I can look at summary(i) after running st_intersection
> and I can see multiple POINTS, GEOMETRYCOLLECTIONS and LINESTRINGS that
> were created in addition to the desired POLYGONS and MULTIPOLYGONS. This is
> not a problem of messy layers per se--we should expect point and line
> intersections when comparing these kinds of administrative boundaries, but
> they shouldn't be relevant for interpolation purposes and could be
> reasonably dropped. The problem is that st_cast doesn't seem set up to make
> allowances.
> 
> I expect the answer to this question is to figure out what parameters I can
> set via ... in st_interpolate_aw to alter the behavior of st_intersection()
> or st_cast() to make this work. Alternatively, there may be a path through
> tricks like st_set_precision() or st_snap() to pre-process my polygons so
> this doesn't happen (see somewhat relevant back and forth here
> <https://github.com/r-spatial/sf/issues/860>). I have already run
> st_is_valid() on my polygons and they are fine. Also, I fundamentally think
> that st_intersection is behaving appropriately with regards to the data,
> and that it is st_cast not being willing to drop geometries with no area
> that causes the problem.
> 
> The example below doesn't get me all the way to the solution as it deals
> with the code internal to st_interpolate_aw not the function itself, but it
> is small, draws on polygons from the vignette
> <https://cran.r-project.org/web/packages/sf/vignettes/sf1.html>, and
> generates the same error as my code when running st_interpolate_aw(). If
> the example below could be run with a parameter change to st_intersection
> or st_cast so that st_cast doesn't return an error I
> think I could get st_interpolate_aw to work
> 
> a <- st_polygon(list(cbind(c(0,0,7.5,7.5,0),c(0,-1,-1,0,0))))
> b <- st_polygon(list(cbind(c(0,1,2,3,4,5,6,7,7,0),c(1,0,.5,0,0,0.5,
>         -0.5,-0.5,1,1))))
> c<-st_sfc(a,b)
> plot(c)
> i <- st_intersection(c[[1]],c[[2]])
> i=st_cast(i,"MULTIPOLYGON")
> 
> Thanks,
> Chris
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

-------------- next part --------------
A non-text attachment was scrubbed...
Name: pEpkey.asc
Type: application/pgp-keys
Size: 3110 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20190904/5159acf9/attachment.bin>