Sunday 23 September 2018

arcpy - Python script to iterate through geodatabase and assign field value based on spatial location


I have two types of spatial data. One is a geodatabase of parcel data by county in Florida, and the other is a shapefile containing wind_zones through out Florida (used for enforcing certain building codes). I want to create and populate two fields in the parcel data based on which county and which wind_zone they are in. I have found a post here: Assigning field value based on location without join that answers a similar question and is what I am basing my code off of, but I need to iterate this process through an entire geodatabase that contains ~70 counties. Is there a line of code that can be added to iterate through the geodatabase so I do not have to manually run the code for each county?


Here is the code I have so far:


import arcpy


# Set overwrite option
arcpy.env.overwriteOutput = True

# KEEPING ORIGINAL LAYERS (NOT CREATING ADDITIONAL "JOINED" LAYER)
# Create FeatureLayers
arcpy.MakeFeatureLayer_management("C:/Data/HAAS/parc_join/parcels_1.gdb", "lyr_parcels")
arcpy.MakeFeatureLayer_management("C:/wind_zones/All_WZs.shp", "lyr_windZones")

# Add a "LOCATED_WZ" & "WZ_CAT" field
arcpy.AddField_management("lyr_parcels", "Located_WZ", "TEXT", "", "", "20")

arcpy.AddField_management("lyr_parcels", "WZ_Cat", "SHORT")


# Create a search cursor for the WZs
rows = arcpy.SearchCursor("lyr_windZones")
for row in rows:
# What you'll do is select each WZ one at a time, and then select all the parcels in that WZ and calculate the "LOCATED_WZ" field
# NOTE: If you are using not using shapefiles, then you'll have to change the FID in the line below to OBJECTID (or similar)
arcpy.SelectLayerByAttribute_management("lyr_windZones", "NEW_SELECTION", "\"FID\" = " + str(row.getValue("FID")))
arcpy.SelectLayerByLocation_management("lyr_Parcels", "INTERSECT", "lyr_windZones", "", "NEW_SELECTION")

arcpy.CalculateField_management("lyr_parcels", "Located_WZ", "'{0}'".format(str(row.getValue("COUNTY"))), "PYTHON_9.3", "")
print "Finished processing " + str(row.getValue("COUNTY"))

# Create a search cursor for the WZs
rows = arcpy.SearchCursor("lyr_windZones")
for row in rows:
# What you'll do is select each WZ one at a time, and then select all the parcels in that WZ and calculate the fields
# NOTE: If you are using not using shapefiles, then you'll have to change the FID in the line below to OBJECTID (or similar)
arcpy.SelectLayerByAttribute_management("lyr_windZones", "NEW_SELECTION", "\"FID\" = " + str(row.getValue("FID")))
arcpy.SelectLayerByLocation_management("lyr_Parcels", "INTERSECT", "lyr_windZones", "", "NEW_SELECTION")

arcpy.CalculateField_management("lyr_parcels", "WZ_Cat", "'{0}'".format(str(row.getValue("DESCRIPT"))), "PYTHON_9.3", "")
print "Finished processing " + str(row.getValue("DESCRIPT"))

Because each county is it's own feature class, I suppose I don't need to add a new field to tell me county name, because I will already know which county it is in.



Answer



Regardless of the decision you take to choose how you'll add the information from the wind zones to your parcels in case a parcel intrsects more than one wind zone: you can use arcpy.ListFeatureClasses() to list all your parcels feature classes, then iterate over the feature classes with a for loop and run your code for each of them:


import arcpy

# Set overwrite option
arcpy.env.overwriteOutput = True


# Define your workspace. This should be your database with parcels feature class
arcpy.env.workspace = r"C\Data\HAAS\parc_join\parcels.gdb"

# This can be done once for all iterations
arcpy.MakeFeatureLayer_management("C:/wind_zones/All_WZs.shp", "lyr_windZones")

for fc in arcpy.ListFeatureClasses():

# KEEPING ORIGINAL LAYERS (NOT CREATING ADDITIONAL "JOINED" LAYER)

# Create FeatureLayers
arcpy.MakeFeatureLayer_management(fc, "lyr_parcels")
...

Be careful, there is an error in your code in this line:


arcpy.MakeFeatureLayer_management("C:/Data/HAAS/parc_join/parcels_1.gdb",    "lyr_parcels")

The input of MakeFeatureLayer is a feature class, not a gdb.


Other attention points:




  • use arcpy.da.Cursors with a with statement, not the old arcpy.Cursors, especially if you're overwriting them (using with deletes the cursor automatically when they've run).

  • you can calculate both fields with one cursor, no need to use 2.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...