Create polygon of fixed area based on cost/friction surface [closed]

Question

I am using ArcGIS Desktop 10.4 and have all extensions. I have a shapefile of 494 points. I also have a cost surface/friction surface raster. I want to create a polygon of 13.4 hectares around each point based on the friction surface. In other words I don't want a simple buffer around each point, The polygon should look irregular as it should follow the friction values, so for example if there is a steep hill to the east the friction value will be higher.

I also want to produce another polygon with the same points and friction surface that will outline one hours walk from the point.

Answer 1

I've checked and not sure anymore this is doable in model, because it needs computation of cumulative total in the table. Anyway script below solves 1st part of your question.

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")

Workflow:

• select 1st point and use it to compute path distance using Toblers' hiking function
• Convert distance to finite intervals
• Compute cumulative area in above raster, interpolates curve (CumArea/Value) at target area (200 000 cells in above script) => threshold.
• Extract area below threshold from path distance raster and convert it to polygon.
• proceed with next point.

Note: polygons labelled by their area. As one can see they are wee bit short of 200k cells.

Any questions, ask before post get closed.

UPDATE, using numpy, might not work with very large DEM. Produces both sets of polygons.

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')

Note much larger distance humans can walk on a flatter terrain and greater accuracy of 20 ha polygon estimates. Straight line on the right is an edge of actual DEM I've used.

Also note, that @Vince suggestion might work for time polygons. They can be computed in one go, because Allocation tool has an option of maximum distance. Just set it to 3600 (seconds). However if you have clusters of points closer to each other than 2 hours walk, you'll get a single time polygon per cluster.