library - ¿Cómo puedo hacer una unión espacial con el paquete sf usando st_join()

sf r

Question

Aquí hay un ejemplo de juguete con el que he estado luchando.

# Make points
point1 <- c(.5, .5)
point2 <- c(.6, .6)
point3 <- c(3, 3)
mpt <- st_multipoint(rbind(point1, point2, point3)) # create multipoint
# Make polygons
square1 <- rbind(c(0, 0), c(1, 0), c(1,1), c(0, 1), c(0, 0))
square2 <- rbind(c(0, 0), c(2, 0), c(2,2), c(0, 2), c(0, 0))
square3 <- rbind(c(0, 0), c(-1, 0), c(-1,-1), c(0, -1), c(0, 0))
mpol <- st_multipolygon(list(list(square1), list(square2), list(square2))) # create multipolygon
# Convert to class'' sf''
pts <- st_sf(st_sfc(mpt))
polys <- st_sf(st_sfc(mpol))
# Determine which points fall inside which polygons
st_join(pts, polys, join = st_contains)

La última línea produce

Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) :
cannot coerce class "c("sfc_MULTIPOINT", "sfc")" to a data.frame

¿Cómo puedo hacer una unión espacial para determinar qué puntos caen dentro de qué polígonos?

Answer 1

También estoy trabajando en las características del paquete sf , así que pido disculpas si esto no es correcto o si hay mejores maneras. Creo que un problema aquí es que si construyes geometrías como en tu ejemplo, no estás obteniendo lo que piensas:

> pts
Simple feature collection with 1 feature and 0 fields
geometry type: MULTIPOINT
dimension: XY
bbox: xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID): NA
proj4string: NA
st_sfc.mpt.
1 MULTIPOINT(0.5 0.5, 0.6 0.6...
> polys
Simple feature collection with 1 feature and 0 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 0 ymin: 0 xmax: 2 ymax: 2
epsg (SRID): NA
proj4string: NA
st_sfc.mpol.
1 MULTIPOLYGON(((0 0, 1 0, 1 ...

Puedes ver que solo tienes una "característica" tanto en pts como en polys . Esto significa que está creando una característica de "multipolígono" (es decir, un polígono constituido por 3 partes), en lugar de tres polígonos diferentes. Lo mismo ocurre con los puntos.

Después de cavar un poco, encontré esta forma diferente (y en mi opinión más fácil) de construir las geometrías, usando la notación WKT:

polys <- st_as_sfc(c("POLYGON((0 0 , 0 1 , 1 1 , 1 0, 0 0))",
"POLYGON((0 0 , 0 2 , 2 2 , 2 0, 0 0 ))",
"POLYGON((0 0 , 0 -1 , -1 -1 , -1 0, 0 0))")) %>%
st_sf(ID = paste0("poly", 1:3))
pts <- st_as_sfc(c("POINT(0.5 0.5)",
"POINT(0.6 0.6)",
"POINT(3 3)")) %>%
st_sf(ID = paste0("point", 1:3))
> polys
Simple feature collection with 3 features and 1 field
geometry type: POLYGON
dimension: XY
bbox: xmin: -1 ymin: -1 xmax: 2 ymax: 2
epsg (SRID): NA
proj4string: NA
ID .
1 poly1 POLYGON((0 0, 0 1, 1 1, 1 0...
2 poly2 POLYGON((0 0, 0 2, 2 2, 2 0...
3 poly3 POLYGON((0 0, 0 -1, -1 -1, ...
> pts
Simple feature collection with 3 features and 1 field
geometry type: POINT
dimension: XY
bbox: xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID): NA
proj4string: NA
ID .
1 point1 POINT(0.5 0.5)
2 point2 POINT(0.6 0.6)
3 point3 POINT(3 3)

Puedes ver que ahora tanto polys como pts tienen tres características.

Ahora podemos encontrar la "matriz de intersección" usando:

# Determine which points fall inside which polygons
pi <- st_contains(polys,pts, sparse = F) %>%
as.data.frame() %>%
mutate(polys = polys$ID) %>%
select(dim(pi)[2],1:dim(pi)[1])
colnames(pi)[2:dim(pi)[2]] = levels(pts$ID)
> pi
polys point1 point2 point3
1 poly1 TRUE TRUE FALSE
2 poly2 TRUE TRUE FALSE
3 poly3 FALSE FALSE FALSE

es decir (como se señaló @symbolixau en los comentarios) que los polígonos 1 y 2 contienen los puntos 1 y 2, mientras que el polígono 3 no contiene ningún punto. El punto 3 no está contenido en ningún polígono.

HTH.

Answer 2

Tenga en cuenta que el conjunto original de multipunto y multipolígono se puede ''convertir'' a punto y polígono, sin crear nuevos objetos:

st_contains(polys %>% st_cast("POLYGON"), pts %>% st_cast("POINT"), sparse = F)
# [,1] [,2] [,3]
#[1,] TRUE TRUE FALSE
#[2,] TRUE TRUE FALSE
#[3,] FALSE FALSE FALSE

Answer 3

Veo una salida diferente:

> # Determine which points fall inside which polygons
> st_join(pts, polys, join = st_contains)
Simple feature collection with 1 feature and 0 fields
geometry type: MULTIPOINT
dimension: XY
bbox: xmin: 0.5 ymin: 0.5 xmax: 3 ymax: 3
epsg (SRID): NA
proj4string: NA
geometry
1 MULTIPOINT(0.5 0.5, 0.6 0.6...

Fue esto con la versión más reciente de CRAN de sf ?