How can I manipulate the parameters of st_interpolate_aw so it doesn't fail as a result of point and line intersections?

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

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
 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: How can I manipulate the parameters of st_interpolate_aw so it doesn't fail as a result of point and line intersections?

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