Más

Recortando un SpatialPoygonsDataFrame a otro SpatialPolygonsDataFrame

Recortando un SpatialPoygonsDataFrame a otro SpatialPolygonsDataFrame


¿Cómo restrinjo unSpatialPolygonsDataFramepara que ninguno de los polígonos quede fuera de un "subyacente"SpatialPolygonsDataFrame?

La aplicación que tengo en mente involucra Voronoi Polygons a la Carson Farmer: el resultado de dicha descomposición es (normalmente) rectangular por defecto, pero a menudo es más atractivo visualmente adaptar el resultado a una geometría subyacente (ciudad / condado / frontera del país, etc.). Otra aplicación común es eliminar Alaska y Hawái de las geometrías de todo EE. UU. Cuando solo los 48 inferiores son de interés.

Esto es similar a esta pregunta, pero más complicado, ya que ambas geometrías de entrada constan de varios polígonos: un simplegIntersecciónresultará en una disolución innecesaria de la geometría principal de interés.

También es similar a esta pregunta, pero no pude seguir exactamente lo que estaba sucediendo allí porque no había un ejemplo reproducible. Creo que el siguiente ejercicio probablemente subsume esta pregunta.

Aquí hay un ejemplo simple con el que trabajar:

biblioteca (sp) biblioteca (rgeos) biblioteca (maptools) puntos <-list ("A" = c (0, 0), "B" = c (1/2, 0), "C" = c (1, 0 ), "D" = c (0, 1/2), "E" = c (1/2, 1/2), "F" = c (1, 1/2), "G" = c (0 , 1), "H" = c (1/2, 1), "I" = c (1, 1), "J" = c (0, -1/2), "K" = c (3 / 2, -1/2), "L" = c (7/24, 3/8), "M" = c (5/8, 3/8), "N" = c (3/8, 5 / 8), "O" = c (5/8, 17/24), "P" = c (3/2, 1)) pts_poly <- función (x) Reducir (rbind, lapply (x, function (y) unlist (puntos [y]))) poly1_list <- list ("poly11" = list (ID = "A", coords = pts_poly (c ("A", "B", "E", "D"))) , "poly12" = lista (ID = "B", coords = pts_poly (c ("B", "C", "F", "E"))), "poly13" = lista (ID = "C", coords = pts_poly (c ("E", "F", "I", "H"))), "poly14" = lista (ID = "D", coords = pts_poly (c ("D", "E" , "H", "G")))) poly2_list <- list ("poly11" = list (ID = "a", coords = pts_poly (c ("J", "M", "L"))), "poly12" = lista (ID = "b", coords = pts_poly (c ("J", "K", "M"))), "poly13" = lista (ID = "c", coords = pts_poly (c ("M", "K", "P"))), "poly14" = lista (ID = "d", coords = pts_poly (c ("M", "P", "O"))), " poly15 "= lista (ID = "e", coords = pts_poly (c ("L", "M", "O", "N")))) poly1 <-SpatialPolygons (lapply (poly1_list, function (pl) {with (pl, Polygons (list (Polygon (coords)), ID = ID))})) df1 <- data.frame (var1 = rnorm (length (poly1)), var2 = rnorm (length (poly1)), row.names = sapply (poly1 @ polygons, function (x) x @ ID)) polydf1 <- SpatialPolygonsDataFrame (poly1, df1) poly2 <-SpatialPolygons (lapply (poly2_list, function (pl) {with (pl, Polygons (list (Polygon (coords)) , ID = ID))})) df2 <- data.frame (var1 = rnorm (length (poly2)), var2 = rnorm (length (poly2)), row.names = sapply (poly2 @ polygons, function (x) x @ ID)) polydf2 <- SpatialPolygonsDataFrame (poly2, df2)

Mucho código. Esto es lo que tenemos visualmente:

(código para parcelas :)

par (mfrow = c (1, 3)) par (mar = c (0,0,0,0)) par (oma = c (0,0,0,0)) plot (polydf1, xlim = c (0 , 1.1), ylim = c (0, 1.1)) título ("Polígonos subyacentes", línea = -15, cex.main = 1.8) texto (Reducir (rbind, puntos [1: 9]) + .03, etiquetas = nombres (puntos) [1: 9], col = "rojo") texto (matriz (c (.25, .25, .75, .25, .75, .75, .25, .75), ncol = 2 , byrow = TRUE), cex = 2, labels = sapply (polydf1 @ polygons, slot, name = "ID")) plot (polydf2, xlim = c (0, 1.6), ylim = c (-. 5, 1.1) ) título ("Polígonos que se recortarán", línea = -15, cex.main = 1.8) texto (Reducir (rbind, puntos [- (1: 9)]) + .03, etiquetas = nombres (puntos) [- ( 1: 9)], col = "red") texto (coordenadas (polydf2), cex = 2, etiquetas = sapply (polydf2 @ polygons, slot, name = "ID")) plot (polydf2) plot (polydf1, add = TRUE) title ("Superposición de polígonos", línea = -15, cex.main = 1.8)

Lo que quiero es eliminar las partes de los triángulos E, F, G y H que se encuentran fuera del primer conjunto de polígonos, A, B, C y D.

Una simple intersección es incorrecta:

plot (gIntersection (polydf1, polydf2, byid = TRUE), main = "Intersección simple: Overkill")

Así es como se ve mi salida deseada (básicamente eliminando los bordes de A, B, C y D):


El enfoque más simple para hacer esto es

biblioteca (raster) x <- recorte (polydf2, polydf1)

Esto hará la intersección geométrica, pero también asegurará que los atributos estén en orden. Ver? 'paquete de trama'(sección XIV) para obtener más funciones que se ocupan de la superposición de polígonos.


La respuesta es bastante simple una vez que la ve, pero pasé medio día golpeándome la cabeza contra el escritorio antes de darme cuenta del simple ajuste que debemos hacer; No pude encontrarlo en ningún otro lugar en línea, así que quería dejar esto aquí para referencia de cualquier otra persona que luche como yo.

La clave es eliminar los bordes del primer conjunto de polígonos antes de la intersección con los polígonos importantes:

clipped_polys <- gIntersection (polydf2, gUnaryUnion (polydf1), byid = TRUE, id = sapply (polydf2 @ polygons, slot, name = "ID"))

Para asegurarnos de que los polígonos conservan sus nombres originales, también especificamos la ranura de identificación nosotros mismos. Aquí está el resultado:

Y, por supuesto, podemos volver a agregar cualquier dato asociado conpoli2fácilmente como se señala en esta respuesta:

clipped_polys <- SpatialPolygonsDataFrame (clipped_polys, polydf2 @ data)

Ver el vídeo: R studio: extract raster values by XY coordinate and line