Tengo un conjunto de polígonos que representan grandes áreas, por ejemplo, barrios de la ciudad. Quiero identificar las grandes áreas superpuestas entre ellos.
Pero hay un problema: a veces estos polígonos se superponen a lo largo de sus perímetros (porque fueron dibujados con poca precisión). Esto generará superposiciones largas y estrechas que no me importan.
Pero otras veces habrá grandes superposiciones de polígonos robustos, lo que significa grandes áreas donde el polígono de un vecindario se superpone a otro. Quiero seleccionar solo estos.
Vea la imagen a continuación de solo las superposiciones. Imagine que quisiera seleccionar solo el polígono azul en la esquina inferior izquierda.
Podría mirar áreas, pero a veces las angostas son tan largas que terminan teniendo áreas tan grandes como el polígono azul. Intenté hacer una relación de área / perímetro, pero eso también arrojó resultados mixtos.
Incluso he intentado usar ST_MinimumClearance
, pero a veces las áreas grandes tendrán una parte estrecha unida o dos vértices muy cercanos.
¿Alguna idea de otros enfoques?
Al final, lo que mejor funcionó para mí fue usar un búfer negativo, como lo sugieren @Cyril y @FGreg a continuación.
Usé algo como:
ST_Area(ST_Buffer(geom, -10)) as neg_buffer_area
En mi caso, las unidades eran metros, entonces 10 m de amortiguación negativa.
Para polígonos estrechos, esta área devolvió cero (también, la geometría estaría vacía). Luego usé esta columna para filtrar los polígonos estrechos.
Respuestas:
Intentaría crear un búfer negativo, si come polígonos delgados, entonces es bueno, si no come el polígono, entonces es mío ... :-)
ejecute este script, habiendo establecido previamente 2/3 del ancho de los polígonos lineales ...
OS: -) ...
fuente
ST_Area(ST_Buffer(geom, -10))
, -10 siendo -10 metros en mi caso. Si algo devuelve 0 de esa expresión, entonces podría filtrarlo.En lugar de área / perímetro, es mejor usar el área dividida por el cuadrado del perímetro (o su inverso).
Esto también se llama "índice de forma". El cuadrado del perímetro dividido por el área tiene un valor mínimo de 4 * Pi () (en el caso de un disco, que es la geometría 2D más compacta), por lo que puede normalizarse por 4 * Pi () para un fácil interpretación (los valores normalizados cercanos a 1 significan que tiene objetos muy compactos y los cuadrados tienen valores de aproximadamente 1.27).
EDITAR: Un umbral en el área sería útil para eliminar los artefactos muy pequeños, que podrían ser compactos. Entonces el índice de forma mostraría un mejor contraste. EDITAR: además de esta respuesta, el uso de ST_Snap podría ayudarlo a resolver el problema antes de que ocurra.
fuente
(o.overlap_perimeter^2 / o.overlap_area) / (4 * Pi()) as overlap_ratio
? Esto está teniendo peores resultados para mí que solo el área / perímetro.o.overlap_perimeter / (4 * sqrt(o.overlap_area)) as overlap_ratio
acuerdo con este documento, pero aún peores resultados (aunque es difícil cuantificar lo que quiero decir con peor) isprs-ann-photogramm-remote-sens-spatial-inf-sci.net/I-7/135/… , página 183.Una opción sería usar la razón del área del polígono a la línea más larga que se puede dibujar usando sus extremidades. Identificando polígonos largos y estrechos.
select * from polygons where ST_Length(ST_LongestLine(geom, geom)) < ST_Area(geom) * 4
Esto funciona bastante bien para los polígonos astillados. Puede ajustar la proporción (con qué multiplica el área) para satisfacer sus necesidades y proyección.
fuente
Parece que esto podría coincidir con su caso de uso: elimine los polígonos seleccionados
Parece que le gustaría probar la opción "Límite común más grande".
fuente
Esto me parece un caso de uso perfecto para la extensión de topología PostGIS . El parámetro de tolerancia de la topología determinará hasta qué punto permite que los vértices se ajusten a otros polígonos existentes, para hacer frente a la baja precisión de los datos de origen y limpiarlos.
En resumen, la estrategia es:
1. Habilite la extensión de topología
2. Crear una nueva topología vacía
El tercer parámetro es la tolerancia, en las unidades del CRS; elígelo sabiamente. Idealmente, desea un CRS donde la unidad es metros. Si la unidad CRS no tiene medidores, como con WGS 84, también conocido como 4326, úselo
ST_Transform
para reproyectar sus polígonos.3. Agregue una columna TopoGeometry a la tabla de polígonos
Esto devuelve un nuevo
layer_id
. Guárdelo, será necesario más tarde. Será una capa1
si comienza desde cero y se incrementará en cada nueva llamada.4. Agregue todos los polígonos a la topología.
Esto puede llevar varias horas para un gran conjunto de datos, tenga paciencia.
1
es el layer_id devuelto anteriormente.5. Encuentra caras que aparecen en varios vecindarios
Encuentre todas las caras de la topología que están presentes en 2 o más topogeometrías. Dejaré la consulta como ejercicio. Lo más fácil es probablemente con la
GetTopoGeomElements
función, luego agrupar por identificación facial y mirar las que cuentan 2 o más. Alternativamente, puede crear una nueva tabla con la geometría limpia de la columna topogeom, simplemente convertirla en geometría estándartopogeom::geometry
y repetir lo que ya tiene ahora, pero ahora con un conjunto de datos limpio sin superposiciones de la astilla.fuente