Porcentaje de polígono en un archivo de forma dentro del polígono de otro

13

Soy un novato, disculpas si esto es obvio / ya me han preguntado y respondido pero no pude encontrar nada.

Tengo dos shapefiles: 1. una capa límite administrativa para un condado en el Reino Unido conocida como límite LSOA que tiene 500 pequeñas zonas en ella 2. una zona de inundación.

Idealmente, quiero averiguar cuáles de las pequeñas zonas de LSOA están ≥50% dentro de la zona de inundación y terminan con un sí / no o un 1/0 para cada una de las 500 zonas de LSOA.

Pero no sé cómo hacer esto. Pensé que podría unir los dos archivos de forma, pero no hay un atributo común entre ellos. Luego pensé que podría usar la función Unir atributo por ubicación, que funcionó y me muestra qué LSOA están en la zona de inundación, pero eso es casi todo (ver imagen 2).

Creo que este es un problema de SQL pero no lo sé. Soy nuevo en QGIS y nunca he usado PostgreSQL.

Cualquier ayuda sería muy apreciada. Puedo proporcionar cualquier información que necesiten ustedes para ayudarme.

ingrese la descripción de la imagen aquí

ingrese la descripción de la imagen aquí

KJGarbutt
fuente

Respuestas:

11

Esta es una tarea relativamente simple utilizando las herramientas de geoprocesamiento incluidas en QGIS.

  1. Calcule el área de sus zonas LSOA.

    • Abra la tabla de atributos de la capa LSOA.
    • Habilitar el modo de edición.
    • Abre la calculadora de campo.
    • Cree un nuevo campo de tipo "Número decimal (real)" con la expresión "$ area".
    • Desactiva el modo de edición (guardar ediciones).
  2. Combine la capa de la zona de inundación en una única entidad de varias partes.

    • Vector > Geometry Tools > Singleparts to Multipart.
    • Seleccione "--- Fusionar todo ---" para el campo ID único.
  3. Intersecte la capa de zona LSOA con la capa de zona de inundación multiparte.

    • Vector > Geoprocessing Tools > Intersect.
    • La capa de entrada son las zonas LSOA, la capa de intersección son las zonas de inundación.
  4. La capa resultante serán las partes de las zonas LSOA (con los atributos de la capa de zonas LSOA) que se superponen con la capa de zonas de inundación. Para calcular la proporción de cada zona LSOA dentro de una zona de inundación:

    • Calcule el área de las entidades intersectadas (como en el paso 1), luego
    • Agregue otro campo, dividiendo el área original (total) por el área intersectada. El resultado es un decimal entre 0 y 1. Multiplica por 100 para dar un porcentaje.
  5. Unir la capa LSOA original a la capa intersectada, utilizando la ID única compartida por ambas capas.

  6. Exporte la capa unida como un nuevo archivo de forma.

  7. Eliminar los atributos duplicados.

Et voilà!

Sin el paso 2, se crearía una característica individual para cada característica de zona de inundación diferente para cada característica de LSOA. Probablemente esto no sea lo que desea si solo está interesado en la cobertura total para cada zona LSOA. Si desea diferenciar entre inundaciones fluviales / de marea / pluviales (y los datos de la zona de inundación lo admiten), puede convertir una sola parte en multiparte especificando el campo "TIPO" como el campo ID único.

