Disolver/Unificador Mal Comportado/Polígonos Irregulares en R

Question

Estoy trabajando en un proyecto en el que participa, que reúne varios temas dando vueltas difícil encontrar partes de una imagen. La salida de los usuarios es bastante impresionante.

Genial, derecho?

Sin embargo, lo que quiero es obtener el área de la imagen que tiene 1, 2, 3, 4, etc. los usuarios de la selección. Y estoy usando R para el análisis con el sp, raster, y maptools paquetes (con un poco de rgeos tiro en cuando es necesario). Esencialmente, mi flujo de trabajo es 1) crear un único polígono de cada entrada del usuario (de hecho tengo la pista de cada selección individual por usuario, pero parecía más fácil de combinar en un único SpatialPolygon por usuario) 2) combinación de todos ellos juntos en un solo Objeto SpatialPolygonsDataFrame 3) Rasterizar el objeto mediante la función contar 4) Evaluar el área de la trama con n número de usuarios de seleccionarlo

Aquí está la salida de SpatialPolygonsDataFrame a jugar con

Esto funciona muy bien, pero es slooow debido a la rasterización. Siento que debería ser capaz de hacer algo con el SpatialPolygons objeto. Y he probado un par de cosas con la unión, gIntersects, gUnion, y se cruzan. Sin embargo, me siguen saliendo errores como los siguientes

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, "rgeos_difference") : 
TopologyException: Input geom 1 is invalid: Self-intersection at or near point

Justo lo suficiente. A menudo un determinado usuarios SpatialPolygons no son válidos cuando me probarlos mediante gIsValid. Mirando, por ejemplo, sólo los dos primeros, esto se vuelve obvio que hay muchos superposición de puntos.

Usted puede ver los puntos de intersección de una manera bastante clara. Además, en poly2 (el azul) no parece ser un poco de auto-intersección, causando gIsValid para devolver false.

También he intentado unionSpatialPolygons (con avoidRGEOS=T), lo que crea un objeto único unionSpatialPolygons(polysData[1:2,], IDs=names(polysData[1:2,]@polygons)) Pero las fronteras todavía se superponen como en la imagen de arriba, en lugar de ser una suave polígono exterior de la frontera - es decir, los polígonos no son disueltos juntos. A continuación, el uso de nuevas funciones en este nuevo SpatialPolygons objeto, tengo el mismo problemas.

Así que, ¿hay una manera de utilizar los polígonos en lugar de ir a los rásteres? Seguro que haría mi vida más rápido, lo cual sería excelente. Me siento como este debe ser un problema común con una rueda inventado que todavía tengo que encontrar. Los pensamientos?

Answer 1

Sin los datos originales, no puedo estar seguro de que esto funcionará, pero pensé que podría ayudarte. Yo no lo trajo todo el camino hasta allí, esta solución todavía, es probable que necesita un cierto nivel de automatización, pero podría dar un camino a seguir

En primer lugar, creo que algunos espacial de los polígonos

polypoints1 <- matrix(c(1,2,2,1,1,2,2,1,1,2),ncol=2)
polypoints2 <- matrix(c(1,3,3,1,1,3,3,1,1,3),ncol=2)
polypoints3 <- matrix(c(1,2,2,1,1,2,2,1,1,2)+1.1,ncol=2)
polypoints4 <- matrix(c(1,2,2,1,1,2,2,1,1,2)+0.5,ncol=2)

p1 <- Polygon(polypoints1)
ps1 <- Polygons(list(p1),1)
sps1 <- SpatialPolygons(list(ps1))

p2 <- Polygon(polypoints2)
ps2 <- Polygons(list(p2),2)
sps2 <- SpatialPolygons(list(ps2))

p3 <- Polygon(polypoints3)
ps3 <- Polygons(list(p3),3)
sps3 <- SpatialPolygons(list(ps3))

p4 <- Polygon(polypoints4)
ps4 <- Polygons(list(p4),4)
sps4 <- SpatialPolygons(list(ps4))

Yo les trazan sólo para ver

plot(sps2,col='green')
plot(sps1,add=T,col='blue')
plot(sps3,add=T,col='yellow')
plot(sps4,add=T,col='purple')

Me fundí en un spdf

