Wednesday 29 June 2016

arcgis desktop - How to use iterators in ModelBuilder to create new feature?


I am attempting to use ModelBuilder to automate the process of plotting parcels of land that are related to mineral leases.


The existing lease feature class contains information related to reference document information, such as volume and page.


I am using a database populated with information pertaining to new leases, containing lease ID and reference deed volume and page.


Many of the new leases I plot refer to polygons that have already been plotted and refer to the same deed volume and page.


What I would like to do is create a model that will select features from the existing lease feature class based on the VOLUME and PAGE attributes in my database, copy/paste that feature into a separate feature class (editing layer from which features are loaded into the lease dataset), and populate the fields of the newly copied feature with the lease ID from my database.


I believe iterate row selection could be the best first step, but I am unsure where to go from there.


This is the most complicated ModelBuilder project I have worked on so any information or suggestions would be a huge help.



Answer




The solution is python using a cursor and dictionary to do the matching look-up and then output the result with an insert cursor to an empty feature class that has the same feature type and data structure as the original mineral lease with the additional Lease ID field from the table. The code to do the look-up is incredibly fast this way.


Here is some code that just needs to be adapted to your data (workspace, input feature class, related table, output feature class, and field names of the VOLUME, PAGE, and LEASEID). I assumed your data has a many-to-many relationship, where many parcels can have the same VOLUME and PAGE and many LeaseIDs can have the same VOLUME and PAGE, that way it will work even if you are dealing with a simpler relationship between the two data sources. The output will include all Mineral Lease features, including those that have no matching LeaseID from your table, and all matching features will be duplicated to create as many features as needed to hold every associated LeaseID.


I assumed you wanted to output all of the matched lease features to a single feature class, but if you actually wanted separate feature classes created for each set of Volume and Page values, then that could be done with some minor code modifications, which I could provide. I also could easily modify the code to eliminate the set of Mineral Lease features that have no matching LeaseIDs from being created if you don't actually want them.


I tested this code on my own data using a feature class with over 100,000 features and a many-to-many related table with over 130,000 features and it produced nearly 2,074,000 features that matched on the two relate fields I used. It took just 12 minutes and 17 seconds to create that many features (a rate of about 10,000 features every 3 seconds). Using a ModelBuilder iterator and repeated select queries would have taken hours or days to do the same thing. My code does not even require the related fields involved to be indexed (but you dare not do your approach in ModelBuilder without first indexing the related fields). If the unique values of my two fields had instead resulted in a One-to-One or Many-to-One relationship between the two data sources, then the code would have finished in about 2 minutes.


Open Python Idle from your ArcGIS Program Folder in the Python 2.7 group, open a New File, and paste this code into the new file then save the document with the name you want with the .py extension. Then modify the inputs and field names to match your data and run the code.


import datetime

def hms_string(sec_elapsed):
h = int(sec_elapsed / (60 * 60))
m = int((sec_elapsed % (60 * 60)) / 60)

s = sec_elapsed % 60.
return "{}:{:>02}:{:>05.2f}".format(h, m, s)
# End hms_string

def timedif(end_datetime, start_datetime):
seconds_elapsed = (end_datetime - start_datetime).total_seconds()
return hms_string(seconds_elapsed)
# End timedif

start_time = datetime.datetime.now()

print "Script start: {}".format(start_time)

import arcpy
from arcpy import env

print "Imports loaded. Time elapsed: {}".format(timedif(datetime.datetime.now(), start_time))


# ---Variables that you need to customize to match your data---



# Customize the workspace path to fit your data, for example a geodatabase
workspace = r"C:\Users\CRID\Documents\ArcGIS\YourGeodatabase.gdb"
env.workspace = workspace

# Customize this to match your feature class name (assumes this is in workspace)
inputFC = "MineralLeasesFC"
inputFCFull = workspace + '\\' + inputFC
print inputFCFull


# Customize this to match your lease ID table (assumes this is in workspace)
relatedTable = "LeaseIDsTable"
relatedTableFull = workspace + '\\' + relatedTable
print relatedTableFull

# Customize this to be the output FC name you want (assumes output to workspace)
outputFC = "MineralLeasesWithIDs"
outputFCFull = workspace + '\\' + outputFC
print outputFCFull


