5 votos

¿Cómo colocar automáticamente paquetes en el mapa de adecuación?

Necesito hacer un mapa en el que se coloca al azar en parcelas (grupos de píxeles o polígonos) en un mapa de idoneidad por primera llenado de las áreas más adecuadas, entonces el siguiente más adecuado y así sucesivamente, hasta que un pre-definida número de parcelas que se han colocado en el mapa.

Este es probablemente el mejor explicado con un ejemplo. Tengo un shapefile con polígonos (condados) y cada polígono tiene un valor (número de hectáreas de cultivo). También tengo un mapa de idoneidad (0-9 idoneidad para el cultivo). Necesito colocar el número de hectáreas definidas en cada polígono adecuadamente en el mapa. Primero quiero llenar las tierras más aptas (9) y si todavía hay tierra a la izquierda rellenar el siguiente más adecuado de la tierra (8) y así sucesivamente hasta que el número de hectáreas para que el condado se distribuyen. Cada polígono tiene un número diferente de hectáreas y un rango diferente de idoneidad, por ejemplo, algunos se han 0-9, otros sólo pueden tener de 3 a 7. Dentro de cada idoneidad de la capa de la colocación de las parcelas de que se puede ser al azar.

He empezado a buscar en secuencias de comandos de Python para ArcGIS. Parece que podría ser capaz de hacer esto con Python, pero no estoy muy seguro de por dónde empezar. Yo estaría interesado en escuchar cualquier sugerencia o consejo sobre cómo lograr esto.

1voto

Ankit Ghate Puntos 86

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
 

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X