Identifying internal polygons using Arcpy

This is just a small python script that will identify any polygons that are within other polygons in a shapefile. It has come handy for me for a project I was working on where I had to identify the “internal” polygons. The main trick here is to compare each polygon with the outer boundary of every other

# -*- coding: utf-8 -*-
#!/usr/bin/env python
import arcpy
#-------------------------------------------------------------------------------
# Name:        flag_int_polygons
# Purpose:     Identifies interior polygons
#
#-------------------------------------------------------------------------------
def main():

    output_folder = u"C:\temp"
    inFC=u"C:\temp\polys.shp"

    # Check if the compare shapefile already exists and if so, recreate it
    if arcpy.Exists(output_folder+"/tmpCompare.shp")==True:
        arcpy.Delete_management(output_folder+"/tmpCompare.shp")
    arcpy.CopyFeatures_management(inFC, output_folder+"/tmpCompare.shp")

    try:
        arcpy.AddField_management(output_folder+"/tmpCompare.shp", "LABEL", "TEXT", "", "", 2, "", "", "REQUIRED")
    except:
        pass

     # Get shape column
    shapeName = arcpy.Describe(inFC).shapeFieldName
    shapeName1 = arcpy.Describe(output_folder+"/tmpCompare.shp").shapeFieldName

     # Loop through input features
    rows=arcpy.SearchCursor(inFC)

    # Start loop through input fc
    r=1
    for row in rows:
        partnum = 0
        parcel=row.getValue(shapeName)
        exBoundary=getPolyBoundary(parcel)
        fid=row.getValue("FID")
        rows1=arcpy.UpdateCursor(output_folder+"/tmpCompare.shp")
        for row1 in rows1:
            fid1=row1.getValue("FID")
            parcel1=row1.getValue(shapeName1)
            if fid1 <> fid:
                if parcel1.within(exBoundary)==True:
                    print "Feature " + str(fid1) + " is within feature " + str(fid)
            rows1.updateRow(row1)

        del row1
        del rows1
    del row
    del rows
#-------------------------------------------------------------------------------
# Name:        getPolyBoundary
# Purpose:     Returns the external polygon boundary
#
#-------------------------------------------------------------------------------
def getPolyBoundary(geom):
    partno=0
    pntObj = arcpy.Point()
    arrayObj = arcpy.Array()

    for part in geom:
        pntCounter=0;
        # Loop through points in part
        for pnt in geom.getPart(partno):
            if pnt:
                pntObj.X=pnt.X
                pntObj.Y=pnt.Y
                arrayObj.add(pntObj)
            else:
                break
        pntCounter=pntCounter+1
    bound=arcpy.Polygon(arrayObj)
    #print str(bound)
    return bound
if __name__ == '__main__':
    main()
Advertisements

Convert a Polyline shapefile to PolylineM using arcpy

I am sure I am re-inventing the wheel here but I couldn’t find any ready-made scripts to do this. Plus, I thought, it was high time I get my feet wet with Python. So- behold of my first python script(!). Its quite simple- it basically takes an input linear shapefile and converts it to a new shapefile with M values. If a field containing the measured line length is given, it will use it to calibrate the M values – otherwise it will use the map (geometry) length so the M value of the last vertex of each line will equal the length of the line.

This is the first of a bunch of ‘useful’ Avenue scripts that I used quite often in the past and converted to Python. Will eventually convert the others and learn more about Python and arcpy in the process.

The script is listed below for your viewing pleasure. It is heavily commented so it should be self explanatory. Or, you can simply download it from here

import arcpy, math, os

# Input feature class
#
inFeatures="C:/mapdata/osm_greek_data/projected/road_sample.shp"

# Fieldname containing the measured length for each line
# if null, then the map length will be used
mlField="meas_len"

# Output feature class settings
out_path = "C:/mapdata/osm_greek_data/projected"
out_name = "road_sampleM.shp"
geometry_type = "POLYLINE"
template = inFeatures
spatial_reference=inFeatures
has_m = "ENABLED"
has_z = "DISABLED"

# Create the new shapefile
arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type, template, has_m, has_z, spatial_reference)

#
# Create a list of fields to use when copying field values from input to output
#
fields = arcpy.ListFields(inFeatures)
fieldList = []
for field in fields:
  fieldList.append(field.name)

# Loop through input features
rows=arcpy.SearchCursor(inFeatures)

# Output features
outFeatures=out_path+"/"+out_name

# Open insertcursor
#
outRows = arcpy.InsertCursor(outFeatures)

# Create point and array objects
#
pntObj = arcpy.Point()
arrayObj = arcpy.Array()

shapeName=arcpy.Describe(inFeatures).shapeFieldName #Get shape field name
featType=arcpy.Describe(inFeatures).shapeType       # Get feature type
isMeasured=arcpy.Describe(inFeatures).hasM  # Check if its a polylineM feature already

# Record counter
i=1

# Loop through input shapefile
for row in rows:
    print "Processing row %d " % i
    feat=row.getValue(shapeName)
    noparts=feat.partCount
    partno = 0
    currentM=0
    # Loop through geometry parts
    for part in feat:
        pntCounter=0;
        segmentDistance=0
        prevX=0
        prevY=0
        # Loop through points in part
        for pnt in feat.getPart(partno):
          if pntCounter == 0: # First point
            prevX=pnt.X
            prevY=pnt.Y
            pntObj.X=pnt.X
            pntObj.Y=pnt.Y
            pntObj.M=0
          else:
            curX=pnt.X
            curY=pnt.Y
            # Calculate segment length
            segmentDistance = math.sqrt(math.pow((curX-prevX),2) + math.pow((curY-prevY),2))
            # Use the measured distance
            if mlField <> "":
                segmentDistance=segmentDistance/feat.length*row.getValue(mlField)
            # Add previous M to segment length
            currentM=currentM + segmentDistance
            pntObj.X=pnt.X
            pntObj.Y=pnt.Y
            pntObj.M=currentM
            prevX=pnt.X
            prevY=pnt.Y
          # Add measured points to array
          arrayObj.add(pntObj)
          pntCounter=pntCounter+1
    # Create new row
    #
    outFeat=outRows.newRow()

    # Copy attributes to new row
    #
    for fieldName in fieldList:
      if fieldName not in  ("FID", "Shape"): # ignore FID and Shape fields
        outFeat.setValue(fieldName, row.getValue(fieldName))

    # Assign array of PointMs to new feature
    outFeat.shape=arrayObj

    # Insert feature
    outRows.insertRow(outFeat)

    # Clear array of points
    #
    arrayObj.removeAll()
    i=i+1

# Remove cursor from memory
del outRows
del rows