[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?

Chris Fowler chr|@ @end|ng |rom c@|ow|er@com
Wed Sep 4 08:32:22 CEST 2019


 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]]



More information about the R-sig-Geo mailing list