Hi,
I am using arcgis 10.5.1 64 bit.
I have a set of 10 point features that I am trying to connect using PathDistance. Each point has an id, say 0...9. Effectively I want to generate the path from 0 to 1, 1 to 2, and so on.
I am using both my own tables to define vertical and horizontal factors. I have created these paths manually using the exact same factors. So I know those work.
The two main files I am using are: 1) a sample DEM (5m resolution, tif) called test_dem.tif 2) a point shapefile test_loc.shp which contains 11 points with several columns but the main one used here is 'Id' column. These files and others used for the script can be found below.
The files used for the vertical and horizontal factors are the following,
hsine.txt | vertical3.txt |
---|---|
0 1 | -90 20000 -80 5000 -70 1500 -60 300 -50 100 -40 40 -30 12 -20 3 -10 1 0 3 10 7 20 14 30 24 40 40 50 75 60 180 70 750 80 1500 90 5000 |
I have the following code:
import arcpy
import arcpy.sa as sa
import arcpy.management as mng
from math import pi
# Checking out SA extension
class SAerror(Exception):
pass
try:
if arcpy.CheckExtension('Spatial') == 'Available':
arcpy.CheckOutExtension('Spatial')
else:
raise SAerror
except SAerror:
print('Spatial Analyst license is unavailable!')
except arcpy.ExecuteError:
print(arcpy.GetMessages(2))
# ArcGIS directories
arcpy.env.scratchWorkspace = r'C:\GIS\projects\network\potential\SCRATCH'
arcpy.env.workspace = r'C:\GIS\projects\network\potential'
# input rasters
rDEM = arcpy.Raster(r'C:\GIS\projects\network\potential\tst_dem.tif')
shpFeat = r'C:\GIS\projects\network\potential\testloc.shp'
# PARAMETERS
deg2rad = pi/180.0
rad2deg = 180.0/pi
wFactor = 0.2
I = 4.0
wPthPot = 0.75
Gmax = 5.0
# path distance parameters
rCost = "" # use unit
rSurface = rDEM
rVertical = rDEM
hCostfn = r'C:\GIS\projects\network\potential\hsine.txt'
HF = sa.HfTable(hCostfn)
vCostfn = r'C:\GIS\projects\network\potential\vertical3.txt'
VF = sa.VfTable(vCostfn)
maxDist = ""
nbrfn = r'C:\GIS\projects\network\potential\expd.txt'
#Set geoprocessing environmental preferences
arcpy.env.extent = rDEM.extent
arcpy.env.cellSize = rDEM.meanCellHeight
arcpy.env.overwriteOutput= True
# INITIAL CONDITIONS
rTPot = rDEM / rDEM.maximum
rDirRel = rTPot
rG0 = sa.CreateConstantRaster(0.0, 'FLOAT', arcpy.env.cellSize, arcpy.env.extent)
rGmax = sa.CreateConstantRaster(Gmax, 'FLOAT', arcpy.env.cellSize, arcpy.env.extent)
rGt_1 = sa.CreateConstantRaster(0.0, 'FLOAT', arcpy.env.cellSize, arcpy.env.extent)
# sample sequence
pthseq = list(range(11))
# create a feature layer
mng.MakeFeatureLayer(shpFeat, 'lyr')
# MAIN
for o,d in zip(pthseq[:-1], pthseq[1:]):
# lambda direction
rDir = sa.Aspect(rDirRel)
# CALCULATE PATH
mng.SelectLayerByAttribute('lyr', "NEW_SELECTION", ' "Id" = '+str(o)) # origin
# accumulated cost surface
rAcs= sa.PathDistance('lyr', "", rSurface, rDir, HF, rVertical, VF, "", 'rBL.tif')
mng.SelectLayerByAttribute('lyr', "NEW_SELECTION", ' "Id" = '+str(d)) # destination
# path
rPath = sa.CostPath('lyr', rAcs, 'rBL.tif', "EACH_CELL")
rPath = sa.Con(sa.IsNull(rPath), 0, 1)
pathFn = 'path'+str(o)+'_'+str(d)+'.tif'
rPath.save(pathFn)
rW = (rG0 - rGt_1) * wFactor
rMark = rPath * (1.0 - (rGt_1 / rGmax)) * I
rGt = rGt_1 + rW + rMark
rGt.save('rGt.tif')
# DECAY
rPthPot = sa.FocalStatistics(rGt, sa.NbrWeight(nbrfn),'SUM')
rPthPot = rPthPot.maximum - rPthPot
rPthPot = rPthPot / rPthPot.maximum
# Update directional relation
rDirRel = rPthPot * wPthPot + rTPot * (1.0 - wPthPot)
# Update rGt_1 *** missing earlier ***
rGt_1 = rGt
Overall, I get very weird behavior. Some times, it will fail when trying to generate the backlink. In other occasions, it will calculate the first path and then fail. It will calculate three paths and fail. Only once it has ran through the entire sequence. After crashing I make sure to delete all possible files (usually the backlink file is created but it is corrupted).
I have tried many things:
1. Creating a fresh feature layer and erasing it at the end of each iteration.
2. Generating with each iteration different feature layers, one for the origin and another for the destination, and deleting them after each iteration.
3. I have tried saving to the disk the layer feature as well as all the other rasters (I set arcpy to overwrite).
4. I have tried saving to the disk the backlink
5. I have used relative and absolute filenames.
6. I am using only .tif files
7. I have use try...except to catch information about the error with no success
I have spent much time reading through the different forums, etc. I know there were problems with this method in previous iterations of the software. The lack of any feedback relating to possible errors is a real problem! Often it will be python that crashes all together. It is very finicky, there are certain Horizontal Factor settings (even those found in the help pages) that will not generate a complete cost accumulation raster. This is not my case though.
My question is whether there is something that I am missing, something that needs to be reset or deleted? or perhaps someone can suggest how to find more information after the procedure fails???
thank you,
Marc
Having just done this myself...
I pointed PyCharm to the python.exe in here :
C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3
Thanks Neil,... I did find out at the end (only two options). Still running
into a non-descrip error with no clue as to why things do not work.
Marc
Xander Bakker Dan Patterson Incidentally perusing the release notes for ArcGis pro 2.0 I found the following under Issues Addressed:
NIM099397 | The Path Distance tool generates an output raster with cell values of all zero when the Table option is used for the Vertical Factor parameter. |
Not certain this might be the cause of failure. Thought I would mention it...
Thanks mllo1 for mentioning. Not sure if this explains everything, but it certainly does not help the get the result you are after.
I did manage to run the script in arcgis pro 2.0. It failed once again in the Path Distance but it did not blow up. I managed to get this output,
"C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\python.exe" C:/Users/Marcos/PycharmProjects/Network/network.py
Traceback (most recent call last):
File "C:/Users/Marcos/PycharmProjects/Network/network.py", line 92, in <module>
rAcs = sa.PathDistance('lyr', "", rSurface, rDir, HF, rVertical, VF, "", "rBL.tif")
File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\sa\Functions.py", line 1314, in PathDistance
source_direction)
File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
result = wrapper(*args, **kwargs)
File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\sa\Functions.py", line 1298, in Wrapper
source_direction)
File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\_base.py", line 506, in <lambda>
return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError: ERROR 999999: Error executing function.
Failed to create raster dataset
Failed to execute (PathDistance).
Still somewhat cryptic! 😞
Yep, ERROR 999999 are the best there are... Good thing it didn't blow up, so no need to include a message to wear protective eye-wear when running this script.
Hi xander_bakker and Dan_Patterson,
Still have not resolved this...
When I run the script (as is) using arcgis pro 2.01, it stops the first time it encounters PathDistance with an error. If, however, I run it through the desktop application (with the exact same parameters) it works fine. Can any of you suggest ways in which I could explore this behavior? By that I mean... what are the possible differences when running it as a script or through the application? Perhaps there lies the problems....
thanks,
Marc
I guess in that case it would be best to contact support and see if they know why the error occurs.
Edit: in case you have a tool that executes correctly in Pro and when you run script it does not, you should try to isolate the problem. When you execute a tool you can copy the Python snippet from the Results window. In case you execute that snippet does that work (generating a new result)?