topic Re: arcpy, Spatial Analyst, and Exp in Python Questions

arcpy, Spatial Analyst, and Exp

JenniferHorsman — Sat, 11 Dec 2021 09:38:21 GMT

Hello,

I am trying to convert some 9.3 python scripts that use the geoprocessor object gp to use arcpy instead. I am struggling with the conversion of Single Output Map Algebra (SOMA) that uses exp. Also, the script first calculated IDW on a set of points and that alone was quite a struggle to figure out since the 10.2 help for Spatial Analyst IDW does not show the correct input parameters! Add to the confusion the decision of when to use arcpy.Raster and when not to.

Anyway, I am very close, but for some reason when I run arcpy.Exp_sa on my raster as such:

ExpOut = arcpy.Exp_sa(IDWOut)

another raster is automatically output by Exp with the same name as IDWOut, but with Exp_ at the beginning and _1 at the end. For example, if IDWOut is "IDW_PAH_filtered", then a raster called "Exp_IDW_PAH_filtered_1" is automatically created. According to the help page http://resources.arcgis.com/en/help/main/10.2/index.html#//009z0000008s000000 the output should be saved to ExpOut, but when I look at the values of that raster, I can see that they have not had the exponential applied.

So confused!!

Here are the actual code snippets:

geoDB = sys.argv[1]
COCLayer = sys.argv[2]
COCField = sys.argv[3]
PerfLog = sys.argv[4]
IDWCellSize = sys.argv[5]
IDWPower = sys.argv[6]
IDWRadius = sys.argv[7]
IDWMask = sys.argv[8]
IDWBarrier = sys.argv[9

and

IDWOut = geoDB + "\\IDW_" + COCLayer
arcpy.Idw_sa(geoDB + "\\temp_data", COCField, IDWOut, IDWCellSize, IDWPower, IDWRadius, IDWBarrier)
ExpOut = arcpy.Exp_sa(IDWOut)
IDWOut = arcpy.Raster(ExpOut) / 1000

And the temp_data is a temporary point shapefile that is created correctly by the script. I started wondering if Exp has the wrong input parameters on the help page just like Idw and maybe needs an output raster as one of the parameters, like:
arcpy.Exp_sa(IDWOut, ExpOut)
but when I tried that, I got an error saying "name 'ExpOut' is not defined".

If anyone has any insights, I will welcome them. Yes, I know I can still use the old gp method, but I really want to avoid that. By the way, here is what it used to look like:

IDWOut = geoDB + "\\IDW_temp"
outRaster = geoDB + "\\IDW_" + COCLayer


gp.Idw_sa(geoDB + "\\temp_data", COCField, IDWOut, IDWCellSize, IDWPower, IDWRadius, IDWBarrier)       
inExpress = "( exp (" + IDWOut + ") ) / 1000"
gp.SingleOutputMapAlgebra_sa(inExpress, outRaster)

The inputs sys.argv[] were the same. In that case, IDWOut was a temporary raster that was deleted later and the final output was outRaster, but doing something similar with arcpy was creating even more problems.

Re: arcpy, Spatial Analyst, and Exp

XanderBakker — Sat, 11 Dec 2021 09:38:24 GMT

Hi Jennifer,

There are a few suggestion I have when I look at your second snippet of code:

use os.path.join(pathname, filename) to create a valid output (requires "import os")

use "IDW_{0}".format(COCLayer) to create the name

access the Idw tool through arcpy.sa.Idw

note that the output is not a parameter of IDW

the result of a arcpy.sa tool is normally a raster object (in memory)

use the .save methode on the object to save the raster to a location

import arcpy, os

geoDB = r'C:\Path\to\Folder\With\GeoDB.gdb'
COCLayer = 'aLayerName'
COCField = 'aFieldName'
IDWCellSize = 10
IDWPower = 2
IDWRadius = 100
IDWBarrier = r'C:\Path\to\Folder\With\GeoDB.gdb\myBarriers'

# use os.path.join to combine path and file names and use string.format to format names
idw_out = os.path.join(geoDB, "IDW_{0}".format(COCLayer))

# Idw (in_point_features, z_field, {cell_size}, {power}, {search_radius}, {in_barrier_polyline_features})
fc_points = os.path.join(geoDB, 'temp_data')
myIDWraster = arcpy.sa.Idw(fc_points, COCField, IDWCellSize, IDWPower, IDWRadius, IDWBarrier)
myEXPraster = arcpy.sa.Exp(myIDWraster)
myFinalRaster = myEXPraster / 1000 # arcpy.Raster(myEXPraster) / 1000
myFinalRaster.save(idw_out)

Kind regards,

Xander

BTW see samples in Help:
Exp (Spatial Analyst)
IDW (Spatial Analyst)

Re: arcpy, Spatial Analyst, and Exp

JenniferHorsman — Sat, 11 Dec 2021 09:38:27 GMT

Hi Jennifer,

There are a few suggestion I have when I look at your second snippet of code:

use os.path.join(pathname, filename) to create a valid output (requires "import os")

use "IDW_{0}".format(COCLayer) to create the name

access the Idw tool through arcpy.sa.Idw

note that the output is not a parameter of IDW

the result of a arcpy.sa tool is normally a raster object (in memory)

use the .save methode on the object to save the raster to a location

Thank you for your reply. I was persistent yesterday and found a solution. By the way, as I mentioned in my original post, I did follow the suggestions on the help pages for Exp and Idw, and unfortunately, both appear to be incorrect for ArcGIS 10.x versions. Both of them can take an output file name as a parameter. And if an output file name is not supplied, then both will write output to the workspace (in my case, a personal geodatabase) with an automatically generated name.

Here is how I implemented Idw and Exp:

The script arguments (inputs) remain the same. I added another variable back in, a string for the name of the final result output raster, outRaster. I added parameters for the temporary outputs for Exp and Idw which I deleted later. And, I used another raster variable, MARaster, for some math algebra output. This step was probably not necessary, but it helped me sort out the debugging process.

IDWOut = geoDB + "\\IDW_temp"
outRaster = geoDB + "\\IDW_" + COCLayerBase


arcpy.Idw_sa(geoDB + "\\temp_data", COCField, IDWOut, IDWCellSize, IDWPower, IDWRadius, IDWBarrier)
arcpy.Exp_sa(IDWOut, ExpOut)
MARaster = arcpy.Raster(ExpOut) / 1000
MARaster.save(outRaster)


arcpy.Delete_management(geoDB + "\\temp_data")
arcpy.Delete_management(IDWOut)
arcpy.Delete_management(ExpOut)

Re: arcpy, Spatial Analyst, and Exp

XanderBakker — Sat, 11 Dec 2021 09:38:29 GMT

Hi Jennifer,

First of all, I'm glad that things work for you, but somehow there is something that doesn't make sense to me. The IDW interpolation is available in Spatial Analyst, 3D Analyst and Geostatistical Analyst. The IWD in ga is different, since it also requires a layer, so let's take this one out of the equation.

If you manually start the tool (in 10.2) and copy the Python snippet after a successful run, you can see the code generated. In case of 3D Analyst it call the tool from;

arcpy.Idw_3d

In case of Spatial Analyst it does this:

arcpy.gp.Idw_sa

Apparently it calls it from the old (?) gp module. That's odd... If you go the the Help page on Idw in Spatial Analyst and scroll down to the code sample you see they are using;

arcpy.sa.Idw

The 3D Analyst IDW code sample in the Help is using the same arcpy.Idw_3d. You are using arcpy.Idw_sa. In Python you can print the help from any function like this:

help(arcpy.sa.Idw)

It prints the Help on the function. Just have a closer look at what happens:

arcpy.Idw_3d
Help on function Idw in module arcpy.ddd:
Idw(in_point_features=None, z_field=None, out_raster=None, cell_size=None, power=None, search_radius=None, in_barrier_polyline_features=None)

arcpy.ddd.Idw
Help on function Idw in module arcpy.ddd:
Idw(in_point_features=None, z_field=None, out_raster=None, cell_size=None, power=None, search_radius=None, in_barrier_polyline_features=None)

arcpy.sa.Idw
Help on function Idw in module arcpy.sa.Functions:
Idw(in_point_features, z_field, cell_size='#', power='#', search_radius='#', in_barrier_polyline_features='#')

arcpy.gp.Idw_sa
Help on function <lambda> in module arcpy.geoprocessing._base:
<lambda> lambda *args

arcpy.Idw_sa
AttributeError: 'module' object has no attribute 'Idw_sa'

So, the 3D IDW does have an output raster parameter. The Idw_sa (inside gp or in directly from arcpy) do not show any help. However from the python snippet obtained from manually executing the Spatial Analyst IDW we know this:

arcpy.gp.Idw_sa("myPointFeatureClass","FieldName","myOutputRaster","cellSize","power","some other parameters","#")

It does have an output raster parameter.

Many tools have different ways to execute them:

#this:
arcpy.management.AddField
# is the same as:
arcpy.AddField_management

What you should take into account for future development is that within Spatial Analyst there is a trend to use the structure:

outputraster = tool(parameters)

This creates a temporal raster object which can be saved by calling the save methode on the object;

outputraster.save(outPathAndName)

So use arcpy.sa.Exp and arcpy.sa.Idw (instead of arcpy.Exp_sa and arcpy.Idw_sa).

It would be nice if Jason Scheirer (from the Python development team at Esri) could jump in the conversation and shed his light on this, but I suppose he is too busy with the DevSummit preps which is starting today...

Kind regards,

Xander