Crear un polígono de área fija basado en la superficie de coste/fricción

Question

Estoy utilizando ArcGIS Desktop 10.4 y tengo todas las extensiones. Tengo un shapefile de 494 puntos. También tengo un ráster de superficie de coste/superficie de fricción. Quiero crear un polígono de 13,4 hectáreas alrededor de cada punto basado en la superficie de fricción. En otras palabras, no quiero un simple buffer alrededor de cada punto. El polígono debe tener un aspecto irregular, ya que debe seguir los valores de fricción, así que, por ejemplo, si hay una colina empinada al este, el valor de fricción será mayor.

También quiero producir otro polígono con los mismos puntos y la superficie de fricción que contorneará una hora de camino desde el punto.

Answer 1

Lo he comprobado y no estoy seguro de que esto se pueda hacer en el modelo, porque necesita el cálculo del total acumulado en la tabla. De todas formas el script de abajo resuelve la primera parte de tu pregunta.

import arcpy,scipy
from arcpy.sa import *
import pandas as pd
from scipy import interpolate

dem = arcpy.Raster("DEM")
vf = VfTable(r'c:\SCRATCH\vf.txt')
aStep=5
target=200000
mxd = arcpy.mapping.MapDocument("CURRENT")
lyr = arcpy.mapping.ListLayers(mxd,"POINTS")[0]
g=arcpy.Geometry()

PGONS=[]
tbl = arcpy.da.TableToNumPyArray(lyr,"OID@")
for i,fid in enumerate(tbl):
    arcpy.AddMessage('%i out of %i' %(i+1,len(tbl)))
    lyr.setSelectionSet ("NEW",list(fid))
    outPD = PathDistance(lyr, "", "", "", "", dem, vf)
    outPD.save("C:/scratch/pd")
    intR=Int(outPD/aStep)
    vrt=arcpy.management.BuildRasterAttributeTable(intR, "Overwrite")
    areas = arcpy.da.TableToNumPyArray(vrt,("Value","COUNT"))
    df=pd.DataFrame(areas)
    cumSum=df["COUNT"].cumsum()
    f = interpolate.interp1d(cumSum,range(len(df)))
    threshold = f(target)*aStep
    pGon=Con(outPD<=threshold,1)
    areas=arcpy.conversion.RasterToPolygon(pGon, "in_memory/pgon", "NO_SIMPLIFY")
    lp=0
    gList=arcpy.CopyFeatures_management(areas,g)
    for item in gList:
        if item.area<lp:continue
        lp=item.area
        biggie=item
    PGONS.append(biggie)
arcpy.CopyFeatures_management(PGONS,"C:/scratch/pgons.shp")

Flujo de trabajo:

• seleccione el primer punto y utilícelo para calcular la distancia de la trayectoria utilizando Función de senderismo de Toblers
• Convertir la distancia en intervalos finitos
• Calcula el área acumulada en la trama anterior, interpola la curva (CumArea/Valor) en el área objetivo (200 000 celdas en el script anterior) => umbral.
• Extraer el área por debajo del umbral del ráster de distancia de la ruta y convertirlo en polígono.
• Proceda con el siguiente punto.

Nota: los polígonos están etiquetados por su superficie. Como se puede ver, son un poco menos de 200k celdas.

Cualquier duda, pregunta antes de que se cierre el post.

UPDATE, usando numpy, podría no funcionar con DEM muy grandes. Produce ambos conjuntos de polígonos.

import arcpy
from arcpy.sa import *
import numpy as np
''' PARAMETERS '''
targetHa = 20
minutes = 7
''' ONE OFF OPERATIONS '''
dem = arcpy.Raster("DEM")
target = targetHa*10000/dem.meanCellHeight**2
seconds = minutes*60
vf = VfTable(r'c:\SCRATCH\vf.txt')
mxd = arcpy.mapping.MapDocument("CURRENT")
lyr = arcpy.mapping.ListLayers(mxd,"POINTS")[0]
g=arcpy.Geometry()
areaPGONS=[];timePGONS=[]
tbl = arcpy.da.TableToNumPyArray(lyr,"OID@")

def queryRaster(R,V,appendTo):
    pGon = Con(R <= V, 1)
    areas=arcpy.conversion.RasterToPolygon(pGon, g, "NO_SIMPLIFY")
    gList=arcpy.CopyFeatures_management(areas,g)
    appendTo+=gList
    return appendTo
''' SHUFFLE THROUGH POINTS '''
for i,fid in enumerate(tbl,1):
    arcpy.AddMessage('Processing %i out of %i' %(i,len(tbl)))
    lyr.setSelectionSet ("NEW",list(fid))
    outPD = PathDistance(lyr, "", "", "", "", dem, vf)
    np2d = arcpy.RasterToNumPyArray(outPD,"","","",np.nan)
    threshold = float(np.sort(np2d, axis=None)[target]); del np2d
    areaPGONS = queryRaster(outPD,threshold,areaPGONS)
    timePGONS = queryRaster(outPD,seconds,timePGONS)
##  REMOVE LINE BELOW AFTER TEST RUN
    if i == 3:break
arcpy.CopyFeatures_management(areaPGONS,"C:/scratch/area_pgons.shp")
arcpy.CopyFeatures_management(timePGONS,"C:/scratch/time_pgons.shp")
arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")
arcpy.AddMessage('\nRemove polygons equal cell size from both sets\n')

Obsérvese la mayor distancia que pueden recorrer los humanos en un terreno más llano y la mayor precisión de las estimaciones de los polígonos de 20 ha. La línea recta de la derecha es un borde del DEM real que he utilizado.

También hay que tener en cuenta que la sugerencia de @Vince podría funcionar para los polígonos temporales. Se pueden calcular de una sola vez, porque la herramienta de Asignación tiene una opción de distancia máxima. Sólo hay que ponerla en 3600 (segundos). Sin embargo, si tiene grupos de puntos más cercanos entre sí que 2 horas de camino, obtendrá un solo polígono de tiempo por grupo.