Export Kriging semivariogram averaged values as table using Python

[LuísPaixão]

Hello all.

I'm a student currently in a project to analyze agricultural parcels in order to obtain relations in biomass, vegetation indexes, water content, production, etc...

I´m fairly new to python scripting, only started using python 6 months ago. Nevertheless i´m trying to use arcpy to automatically convert a series of rasters into points and then perform an ordinary kriging "analysis" to each dataset.

I said analysis because i don't need to do the actual kriging interpolation. All i want is to able to save into a table the model parameters (range, sill, nugget) as well as the averaged semivariogram values. I can get the model parameters, using a xml template created in the Geoprocessing Wizard, but can´t seem to find any way to save the averaged values as table. I would be thankful for any input that can help me solve my problem.

Don´t know if it is important but i´m using ArcGIS 10.1. Also, here is my code so far:

import os
import arcpy

arcpy.CheckOutExtension("GeoStats")
path = r"D:\AgroInsider\Palminha_Mapas\IMG\NDVI\1275"
out=r"C:\Downloads\Geostatistics\point"
model=r"C:\Downloads\Geostatistics\NDVI_model.xml"

rasters=[]
shps=[]

walk = arcpy.da.Walk(path, topdown=True, datatype="RasterDataset")
for dirpath, dirnames, filenames in walk:
   for filename in filenames:
      rasters.append(os.path.join(dirpath, filename))

for rast in rasters:
   nome=os.path.splitext(os.path.basename(rast))[0]
   split=nome.split('_')
   split1=split[1]+"_"+split[2]+"_"+split[3]+"_"+split[4]+"_"+split[5]
   rast=arcpy.RasterToPoint_conversion(rast+"\\Band_1", out+"\\"+split1+".shp")

   arcpy.MakeFeatureLayer_management(rast,nome)
   arcpy.SelectLayerByAttribute_management(nome, "NEW_SELECTION", ' "GRID_CODE" = 1 ')
   arcpy.DeleteFeatures_management(nome)

walk = arcpy.da.Walk(out, topdown=True, datatype="FeatureClass")
for dirpath, dirnames, filenames in walk:
   for filename in filenames:
      shps.append(os.path.join(dirpath, filename))

for shp in shps: ####   Geostatistic segment
   nome=os.path.splitext(os.path.basename(shp))[0]
   indata=shp+" GRID_CODE"
   arcpy.GACreateGeostatisticalLayer_ga(model, indata, "nome"+nome)
   xmlPath = "/model[@name = 'Kriging']/model[@name = 'Variogram']/value[@name = 'Nugget']"
   nugget = arcpy.GAGetModelParameter_ga("nome"+nome, xmlPath)
   print "Nugget ", nugget
   xmlPath = "/model[@name = 'Kriging']/model[@name = 'Variogram']/value[@name = 'NumberOfLags']"
   lags = arcpy.GAGetModelParameter_ga("nome"+nome, xmlPath)

   xmlPath = "/model[@name = 'Kriging']/model[@name = 'Variogram']/value[@name = 'LagSize']"
   lagSize = arcpy.GAGetModelParameter_ga("nome"+nome, xmlPath)

   xmlPath = "/model[@name = 'Kriging']/model[@name = 'Variogram']/model[@name = 'VariogramModel']/value[@name = 'Range']"
   rnge = arcpy.GAGetModelParameter_ga("nome"+nome, xmlPath)
   print "Range ", rnge
   xmlPath = "/model[@name = 'Kriging']/model[@name = 'Variogram']/model[@name = 'VariogramModel']/value[@name = 'Sill']"
   sill = arcpy.GAGetModelParameter_ga("nome"+nome, xmlPath)
   print "Sill ", sill

Best Regards

[EricKrause]

Sorry to say this after you've done so much coding already, but the binned and averaged semivariances shown in the Geostatistical Wizard can only be exported through the user interface of the wizard.  They can't be exported with Python. 

If doing this manually is a possibility, you can add a step to save the geostatistical layer as a layer file using the Save To Layer File geoprocessing tool.  You can then add the layers to ArcMap and right-click -> Method Properties.  This will open the Geostatistical Wizard for that layer, and you can export the semivariogram table.

[LuísPaixão]

Thank you for the quick reply Eric. That is a possibility, but not very viable. I have hundreds of rasters to analyze and every week i get a dozen more. This is a two year project, so you see my problem.

Let me ask...

Is there any way i can do the maths myself adding a bit more coding to the script?

Or do you know another language/program that allows me to achieve my goal?

Thanks again.

Best regards

[EricKrause]

I don't want to give specific recommendations of what software and packages to use, but you should do some research into "empirical semivariograms."  The algorithm to compute them is not terribly complicated, and I'm sure you'll be able to find it implemented from a reliable source.

[LuísPaixão]

Thanks Eric, i'll see what i can find.

Best regards