3

I need to perform an intersection between two layers, one is blocks and the other is buffers for each block which are of the same size of the road. Blocks and buffers share the same FID.

What I'm trying to obtain is the compactness of each block. I tried using an iterator but it either takes forever or after hours of processing gives me only the intersection for one pair. Also tried tabulate intersect but the output. A regular intersection doesn't work because some buffers go on neighbouring blocks and I only want the intersection between polygons sharing the same FID.

My python skills are limited but I'm learning.

Just a sample; my dataset is much larger than this

  • just a remark, but why don't you compute your compactness with area/perimeter^2 ? – radouxju Aug 19 '15 at 13:53
  • can you work with sql? – elmo Aug 19 '15 at 14:17
  • elmo, I can... not an expert but I kind of understand the syntax and have been using it; radoyxju, thanks for the suggestion, I'll have it in mind if I can't figure this out. – 1buzz Aug 19 '15 at 14:29
  • can you give more information please, blocks are corcles or rectangles – elmo Aug 19 '15 at 14:37
  • Block are the rectangles; although in some areas of the city I have circular blocks, in which case the compactness I'm trying to compute will be close or equal to 1, as the shape of the block will be the same as the shape of the circle buffer of the same area. – 1buzz Aug 19 '15 at 15:14
3

This solution uses python (w/ arcpy) to pair up the geometry of the blocks and buffers based on their FID, performs an intersection of those pairs, and writes the resultant geometry to a shapefile along with the source features' FID. The last function within the script adds the percentage of each buffer that is made up of its corresponding 'intersection' geometry (i.e. the block's compactness) and writes it to the field compactnes of the intersection shapefile.

To adapt this script to your project just set the blocks and buffers variables to the file paths of your datasets and set the intersection variable the path where you want the output to be written.

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()
1

I am not sure if you are using ArcGIS but if you are, you can:

  1. Join the buffers with the blocks, and choose keep matching records only,
  2. step 1 gives you only the buffers that have matching fid's with the blocks then you can intersect the new buffers with the blocks.

Depending on the error you get when you try to intersect the layers, it might be your datasets are too big to run this in one go.

  • I'm using ArcGIS 10.2.2 and also have QGIS Pisa. What I'm trying to get is the percentage of the block that falls inside a the buffer (circle) of the same area. There's a buffer for each block in my dataset. I'm not getting errors, but intersect just doesn't work because when a buffer is too big (because the block is big) it touches adjacent blocks. What I need is some kind of loop which takes the first FID from blocks, then looks for the same in Buffers, and intersect them; then moves to the second pair, and so on. I have 11000 polygons. – 1buzz Aug 19 '15 at 16:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.