Gracias por los comentarios. La solución que hemos encontrado se publica a continuación.
# This script pulls the number of acres specified for conversion from a shapefile of all counties,
# it then gets the suitability raster for the state and the field polygon layer for the county and
# performs zonal statistics, determining the mean suitability for each farm.
# It then moves through the polygons from most suitable to least suitable and selects them
# until a threshold area set in the county shapefile is reached. It then exports these polygons to a
# new shapefile. Minimum parcel size can be set using the SQL statement in the UpdateCursor command.
# If they have the same suitability score, it selects larger pracels first.
# Import tools, licenses, and declares environment settings
import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
# Variables
County_acreage = r"C"
County_parcels = r"C:"
State_suit = r"C:"
ZonalTable = r"C:"
State_FIPS = ""
# Get number of harvested acres and then loop through all counties. SQL statement can limit extent.
rows = arcpy.SearchCursor(County_acreage,"STATE_FIPS = '27'")
for row in rows:
# Reset variables and find out how many acres are predicted
totalArea = 0
State_FIPS = row.getValue('STATE')+"_"+row.getValue('FIPS')
County_parcels = r"C:"
State_suit = r"C:"
HarvestAcres = float(row.getValue('DOE_45'))
HarvestSqMeters = HarvestAcres * 4046.86 #Converts acerage to square meters
print HarvestSqMeters
# Defines file path for the county being processed and state suitability raster
State_suit = State_suit+str(row.getValue('STATE'))+".tif"
County_parcels = County_parcels+str(row.getValue('FIPS'))
print "Processing: "+State_FIPS
# Makes parcels file a FeatureLayer to enable selecting and other operations
arcpy.MakeFeatureLayer_management(County_parcels,"FeatureLayer")
FeatureLayer = "FeatureLayer"
# Creates table with mean suitability score by field, joins it back to fields featurelayer
ZonalStatisticsAsTable(FeatureLayer, "Id", State_suit, ZonalTable, "DATA", "MEAN")
arcpy.JoinField_management (FeatureLayer, "OBJECTID", ZonalTable, "ID")
# Iterate through the parcles
# The first argument in the UpdateCursor is a SQL where clause, can be used to limit quality and size.
# The last argument in UpdateCursor sorts the parcels by descending mean suitability score, then size.
lines = arcpy.UpdateCursor(FeatureLayer, "COUNT > 5", "","", "MEAN D, COUNT D")
for line in lines:
totalArea = totalArea + line.AREA
line.grid_code = 5 #Used to label parcels that meet the requirements, can be any nubmer
lines.updateRow(line)
if totalArea > HarvestSqMeters:
break
arcpy.SelectLayerByAttribute_management(FeatureLayer, "ADD_TO_SELECTION", 'grid_code = 5')
arcpy.CopyFeatures_management(FeatureLayer, r"C:"+State_FIPS)
arcpy.Delete_management(FeatureLayer)
print "Completed: "+State_FIPS