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@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo