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
