¿Intersección de dos polígonos con el mismo FID?

Question

Necesito realizar una intersección entre dos capas, una son bloques y la otra son buffers para cada bloque que son del mismo tamaño de la carretera. Los bloques y los buffers comparten el mismo FID.

Lo que intento obtener es la compacidad de cada bloque. He intentado utilizar un iterador, pero o bien tarda una eternidad o después de horas de procesamiento sólo me da la intersección de un par. También probé con tabulate intersect pero la salida. Una intersección regular no funciona porque algunos buffers van en bloques vecinos y sólo quiero la intersección entre polígonos que comparten el mismo FID.

Mis conocimientos de python son limitados pero estoy aprendiendo.

Answer 1

Esta solución utiliza python (w/ arcpy ) para emparejar la geometría de los bloques y los búferes en función de su FID realiza una intersección de esos pares y escribe la geometría resultante en un archivo shape junto con las características de origen. FID . La última función del script suma el porcentaje de cada búfer que está formado por su geometría de "intersección" correspondiente (es decir, la compacidad del bloque) y lo escribe en el campo compactnes del shapefile de la intersección.

Para adaptar esta secuencia de comandos a su proyecto, sólo tiene que establecer el parámetro blocks y buffers a las rutas de los archivos de sus conjuntos de datos y establezca las variables intersection la ruta en la que desea que se escriba la salida.

import os
import arcpy
from os import path
from arcpy import da, env, management

env.overwriteOutput = True

blocks = 'P:/ath/to/block/shapefile/block_shp_name.shp'
buffers = 'P:/ath/to/buffer/shapefile/buffer_shp_name.shp'
intersection = 'P:/ath/to/intersection/shapefile/intersection.shp'
compact_field = 'compactnes'

def createIntersectionFc():
    """Create feature class for the intersection features to be written to"""

    # get the spatial reference from the blocks feature class and use
    # it define that attribute for the intersection fc
    desc = arcpy.Describe(blocks)
    srs = desc.spatialReference

    shp_path = path.dirname(intersection)
    shp_name = path.basename(intersection)
    geom_type = 'POLYGON'
    management.CreateFeatureclass(shp_path, shp_name, 
        geom_type, spatial_reference=srs)

    # drop default feature class field and add field to hold the percentage
    # of the buffers that the intersection makes up
    f_type = 'DOUBLE'
    management.AddField(intersection, compact_field, f_type)

    drop_field = 'Id'
    management.DeleteField(intersection, drop_field)

def intersectFeatures():
    """Intersect the block and buffers geometries and write the result to a 
    feature class"""

    geom_mapping = {} 
    fields = ['OID@', 'SHAPE@']
    with da.SearchCursor(blocks, fields) as cursor:
        for oid, geom in cursor:
            geom_mapping[oid] = {'block': geom}

    with da.SearchCursor(buffers, fields) as cursor:
        for oid, geom in cursor:
            if oid in geom_mapping:
                geom_mapping[oid]['buffer'] = geom
            else:
                print 'There is no block match'
                print 'for the buffer with FID: {0}'.format(oid)

    i_cursor = da.InsertCursor(intersection, fields)

    for oid, geom_dict in geom_mapping.iteritems():
        intersect_geom = geom_dict['block'].intersect(geom_dict['buffer'], 4)
        i_cursor.insertRow((oid, intersect_geom))

    del i_cursor

def getBufferPercentage():
    """Get the percentage of each buffer that is made up by its corresponding
    intersection geometry and write it to the intersection fc"""

    area_dict = {}
    s_fields = ['OID@', 'SHAPE@AREA'] 
    with da.SearchCursor(buffers, s_fields) as cursor:
        for oid, area in cursor:
            area_dict[oid] = area

    u_fields = s_fields + [compact_field]
    with da.UpdateCursor(intersection, u_fields) as u_cursor:
        for oid, area, compact in u_cursor:
            compact = area / area_dict[oid]
            u_cursor.updateRow((oid, area, compact))

createIntersectionFc()
intersectFeatures()
getBufferPercentage()

Answer 2

No estoy seguro de si utiliza ArcGIS, pero si es así, puede hacerlo:

1. Une los búferes con los bloques y elige conservar sólo los registros coincidentes,
2. En el paso 1 sólo se obtienen los búferes que tienen fid coincidentes con los bloques y, a continuación, se pueden intersecar los nuevos búferes con los bloques.

Dependiendo del error que obtengas cuando intentes intersecar las capas, puede ser que tus conjuntos de datos sean demasiado grandes para ejecutar esto de una sola vez.