Automating Empirical Bayesian Kriging

[DanaCarstens]

Hello, 

I am trying to write a Python script that will apply empirical bayesian kriging to >120 columns of water quality data. Each column is representative of a different month in either 1985 and 2015. I've gotten the script below to work successfully in ArcMap but it is producing temporary layers that go away once I close the program. How can I modify my script so that I can set a standard symbology for each layer produced and have it exported as a jpeg (or tif.) to some specified folder? Thank you!

import arcpy
from arcpy import env
import arcpy, os
from arcpy.sa import *

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

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

# Set workspace
env.workspace = r"K:/Research/2018 GIS"

# Set local variables
inPointFeatures = r"K:/Research/2018 GIS/DO_data.shp" # input shapefile
zFields = arcpy.ListFields("DO_data.shp",field_type="Double")

# EBK Veriables
outRaster = ""
inPointFeatures = "DO_data.shp"
cellSize = 8.60966128128307E-03
transformation = "NONE"
maxLocalPoints = 100
overlapFactor = 1
numberSemivariograms = 100

# Set variables for search neighborhood
searchNeighbourhood = "NBRTYPE=SmoothCircular RADIUS=0.474897239843659 SMOOTH_FACTOR=0.2"
outputType = "PREDICTION"
quantileValue = 0.5
thresholdType = "EXCEED"
probabilityThreshold = ""
semivariogram = "POWER"

for zField in zFields[3:]: #loop on the fields while ignoring the first three fields
outLayer = str(zField) + ".shp"
OutEBK = arcpy.EmpiricalBayesianKriging_ga(inPointFeatures, zField.name, outLayer, outRaster, cellSize, transformation, maxLocalPoints,overlapFactor, numberSemivariograms, searchNeighbourhood, outputType, quantileValue, thresholdType, probabilityThreshold, semivariogram)

[EricKrause]

Hi Dana,

Your "outLayer" is a geostatistical layer.  These are in_memory layers, and you will need to persist them somehow.

If you ultimately want rasters, you should probably just bypass the geostatistical layer entirely.  Instead of passing an empty string for "outRaster", instead just supply the file path, name, and file format where you want the raster.  Then pass the "outLayer" as an empty string.  

This will create all of the physical rasters on disk.  You can then import them as layers and apply any symbology that you want.

-Eric

Edit: If you really do want the geostatistical layers, you should add a step to save the geostatsitical layer to a layer file with the "Save to Layer File" geoprocessing tool.

[SteveLynch]

outLayer is a geostatistical layer (in memory). To persist it you'll use SaveToLayerFile, or in the script convert it to contours using GALayerToContour or specify an output raster in the call to EBK.

-Steve

[DanaCarstens]

Thank you both for your reply. I tried to use SaveToLayerFile and ran into some issues. Do I need to incorporate it into my script like this?

for zField in zFields[3:]: #loop on the fields while ignoring the first three fields
outLayer = str(zField) + ".shp"
out_layer_file = str(zField) + ".lyr"
OutEBK = arcpy.EmpiricalBayesianKriging_ga(inPointFeatures, zField.name, outLayer, outRaster, cellSize, transformation, maxLocalPoints,overlapFactor, numberSemivariograms, searchNeighbourhood, outputType, quantileValue, thresholdType, probabilityThreshold, semivariogram)
arcpy.SaveToLayerFile_management(outLayer, out_layer_file, "ABSOLUTE")

I would like the new layer to have the same name as zField. 

[SteveLynch]

Dana

I'd change outLayer = str(zField) + ".shp" to outLayer = str(zField) because it is an in-memory geostatistical layer.

Secondly, outRaster  is empty which means that you opted only to output a ga layer and therefore you don't need to specify a cell size as it only pertain to the output raster.

-Steve

[DanaCarstens]

Hi Steve,

I tried removing the .shp and I'm receiving the same error  " ExecuteError: ERROR 999999: Error executing function. Failed to execute (SaveToLayerFile). "

I'm not very familiar with Python. Do I need to run the SaveToLayerFile tool through another loop so that it is applied to all of my EBK output layers? 

for zField in zFields[3:]: #loop on the fields while ignoring the first three fields
outLayer = str(zField)
out_layer_file = str(zField)
OutEBK = arcpy.EmpiricalBayesianKriging_ga(inPointFeatures, zField.name, outLayer, outRaster, cellSize, transformation, maxLocalPoints,overlapFactor, numberSemivariograms, searchNeighbourhood, outputType, quantileValue, thresholdType, probabilityThreshold, semivariogram)
arcpy.SaveToLayerFile_management(outLayer, out_layer_file, "ABSOLUTE")

I appreciate your help. 

Dana 

[SteveLynch]

Dana

outLayer = str(zField.name)

out_layer_file = str(zField.name) + ".lyr" if you are using ArcGIS Desktop (...or ".lyrx" if using ArcGIS Pro)

BTW, ListFields returns a field object so you should rather be using zField.name like you're doing in the call to EBK

-Steve