Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I need to make a map that randomly places parcels (groups of pixels or polygons) on a suitability map by first filling the most suitable areas, then the next most suitable and so on, until a pre-defined number of parcels have been placed on the map.

This is probably best explained with an example. I have a shapefile with polygons on it (counties) and each polygon has a value (number of acres for a crop). I also have a suitability map (0-9 suitability for that crop). I need to place the number of acres defined in each polygon onto the suitably map. I want to first fill the most suitable lands (9) and if there is still land left fill the next most suitable land (8) and so on until the number of acres for that county are distributed. Each polygon has a different number of acres and a different range of suitability, e.g. some will have 0-9, others might only have 3-7. Within each suitability layer the placing of the parcels can be random.

I have just started looking into Python scripting for ArcGIS. It seems like I might be able to do this with Python but I'm not really sure where to start. I would be interested in hearing any suggestion or advice on how to accomplish this.

share|improve this question
2  
you problem is a Land Consolidation map ? => en.wikipedia.org/wiki/Land_consolidation –  geogeek Jul 27 '12 at 20:37
    
when you say you want to "fill" the most suitable lands, do you mean accumulate a collection that fits certain cost requirements? if you're still trying to solve this problem, could you write a more specific explanation, perhaps with images? –  gotanuki Aug 28 '12 at 5:34
    
Well, putting it very simple, all you need to do is order your parcels by "suitability", and iterate over them summing their area until the sum becomes higher than the desired total value of acres. –  Alexandre Neto Jun 21 '13 at 21:45
    
An efficient solution would sort the parcel records by decreasing suitability, compute their cumulative areas until the desired total area is reached, and select all such parcels. This is a purely data summary operation, needing no GIS capabilities at all. Note, though, that it is incorrect to refer to such solutions as "random": they are perfectly determined. It is unclear how you intend to handle parcels that are assigned only ranges of suitability: perhaps that is why no answers have been offered. –  whuber Jul 30 '13 at 18:13
add comment

1 Answer 1

up vote 1 down vote accepted

Thanks for the comments. The solution we came up with is posted below.

# 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
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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