AnsweredAssumed Answered

Export Kriging semivariogram averaged values as table using Python

Question asked by luispaixao on Dec 14, 2017
Latest reply on Dec 15, 2017 by luispaixao

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

Outcomes