Kernel Interpolation with Barrier - Python Issue (bug?)

[BryanNicholson]

Hello all, I have perhaps an odd issue with regards to the Kernel Interpolation with Barriers (KIB) prediction method using Python (2.7) and ArcGIS 10.1.

First, using the Geostatistical Analyst wizard to perform the KIB, or Kernel Smoothing as it is referred to in the wizard, works fine and dandy. However, I want to move the KIB from the wizard into a python script so I can automate it for ~100 maps. 

First, here are some parameters I use and want to apply from the wizard: Function=Exponential; Surface Type=Prediction; Bandwidth=1000; Smoothing Factor=0.2; and Radius=8000.

The link for the method is here: Kernel Interpolation With Barriers—Help | ArcGIS Desktop 

The issue is that the KIB does not have a place to put the Smoothing and Radius (Search Neighborhood) information into the python script. I tried to run the script without it but I get this error:

ERROR 040008: Input value is out of range.
Failed to execute (KernelInterpolationWithBarriers).

I also tried to run the tool from within ArcMap but it also does not have a spot to input the search neighborhood. 

I tried to figure out a way to include the Search Neighborhood environment in the script but I haven't found a way.

Does anyone have a solution to this or any feedback?

Thank you!

[DarrenWiens2]

I can't help you with the fact that there is no parameter for search neighbourhood, but if you post your script leading up to the error, we may be able to figure that out.

[EricKrause]

Hi Bryan,

You are correct that the Kernel Interpolation with Barriers geoprocessing tool does not have a search neighborhood option.  This option is only available in the Geostatistical Wizard and only for Exponential, Gaussian, and Constant kernel types.

