# How can I manipulate the parameters of st_interpolate_aw so it doesn't fail as a result of point and line intersections? Classic List Threaded 2 messages 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) bb = st_sf(value = 2, geom = c) cc = st_sf(value = 3, geom = c + c(.5,.5)) st_interpolate_aw(aa, cc, extensive = TRUE) cc = st_sf(value = 3, geom = c + 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, 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 > ). 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 > , 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[],c[]) > i=st_cast(i,"MULTIPOLYGON") > > Thanks, > Chris > > [[alternative HTML version deleted]] > > _______________________________________________ > R-sig-Geo mailing list > [hidden email] > https://stat.ethz.ch/mailman/listinfo/r-sig-geo> -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo pEpkey.asc (3K) Download Attachment