data=data.frame(c(x=rep(1,4)),row.names=c(1:4))
sps <- SpatialPolygons(list(ps1,ps2,ps3,ps4))
spdf <- SpatialPolygonsDataFrame(sps,data)

Puede identificar qué polígono se superpone que así:

gIntersects(spdf,spdf,byid =T)

Desde el comando de arriba podría crear algún tipo de bucles para hacer la superposición de las combinaciones de abajo (solo estoy ignorando sps4 por razones de brevedad en este punto)

poly2a <- gIntersection(spdf[2,],spdf[1,],drop_lower_td=T)
poly2a <- SpatialPolygonsDataFrame(poly2a,data.frame(c(x=1),row.names=c(1)))
plot(poly2a,add=T,col='red')

Esta vez tenemos que cambiar el ID ya que vamos a rbind estos más tarde

poly2b <- gIntersection(spdf[2,],spdf[3,],drop_lower_td=T)
poly2b <- spChFIDs(poly2b,"2")
poly2b <- SpatialPolygonsDataFrame(poly2b,data.frame(c(x=1),row.names=c(2)))
plot(poly2b,add=T,col='red')

Combinar la superposición de polígonos en otro spdf

spdf_overlaps <- rbind(poly2a,poly2b)
poly2 <- unionSpatialPolygons(spdf_overlaps,rep(1,2))
plot(poly2,add=T,col='blue')

Ahora tenemos poly2 que es donde tenemos 2 capas superpuestas (excepto combinaciones con sps4), a continuación, la figura 3 capas, solo tenemos que ver dónde poly2 y spdf se superponen (en el caso de una forma más automatizada versión de esto, usted necesita para asegurarse de que 'poly2' no incluye sps4 como en este ejemplo)

gIntersects(poly2,spdf[4,],byid =T)

poly3 <- gIntersection(poly2,spdf[4,],drop_lower_td=T)
plot(poly3,add=T,col="red") 

Check it out

gIsValid(poly2)
gIsValid(poly3)

Alternativamente, usted puede siempre hacer un pseudo rasterización, mucho más fácil, pero usted pierde un poco de detalle dependiendo del tamaño de la celda:

Primero hacer la red:

bb <- bbox(spdf)
cs <- c(0.1,0.1)  # cell size
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)


sp_grd <- SpatialGridDataFrame(grd,
                           data=data.frame(id=1:prod(cd)))

A continuación, realice la cuadrícula en un polígono que se utiliza para la superposición

library(Grid2Polygons)
grid <- Grid2Polygons(sp_grd)
plot(grid)

Contar el número de polígonos que se superponen en cada celda de la cuadrícula

count <- apply(gContains(spdf,grid,byid=T),1,sum)

Finalmente, la trama!

plot(grid)
for(i in 1:length(grid)){
    plot(grid[i,],col=rev(heat.colors(3))[count[i]],add=T)
}

Answer 2

Suena como la limpieza de la geometría hasta le permitirá estar en el camino a su no-trama solución ... aquí algo de un parche que no ayuda con la fijación de la geometría mal:

# Load the library and problematic data
library(rgeos)
load("oneImage2_spdf.Rdata")

>gIsValid(polysData, reason = T)
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, "rgeos_union") : 
TopologyException: Input geom 0 is invalid: Self-intersection at or near point -119.84228271000001 34.349193950783807 at -119.84228271000001 34.349193950783807

Que el infractor de datos; tal vez de un solo punto, pero haciendo un búfer de 0 anchura debería solucionar el problema.

# Let's buffer by 0:
polysData <- gBuffer(polysData, width=0, byid = T)

gBuffer no le gusta esto - hay un no proyectados CRS unido a él!

Warning message:
In gBuffer(polysData, width = 0, byid = T) :
  Spatial object is not projected; GEOS expects planar coordinates

Tiempo para engañar a rgeos:

# Set to a projected CRS - could be anything
polysCRS <- polysData@proj4string
polysData@proj4string <- CRS("+init=epsg:32737")

# Now buffer it but preserve individual objects
polysData <- gBuffer(polysData, width = 0, byid = T)

# Reset the original CRS
polysData@proj4string <- polysCRS

Compruebe si la geometría es válido:

>gIsValid(polysData, reason = T)
"Valid Geometry"