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
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-geostatistica...
just to make sure that you have the right inputs.
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
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.
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
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.
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)