Snorfalorpagus
fuente
¡Gracias por tu ayuda! Muy apreciado. Sin embargo, estoy teniendo algunos problemas. He seguido los pasos. El paso 3, el Intersect, tardó 10 horas en completarse y cuando lo hice, todo lo que obtuve fue un archivo de formas vacío: i.imgur.com/QIM6Gtg.png ¿Hay algo que me haya perdido? Intenté completar el proceso y seguir el paso 4, pero no hay datos para calcular el área de intersección.
KJGarbutt
He tenido problemas para hacer intersecciones con capas de inundación antes. Las características son grandes y complicadas. La forma en que he trabajado en el pasado es dividirlos en características más pequeñas, para que el índice espacial pueda hacer más trabajo. Para hacer esto, cree una cuadrícula de vector de la misma extensión que la capa de inundación ( Vector > Research Tools > Vector Grid... Output grid as polygons), luego intersecte la cuadrícula con la capa de inundación. Luego use la salida en lugar de la capa de inundación en el paso 3. Supongo que la razón por la que la capa estaba vacía es porque se bloqueó.
Snorfalorpagus
Gracias de nuevo. El único problema ahora es que QGIS se bloquea cada vez que intento crear la cuadrícula de vectores. He seguido el consejo desde aquí, pero se bloquea cada vez. Cambié los parámetros varias veces e intenté usar solo el archivo de forma de zona de inundación, en lugar de tener todo el archivo de mi proyecto abierto y falla cada vez. ¿Algunas ideas? ! Captura de pantalla aquí .
KJGarbutt
Los parámetros X e Y que ha especificado son demasiado pequeños. Pruebe algo como 1000 x 1000. Incluso puede hacer esto varias veces, es decir, primero hacer 5000 x 5000, usar la salida para crear 500 x 500. Vea la respuesta relacionada aquí: gis.stackexchange.com/a/66319/12420
Snorfalorpagus
¡Casi lo he roto con tu ayuda! Sin embargo, cuando voy a unir la capa LSOA original con la capa intersectada, pierdo muchos datos. Creo que se debe a que algunos de los cuadrados de cuadrícula de vectores creados caen dentro de la misma área LSOA y, por lo tanto, tienen el mismo código LSOA que cada uno. Por lo tanto, termino con más de 2 cifras porcentuales para cada área de LSOA cuando hago la unión y parece que solo obtengo una de ellas. ¿Hay alguna manera de sumar cada porcentaje para cada cuadrado de la cuadrícula de vectores con el mismo LSOA?
KJGarbutt
6

Puede usar spaceialite y algunas funciones espaciales de SQL.

Select t1.geometry, t1.ID, area(t1.geometry), area(t2.geometry) ...... (anything you need to have in the table results)

(area(intersection(t1.geometry,t2.geometry))) as "Commun_AREA"

, ("Commun_AREA"*100/(area(t1.geometry))) as "Percent_AREA"

From lsoa as t1, flood_zone as t2

Where Intersects( t1.geometry,t2.geometry ) = 1
Cirilo
fuente
3

Esto parece algo que podría hacerse mucho más fácilmente que las respuestas enviadas. Usaría un script simple de Python personalmente:

floodName = "the layer name here"
boundryName = "the layer name here"
fieldName = "the name of the field to contain the output 1/0"
minCoverage = 0.5 # the minimum amount of area covered to write 1
updateMap = [] # this will store values to be written    

# get layers
floodLayer = QgsMapLayerRegistry.instance().mapLayersByName(floodName)[0]
boundryLayer = QgsMapLayerRegistry.instance().mapLayersByName(boundryName)[0]
fieldIndex = boundryLayer.dataProvider().fieldNameIndex(fieldName)    

# iterate through boundries
for b in boundryLayer.getFeatures():
    # get only flood features that intersect with this feature's bounding box
    # this will make the script go way faster than it would otherwise
    request = QgsFeatureRequest().setFilterRect(b.geometry().boundingBox())
    floodGeom = geometry()
    floodFeat = QgsFeature()
    iter = floodLayer.getFeatures(request)
    iter.nextFeature(feat)
    while iter.nextFeature(feat):
        floodGeom = floodGeom.combine(feat.geometry())
    intersectGeom = b.geometry().intersection(feat.geometry())
    if intersectGeom.area() > minCoverage * b.geometry().area():
        updateMap[b.id()] = {fieldIndex : 1}
    else:
        updateMap[b.id()] = {fieldIndex : 0}

boundryLayer.dataProvider().changeAttributeValues(updateMap)

esto solo evalúa los polígonos de inundación que se cruzan con el cuadro delimitador de cada capa límite, por lo que debería ser bastante rápido de ejecutar, luego solo actualiza un campo en la capa existente (en lugar de una operación compleja de hacer una capa completamente nueva y copiar valores antiguos luego eliminando)

Jesse McMillan
fuente
2

Tuve el mismo problema que KJ siguiendo las instrucciones de Snorfalorpagus usando el método "Intersecar" en el Paso 3. Me tomó bastante tiempo calcularlo y lo que me quedó quedó en blanco.

Intenté seguir los mismos pasos, excepto usar el método "Clip" en QGIS en lugar de Intersect, por lo que, en su ejemplo, lo que quedaría serían las partes de las áreas que NO están cubiertas por la zona de inundación. Esto pareció funcionar por alguna razón y pude usar el cálculo de campo "Área" del paso anterior, más un nuevo cálculo de "Área" en las partes restantes de cada polígono, para calcular el% de cada área que NO estaba cubierto por la otra capa de polígono.

Eso es técnicamente lo contrario de lo que pediste. Pero a partir de ahí, solo es cuestión de restar cada valor de 1 para obtener lo que cubre la zona de inundación.

escarbador de vidas ajenas
fuente