1

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.

  • 2
    What have you tried? There are multiple solutions for each point (there isnt just one shape that is the solution), what shape do you want? For example a limit on minimal roundness/compactness etc. – BERA Apr 15 at 17:11
  • 2
    It sounds as if CostAllocation would be involved in your solution. – Vince Apr 15 at 18:17
  • 1
    Repeat your cost distance calculation 497 times, I'd use Python, but doable in model. – FelixIP Apr 15 at 20:14
2

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.

enter image description here

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.

enter image description here

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.

| improve this answer | |
  • Hi Felix thats really helpful thanks. I have never used any scripting before but Im sure I can figure out how to use your script. Just for my clarity the only variables used are dem - which is my dem, and Points - which are my points ? Im also not sure where the target=200000 comes from. – user131981 Apr 16 at 11:23
  • Also Felix based on your answer if I can get it working I will probably use it in my PhD, studying settlement in an Early medieval kingdom in Ireland. I would like to cite you in the thesis if thats ok as it is your script I will be using - Im not sure how to get contact details etc on here, or if thats ok with you ? – user131981 Apr 16 at 11:37
  • My Dem cell size is 1m*1m. 200k is target number of cell, i.e. 20 ha. Run this script from mxd. Just call your points "points" in table of content. – FelixIP Apr 16 at 19:44
  • Thanks Felix I have emailed you now. – user131981 Apr 17 at 11:32

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