libreria library r spatial sf

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



sf r (3)

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?


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.


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


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 ?