Polígonos parcelarios en R. Solución GIS Stackexchange
Muchos de los procesos básicos disponibles en los programas típicos (ArcMap/QGIS) tienen sus homólogos en el entorno R. Con un mínimo esfuerzo, se pueden desarrollar herramientas propias para resultados muy satisfactorios. Si eres de los que dicen que desde que se inventó el modelbuilder escribir código es tontería te diré lo que le dijo El Gran Lewosky a Jesús Quintana: Well…that’s your opinion man.
Este segundo episodio surge a raiz de una pregunta sin contestar del maravilloso foro de los freaks geográficos GIS Stack Exchange. Básicamente se busca un proceso donde, poder crear áreas de asignación por “fachada”, de forma que los polígonos internos de una serie de edificaciones tengan sus áreas de influencia recortadas de forma recta. Aquí el esquema del proceso:
En primer lugar cargaremos las librerías necesarias para el proceso y leeremos una capa de polígonos. En este caso unos bonitos polígonos dibujados a mano sobre la fabulosa ciudad de Santiago de Compostela, donde el sol brilla por su ausencia. Sobre ellos se han establecido unos buffers iniciales de 10 metros, por probar.
library(sf)
library(mapview)
library(rgeos)
library(htmlwidgets)
# Leer layer
pol <- st_read('../../../../00_BASE/POLS/Parcels.gpkg')
# Buffer
polbuf <- st_buffer(pol, 10, 0)
Podrían preocuparnos quizá esos abruptos buffers esquineros, pero eso no es lo que nos atañe (tendré que buscar herramientas en R para poder hacer el estirado de polígonos más que un simple buffer). Lo importante es que ya podemos ver las zonas de solape de buffers sobre las que habrá que buscar la forma de dividirlas para asignarlas a una edificación.
Una de las estrategias que podría ser interesante sería dividir los polígonos de solape por los puntos de medios de intersección entre los mismos, para después utilizar esos puntos como un rayo láser de corte de los polígonos. ¿Cómo calculamos los puntos medios de intersección entre polígonos? Pues así…
PUNTO DE INSPIRACIÓN Y O U ROBO 01
p_intersect <- gIntersection(as(as_Spatial(polbuf)[1,], "SpatialLines"), as_Spatial(polbuf))
# extraer coordenadas
line_1 <- st_as_sf(as(as_Spatial(polbuf), "SpatialLines"))
poly_1 <- st_as_sf(as(as_Spatial(polbuf), "SpatialLines"))
pnts <- st_intersection(line_1, st_cast(poly_1, "MULTILINESTRING", group_or_split = FALSE))
# capturar puntos y mantener solo los unicos
ps <- unique(pnts[st_geometry_type(pnts)=='MULTIPOINT',])
Es el momento de conectar los diodos a los puntos y generar una línea entre ellos que sirva para dividir cualquier polígono que se pase por el medio. La línea es tan poderosa que puede destruir cualquier cosa, por eso aquí se muestra en azul que todo el mundo pueda disfrutar de ella y no tener problemas de visión el día de mañana.
PUNTO DE INSPIRACIÓN Y O U ROBO 02
# Crear líneas con puntos
ls <- st_cast(ps, "LINESTRING")
Llegamos al momento más bonito del proceso en el cual dividimos con el laser esos malditos buffers solapados de la discordia.
PUNTO DE INSPIRACIÓN Y O U ROBO 03
# cortar polígonos con líneas y consolidad
lpi <- gIntersection(as_Spatial(polbuf), as_Spatial(ls))
blpi <- gBuffer(lpi, width = 0.000001)
dpi <- gDifference(as_Spatial(polbuf), blpi)
Ya tenemos los polígonos recortados en base a las fachadas o parcelas iniciales. Ahora sólo queda asignarles correctamente la vivienda a la que corresponden. Tan sencillo como eso.
polbuf2 <- st_join(st_cast(st_as_sf(dpi), 'POLYGON'), pol)
Y alguno se preguntará una vez más si es que esto no puede ser realizado en QGIS/ArcGIS…lo desconozco. La pregunta llevaba 5 meses abierta, con 7 deliciosos votos y ninguna respuesta, así que ya alguien la hubiera respondido si fuese tan sencillo.