However, you can automate your process using the Create Geostatistical Layer geoprocessing tool.  This tool takes a model source as input and allows you to apply that model source to a new dataset.  To create the model source, use the Geostatistical Wizard and perform Kernel Smoothing with the parameters that you want (it doesn't matter which dataset you use).  When you click Finish, there will be a "Save..." option on the Method Report dialog.  This will allow you to save an xml file that will serve as the model source in Create Geostatistical Layer.

In your Python script, run Create Geostatistical Layer for each of your datasets, and use the xml from before as the model source for all of them.  Since you are using ArcMap 10.1, you will need to provide your datasets and fields as strings (more recent versions can use an arcpy class).  See the script examples in the tool documentation above, and feel free to post more questions if you run into any problems with the script.

[BryanNicholson]

Thanks for the reply, I think this is what I needed. I'm just having an issue on the implementation. I'll post some code here for you to look at. I've never pasted code here before so I hope I'm doing this right.

First I am getting this error:

Failed to execute. Parameters are not valid.
ERROR 000840: The value is not a Geostatistical Layer.
ERROR 000814: Invalid file type
ERROR 045001: Input dataset(s) error. Invalid number of datasets. Model requires 2 inputs vs provided 1
Failed to execute (GACreateGeostatisticalLayer).

Below is a snippet of my code:

arcpy.CheckOutExtension("GeoStats")

# Is this supposed to be the .xml file? 

KS_Model = "C:/Arc_Workspace/Metals_Workspace.gdb/KernelS_Template.xml"

# I think my main issue is here, the documentation isn't clear on what the X=Shape, F1='OZONE' , etc. is supposed to be in my case. I think I need to leave the zField variable in here somehow so that it knows which field to include in the model for each iteration. BH_Points_export is a feature class, btw, should it be a shp?
inData = "C:/Arc_Workspace/Metals_Workspace.gdb/BH_Points_export " + zField

# zField is the field.aliasname of the field that has the z values for the Kernel prediction
outGL = str(zField)      #str() to make sure it is a string input

OutKSlyr = arcpy.GACreateGeostatisticalLayer_ga(KS_Model, inData, outGL)

# I've also tried using this, but I get the same error:

OutKSlyr = arcpy.GACreateGeostatisticalLayer_ga("C:/Arc_Workspace/Metals_Workspace.gdb/KernelS_Template.xml",
              "C:/Arc_Workspace/Metals_Workspace.gdb/BH_Points_export F1=FieldName",
              "outGL")

Any suggestions?

Thanks again!

[EricKrause]

Hi again,

Your syntax should be something like:

arcpy.GACreateGeostatisticalLayer_ga("C:/temp/myModelSource.xml", "C:/temp/myDataset.shp myField", "myNewGALayer")

The first parameter should be a file path to the xml.  The second parameter will be a string where you first specify the file path to the dataset (shp and fgdb both work) followed by a space, then the name of the field.  The third parameter is the name of the new geostatistical layer that is created by applying the parameter from the model source to the new dataset.

Since geostatistical layers are saved in memory, you probably also want to save them as a layer file on disk so that you can easily open them in maps.  After the code above, use something like:

arcpy.SaveToLayerFile_management("myNewGALayer","C:/temp/galayerfile.lyr")

Regarding the X, Y, and F1 options that appear in the example, these options are not required in your case.  They are used when the X and Y coordinates of the points are stored in fields rather than in the SHAPE attribute of the feature class. If your datasets were, for example, csv files (which don't have a SHAPE attribute), you could pass the X coordinates, the Y coordinates, and the field (F1) using those options.

One final question... are you using any barriers in your interpolation?  If so, I'll need to give a slightly different code example of how to do that.  If you aren't using barriers, the code above should work.

[BryanNicholson]

I am using a barrier actually. The barrier is a polyline that connects to

itself wherein all the points are located and interpolation is meant to

occur.

Thanks!

[EricKrause]

Thanks for the clarification.

You'll need to alter the statement to something like:

arcpy.GACreateGeostatisticalLayer_ga("C:/temp/myModelSource.xml", "C:/temp/myDataset.shp myField ; C:/temp/myBarriers.shp", "myNewGALayer")

The difference is that you need to pass the feature class, the field, and the barriers in a single string.  You separate the two datasets with a semicolon like in the example above.

[BryanNicholson]

Okay, so I have this:

fields = arcpy.ListFields(inFeature)
for field in fields:
      FieldName = field.aliasName
      if TwoAcid in FieldName:      #TwoAcid is just a string I'm scanning for in each field in the table

            zField = str(FieldName)

            arcpy.GACreateGeostatisticalLayer_ga("C:/Arc_Workspace/Metals_Workspace.gdb/KernelS_Template.xml",
                  "C:/Arc_Workspace/Metals_Workspace.gdb/BH_Points_export field;

                   C:/Arc_Workspace/Metals_Workspace.gdb/Water_Boundary",
                  "outGL")

...and I get this error:

Failed to execute. Parameters are not valid.
ERROR 000840: The value is not a Geostatistical Layer.
ERROR 000814: Invalid file type
ERROR 045001: Input dataset(s) error. `Dataset`: Field `field` does not exist
Failed to execute (GACreateGeostatisticalLayer).

I get this same error message above when using zField in place of 'field'.

I also tried placing the anticipated field name in the 'field' spot, to which I get this error: 

Failed to execute. Parameters are not valid.
ERROR 000840: The value is not a Geostatistical Layer.
ERROR 000814: Invalid file type

Am I still missing something in my code, should I convert the feature class's into shapefiles?

-Bryan

[BryanNicholson]

Small update. I used the Feature class to Shapefile tool to convert the input feature class and the barrier into shapefiles. After running the script again I continue to get the same errors related to a non-existent Field.

[EricKrause · Accepted solution]

I think you have a couple problems.  First, use the field name, not the field alias.

Change this line: FieldName = field.aliasName

to

FieldName = field.name

Later in the code, you need to pass the variable into the tool with:

            arcpy.GACreateGeostatisticalLayer_ga("C:/Arc_Workspace/Metals_Workspace.gdb/KernelS_Template.xml",
                  "C:/Arc_Workspace/Metals_Workspace.gdb/BH_Points_export " + FieldName + ";

                   C:/Arc_Workspace/Metals_Workspace.gdb/Water_Boundary",
                  "outGL")

If you pass "field" inside the string of the parameter, it won't be recognized as a variable.  It will look for a field called "field" and fail when it doesn't find one.

You will also want to store the output geostatistical layer (the last parameter) as a variable.  If you keep the same name ("outGL") for each iteration, it will just get overwritten after each iteration.  Or you can save each one as a layer file, but those will also need unique names.

I am not the greatest Python programmer, and I can't directly test this without your data, but I think what I just outlined will work. 

Optionally, you may want to also consider storing all of the parameters as variables in order to make it easier to read and understand.  Something like:

inXML = "C:/Arc_Workspace/Metals_Workspace.gdb/KernelS_Template.xml"

inPoints = "C:/Arc_Workspace/Metals_Workspace.gdb/BH_Points_export"

infield = "<your field>"

inBarrier = "C:/Arc_Workspace/Metals_Workspace.gdb/Water_Boundary"

inDataset = inPoints + " " + infield + ";" + inBarrier

outGAlayer = "outGL"

arcpy.GACreateGeostatisticalLayer_ga(inXML, inDataset, outGALayer)