topic Re: Iterating Kriging with Anisotropy in ArcGIS GeoStatistical Analyst Questions

Iterating Kriging with Anisotropy

CourtneySpencer — Fri, 12 May 2017 17:51:29 GMT

Hi,

I'm trying to write a python script that will run kriging interpolation with anisotropy on 700+ variables in a dataset. I've read a lot already about using the Create Geostatistical Layer tool and was able to code it but cannot get it to automatically iterate through all my fields. It seems like I'd have to manually change the field each time which is just not efficient with the number of columns I have.

I've done an iteration before for a simple IDW interpolation using ListFields but it doesn't seem to work with this interpolation method.

Can anyone shed some light on this? I've pasted my code below. I'm VERY new to python so I apologize if it's a mess.

import os
import arcpy
import arcpy, os
from arcpy import env
from arcpy.sa import *
Sediment1999 = "C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\StableIsotopeData.gdb\\Sediment1999"
arcpy.env.workspace = r"C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\StableIsotopeData.gdb"
arcpy.env.extent = "C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\DetroitRiverBoundary.shp"
arcpy.env.mask = "C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\DetroitRiverBoundary.shp"

arcpy.CheckOutExtension("GeoStats")
fieldList = arcpy.ListFields("Sediment1999", field_type="Double")

fieldNames = [f.name for f in fieldList ]

for name in fieldNames :
    Kriging45 = arcpy.GeostatisticalDatasets("C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\Kriging45.lyr")
    Kriging45.Sediment1999 = "C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\StableIsotopeData.gdb\\Sediment1999"
    Kriging45.Sediment1999name = name
    arcpy.GACreateGeostatisticalLayer_ga("C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\Kriging45.lyr", "C:\\Users\\courtney\\Desktop\\StableIsotopeData.gbd\\Sediment1999 name", "GaLayer"+name)
    arcpy.GALayerToRasters_ga("GALayer"+name, (os.path.join("C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\Anisotropic\\Sediment\\1999", name+".tif")), "PREDICTION")

arcpy.CheckInExtension("GeoStats")

EDIT: Below is the error message I get.

Traceback (most recent call last):
File "C:\Python27\ArcGIS10.4\Lib\site-packages\Pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript
    exec codeObject in __main__.__dict__
File "E:\Work\Miscellaneous\Mundle Stable Isotope\KrigingCodeYearSpecific.py", line 20, in <module>
    arcpy.GACreateGeostatisticalLayer_ga("C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\Kriging45.lyr", "C:\\Users\\courtney\\Desktop\\StableIsotopeData.gbd\\Sediment1999 name", "GaLayer"+name)
File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\ga.py", line 1297, in GACreateGeostatisticalLayer
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 045001: Input dataset(s) error. Table of inputs is not complete.
Failed to execute (GACreateGeostatisticalLayer).

Thanks so much,

Courtney

Re: Iterating Kriging with Anisotropy

DanPatterson_Retired — Fri, 12 May 2017 18:14:55 GMT

you might try cleaning up the script and used the named variables approach in the last example of code in the help topic http://desktop.arcgis.com/en/arcmap/latest/tools/geostatistical-analyst-toolbox/create-geostatistical-layer.htm

just to make sure that you have the right inputs.

Re: Iterating Kriging with Anisotropy

SteveLynch — Fri, 12 May 2017 18:25:26 GMT

How did you create "C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\Kriging45.lyr". ?

Spatial Analyst's kriging output is a plain raster layer and should not be confused with a geostatistical layer.

You have to create a Geostatistical (kriging) layer by using the Geostatistical Wizard.

-Steve

Re: Iterating Kriging with Anisotropy

CourtneySpencer — Fri, 12 May 2017 18:29:49 GMT

I used the Geostatistical Wizard and ran the kriging with anisotropy first to get a result then saved it to a layer file. I initially had an .xml file and tried this method but when it didn't work I tried the layer file instead. Still didn't work.

I've read about using the intitial .xml or . lyr as a template then just re-running the Create Geostatistical Layer tool on each of my variables.

Re: Iterating Kriging with Anisotropy

SteveLynch — Fri, 12 May 2017 18:54:45 GMT

You have Kriging45 = arcpy.GeostatisticalDatasets(ga_layer) now you want to rather use another field, so

Kriging45.dataset1Field = whatever this field is called. dataset1Field is the property that you want to set (or get).

and likewise

Kriging45.dataset1 = name of the featureclass

Re: Iterating Kriging with Anisotropy

CourtneySpencer — Mon, 15 May 2017 13:13:20 GMT

Thanks Steve.

I understand that I can manually change the variable name in the tool parameters to interpolate each time but I was wondering if there was a way to automatically loop through all my variables without having to change the parameters each time? I have over 700 variables per dataset (of which I have about 10). That'll just take too long.

Re: Iterating Kriging with Anisotropy

CourtneySpencer — Tue, 08 Aug 2017 20:02:16 GMT

I just thought I would update my original question with the solution my colleague and I discovered in case anyone else has the same problem.

# Name: FinalKrigingCode.py
# Description: Uses an existing Geostatistical Layer (xml file) as a template
#              to interpolate another variable, then exports the new Geostatistical
#              Layer as a raster layer. Uses the ListFields function to iterate through
#              all fields within the shapefile.
# Requirements: Geostatistical Analyst Extension

# Import system modules
import arcpy, os

# Set environment settings
arcpy.env.workspace = "C:\\Users\\courtney\\Desktop\\MundleStableIsotope"
arcpy.env.extent = "C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\DetroitRiverBoundary.shp"
arcpy.env.mask = "C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\DetroitRiverBoundary.shp"
Sediment1999 = "C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\Sediment1999.shp"

# Check out the ArcGIS Geostatistical Analyst extension license
arcpy.CheckOutExtension("GeoStats")

#Execute FieldList
fieldList = arcpy.ListFields(Sediment1999, field_type="Double")

fieldNames = [f.name for f in fieldList ]

for name in fieldNames:

    # Set local variables - CreateGeostatisticalLayer
    inLayer = "C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\Kriging45.xml"
    inData = "C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\Sediment1999.shp F1="+name
    outLayer = name

# Execute CreateGeostatisticalLayer
arcpy.GACreateGeostatisticalLayer_ga(inLayer, inData, outLayer)

    # Set local variables - GALayerToGrid
    inLayer = name
    outGrid = (os.path.join("C:\\Users\\courtney\\Desktop\\MundleStableIsotope\\Anisotropic\\Sediment\\1999", name+".tif"))
    cellSize = 30
    cellptsHor = 1
    cellptsVer = 1

# Execute GALayerToGrid
arcpy.GALayerToGrid_ga(inLayer, outGrid, cellSize, cellptsHor, cellptsVer)