Path Distance erratic and crashes in a loop

MLLO · ‎10-25-2017

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 10 1.02 20 1.06 30 1.13 40 1.23 50 1.36 60 1.5 70 1.66 80 1.83 90 2 100 2.17 110 2.34 120 2.5 130 2.64 140 2.77 150 2.87 160 2.94 170 2.98 180 3	-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

DanPatterson_Retired · ‎10-26-2017

Interesting... how can the raster contain all Nodata cells... it is a complete barrier to movement. Before this was discovered, were any optional parameters specified?

Perhaps what is needed is a description of how the cost surface was constructed, how slope and/or aspect fits into the picture, an origin and destination point sample, then run a simple Cost Distance on that prior to moving on to Path Distance (see how path distance works discussion) since there is no mention of even using horizontal or vertical factors.

MLLO · ‎10-27-2017

Thanks Xander Bakker,

The error you mention with the focal stats is the first I hear about it. I had not had that occur to me before. I double checked that file in notepad (I wrote a program to generate it) so I am mystified as to why it might be problematic (it seems to follow esri specs).

The script is meant to simulate primarily the effect that an existing path may have on later paths. This effect is accumulative (meaning that previous paths may affect later paths. The bit of code that you include, refers to different aspects of the simulation. rW is a kind of a 'back to initial condition' factor. With each new path, those locations that are not marked again (i.e. cells where the new path does not cross) with eventually restore to their initial condition (g0... zero in this case) the speed at which they do this will depend on the number of paths to be generated (in this case 0.2,... after 5 paths). rMark describes the impact (I) a new path has. This impact, however, is not unbounded, hence, it cannot exceed a certain maximum amount (Gmax). rGt describes the combined effect that all paths may have. It is updated by the presence of a new path. Keep in mind that after a new path is created, the path is reclassified into (1-path, 0-no path), i.e. there are no NoData cells.

Because of some of your comments I did look over the code and realized that I mistakenly omitted one important line (probably because I was trying to figure out some earlier errors). I have now included this at the very end of the script (I also made a comment). This may help to clarify a bit more the simulation.

As for the comments by Dan Patterson. Not certain I understand your comments. There is no cost surface because PathDistance does not require any (by default it uses 1, see How the path distance tools work—Help | ArcGIS Desktop ). As for the vertical and horizontal factors, I have included both files in the zip file and in the explanation (a simple plot XY shows how these are defined). Unfortunately, I cannot use CostDistance as it does not allow you to specify 'in_horizontal_raster' (see Path Distance—Help | ArcGIS Desktop ) .

DanPatterson_Retired · ‎10-27-2017

I am very familiar with how it works, which is why I was wondering what the raster of nodata's were. I perhaps erroneously assumed that you were including other 'costs' into the analysis other than those associated with a dem... my bad

MLLO · ‎10-27-2017

thanks Dan_Patterson‌ the costs are actually generated via the vertical factor (ie. they are simply a function of the slope)

XanderBakker · ‎10-27-2017

HI mllo1 ,

The error you mention with the focal stats is the first I hear about it. I had not had that occur to me before. I double checked that file in notepad (I wrote a program to generate it) so I am mystified as to why it might be problematic (it seems to follow esri specs).

The kernel you shared looks perfectly normal, so maybe the error was due to my version 10.3 that I have at home. Will try later using 10.5, or maybe change something in my regional settings that might be causing it. Thanks for the functional explanation.

MLLO · ‎10-31-2017

@Xander and Dan_Patterson‌ Did any of you manage to get this to work?

XanderBakker · ‎10-31-2017

Hate to tell you, but I did a test on 10.5 and running it from PyScripter crashed my IDE and running it from the Python Window in ArcMap was the fastest way to close my ArcMap session (crashed too). I haven't tried with ArcGIS Pro, but for some reason I haven't had any luck with it. I would probably have to refactor to get anywhere.

MLLO · ‎10-31-2017

Hi Xander Bakker,

Nothing I have not been suffering for the past week and half! I was hoping it was something dumb that I had just forgot and some one would pointed out. I have tried every trick I know with no success. Python just blows up. I am just puzzled by the behavior.

Any recommendations has to how to narrow down possible causes?

Thanks for trying!

DanPatterson_Retired · ‎10-31-2017

If it is python blowing up then that leave running it within ArcGIS Pro using Python 3.5.. although I see nothing specific in the python or the spatial analysis methods used that would suggest this.

MLLO · ‎11-01-2017

Hi Dan Patterson‌,

I am trying to switch to arcgis pro in order to test the script. I am somewhat confused with what I need to do... I am familiar with jupyter notebooks and pycharm. It seems that arcgis pro contains two python installations? In pycharm I pointed the interpreter to go to

C:\Program Files\ArcGIS\Pro\bin\Python\pkgs\python-3.5.4-0

but when I try importing arcpy it did not find it.

I just want to see whether the script works in arcgis pro before abandoning this project.