# Customize this to be the name of the field names for the VOLUME, PAGE and LEASEID
inVolumeFld = "VOLUME"
inPageFld = "PAGE"
relatedVolumeFld = "VOLUME"
relatedPageFld = "PAGE"
leaseIDFld = "LEASEID"
print "Variables loaded. Time elapsed: {}".format(timedif(datetime.datetime.now(), start_time))


# ---Code You Should Not Need to Modify---



# Check to see if the output already exists and delete it if it does
if arcpy.Exists(outputFCFull):
arcpy.Delete_management(outputFCFull, "FeatureClass")
print "Deleted {}. Time elapsed: {}".format(outputFC, timedif(datetime.datetime.now(), start_time))

# Create the output Feature Class based on the input
arcpy.CreateFeatureclass_management(workspace, outputFC, "POLYGON", inputFCFull, "SAME_AS_TEMPLATE", "SAME_AS_TEMPLATE", inputFCFull)
print "Created {}. Time elapsed: {}".format(outputFC, timedif(datetime.datetime.now(), start_time))


# get the field information for the LeaseID field and add it to the output
fields = arcpy.ListFields(relatedTable)
for field in fields:
if field.name.upper() == leaseIDFld.upper():
fldType = field.type
fldPrecision = field.precision
fldScale = field.scale
fldLength = field.length
fldAliasName = field.aliasName

# Add a field for the LeaseID
arcpy.AddField_management(outputFCFull, leaseIDFld, fldType, fldPrecision, fldScale, fldLength, fldAliasName)
print "Added {} Field. Time elapsed: {}".format(leaseIDFld, timedif(datetime.datetime.now(), start_time))

# Set up the relatedTable field list
relatedFieldsList = [relatedVolumeFld, relatedPageFld, leaseIDFld]

# Build a dictionary from a da SearchCursor from the related table
relatedDict = {}
with arcpy.da.SearchCursor(relatedTableFull, relatedFieldsList) as searchRows:

for searchRow in searchRows:
keyValue = '{};{}'.format(searchRow[0], searchRow[1])
if not keyValue in relatedDict:
relatedDict[keyValue] = [searchRow[2]]
else:
relatedDict[keyValue].append(searchRow[2])
del searchRows, searchRow
print "Loaded relatedDict from {}. Time elapsed: {}".format(relatedTable, timedif(datetime.datetime.now(), start_time))

# Build a field list for the editable fields of the input and output

inputFieldsList = [inVolumeFld, inPageFld]
outputFieldsList = []
fields = arcpy.ListFields(inputFCFull)
for field in fields:
if field.editable and field.type != 'Geometry':
inputFieldsList.append(field.name)
outputFieldsList.append(field.name)
elif field.type == 'Geometry':
inputFieldsList.append("SHAPE@")
outputFieldsList.append("SHAPE@")

outputFieldsList.append(leaseIDFld)
print "Field list created for {}. Time elapsed: {}".format(inputFC, timedif(datetime.datetime.now(), start_time))

# Build a dictionary from a da SearchCursor
inputDict = {}
with arcpy.da.SearchCursor(inputFCFull, inputFieldsList) as searchRows:
for searchRow in searchRows:
keyValue = '{};{}'.format(searchRow[0], searchRow[1])
if not keyValue in inputDict:
inputDict[keyValue] = [list(searchRow[2:])]

else:
inputDict[keyValue].append(list(searchRow[2:]))
del searchRows, searchRow
print "Loaded inputDict from {}. Time elapsed: {}".format(inputFC, timedif(datetime.datetime.now(), start_time))

counter = 0
cursor = arcpy.da.InsertCursor(outputFCFull, outputFieldsList)
for key in inputDict:
rows = inputDict[key]
if key in relatedDict:

for row in rows:
for lease in relatedDict[key]:
counter += 1
if counter % 10000 == 0:
print "Has a lease. Counter = {}. Time elapsed: {}".format(counter, timedif(datetime.datetime.now(), start_time))
relatedRow = list(row)
relatedRow.append(lease)
cursor.insertRow(relatedRow)
else:
for row in rows:

counter += 1
if counter % 10000 == 0:
print "Has no lease. Counter = {}. Time elapsed: {}".format(counter, timedif(datetime.datetime.now(), start_time))
row.append(None)
cursor.insertRow(row)
del cursor, inputDict, relatedDict
print "Inserted {} Records. Time elapsed: {}".format(counter, timedif(datetime.datetime.now(), start_time))
print "Script Finished: {}".format(datetime.datetime.now())

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...