nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Need Help Automation of Kriging using Model builder or python - Reg.

AnushuyaRamakrishnan

I am using arcmap 10.2.

My project deals with the prediction of PM and Ozone concentrations for 2004-2006 (1098 days) using known concentrations of them at
99 grids throughout the county.

I am trying to predict the unknown concentrations at mother's residences using ordinary kriging.

This is what I am doing:

1. Use Geostatistical Wizard
2. Select kriging cokriging
3. Input the Data layers and specifications
4. Select Finish
5. Read the model parameters and click OK
6. Generate Kriged predictions
7. Right click on Kriged data and convert to raster
8. Use spatial analyst tool Extraction
9. Select "Extract values to points"
10. Generate Extracted values of kriged predictions
11. Select conversion tool
12. Choose convert From Raster to Ascii
13. Generate the Ascii table.
14. Alternatively use Sample tool and input multiple rasters to generate an output table.

My question:

1. I need some guidance for automating this kriging process using model builder or python?

2. Is there any way to directly copy and paste the raster values to excel?

Can any one provide me a lead into this

Encuentra más publicaciones etiquetadas con

Accepted answers

All comments

XIANWANG

Hello, I am reading these information, which are really useful for me. According to your discussion. I already can do the kriging using python for a given shpfile. I don't know how to do the loop. What you use for the loop. As here we can not use iterators from model bulider. Thanks a lot!!!!

XIANWANG

Hi, I am reading these information, which is very useful for me. But I am using version 10. There is no arcpy.GeostatisticalDatasets in version 10, right? What I should do for these general workflow:

newDataset = arcpy.GeostatisticalDatasets("C:\\temp\\myXML.xml")
newDataset.dataset1 = "C:\\temp\\myFeatureClass1.shp"
newDataset.dataset1Field = "myField1"
arcpy.GACreateGeostatisticalLayer_ga("C:\\temp\\myXML.xml", newDataset, "outGALayer1")

Thanks very much!!!!

AnushuyaRamakrishnan

Hi Steve:

I have not solved the problem that I reported earlier yet.

I confirm that there is problem with the python code. When I try to run kriging for certain datasets that already have kriged predictions generated, script is throwing me an error. The error is as follows:

Traceback (most recent call last):

File "C:\Users\aramakrishnan\Documents\ArcGIS\scratch\Ordinary Kriging_July2014.py", line 41, in <module>

arcpy.gp.ExtractValuesToPoints_sa(Harris_Lat_Long_csv_Events_lyr, Raster_Path + "OZ8_" + fieldname, Ext_Table_Path + "OZ8_" + fieldname, "INTERPOLATE", "VALUE_ONLY")

File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 498, in <lambda>

return lambda *args: val(*gp_fixargs(args, True))

ExecuteError: ERROR 999999: Error executing function.

Failed to execute (ExtractValuesToPoints).

Failed to execute (ScriptJul21).

Failed at Wed Aug 27 19:35:22 2014 (Elapsed Time: 1 minutes 55 seconds)

I tried to run ordinary kriging by manual methods by inputing the lag, sil and nugget values. However, I am not able to even generate the raster to feed it into the script. I attach my complete python code and xml file.

Will wait to hear from you

Thank you

Anushuya

KrigingordJan18.xml.zip

Ordinary Kriging_July2014.py.zip

SteveLynch

What happens if you run the ExtractValuesToPoints geoprocessing tool from the toolbox and feed in your raster and point feature class? I.e. in your script you could print out the raster name and the name of the feature class.

Also, the GALayerToPoints tool exports a geostatistical layer to points. The tool can also be used to predict values at unmeasured locations or to validate predictions made at measured locations.

This route will save you creating a raster.

-Steve

AnushuyaRamakrishnan

Hi,

I need some help with the error I get while using spatial analyst tool in a python script (attached).

Due to hard disk space issues, I deleted most of the old data files in the default.gdb folder and tried to run the script, which was working fine before the data clean-up.

Below is the error message I get.

Please help me identify what is missing and how to fix the issue.

Thanks,

Anushuya

Executing: Script
Start Time: Sun Aug 10 15:47:46 2014
Running script Script...

Traceback (most recent call last):
File "C:\Users\aramakrishnan\Documents\ArcGIS\scratch\Ordinary Kriging_July2014.py", line 41, in <module>
arcpy.gp.ExtractValuesToPoints_sa(Harris_Lat_Long_csv_Events_lyr, Raster_Path + "OZ4" + fieldname, Ext_Table_Path + "OZ4" + fieldname, "INTERPOLATE", "VALUE_ONLY")
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 498, in <lambda>
return lambda *args: val(*gp_fixargs(args, True))
ExecuteError: ERROR 999999: Error executing function.
Failed to execute (ExtractValuesToPoints).

Failed to execute (Script).
Failed at Sun Aug 10 15:49:38 2014 (Elapsed Time: 1 minutes 51 seconds)

Ordinary Kriging_July2014.py.zip

KrigingordJan18.xml.zip

SteveLynch

You have a typo in the first line, viz.,
kriging should be Kriging

The best way to check for this is,
- keep a copy of the xml file that is good
- use the CreateGeostatisticalLayer gp tool in ArcMap to verify that the modified xml file is good
- do it one step at a time.

Steve

MicheleLazzarini

I enclose also the xml file, because I guess this error should be here.

<model name="kriging" optimize = "ByCrossvalidation"><dataset Label="Dataset" dataset-type="DVA"/><dataset Label="Dataset 2" dataset-type="DVA"/><dataset Label="Dataset 3" dataset-type="DVA" optional="true"/><dataset Label="Dataset 4" dataset-type="DVA" optional="true"/><dataset Label="Decluster's Clipping Dataset" dataset-type="Generic" sub-type="polygon" optional="true" visible="false"/><dataset Label="Decluster's Clipping Dataset 2" dataset-type="Generic" sub-type="polygon" optional="true" visible="false"/><dataset Label="Decluster's Clipping Dataset 3" dataset-type="Generic" sub-type="polygon" optional="true" visible="false"/><dataset Label="Decluster's Clipping Dataset 4" dataset-type="Generic" sub-type="polygon" optional="true" visible="false"/><enum name="KrigingMethodType">Ordinary</enum><enum name="KrigingResultType">Prediction</enum><items name="Datasets"><item name="Dataset"><enum name="TrendType">None</enum><model name="NeighbourSearch" options="CopyFromVariogram"><enum name="Type">Standard</enum><value name="NeighboursMax" auto="true">5</value><value name="NeighboursMin" auto="true">2</value><enum name="SectorType">Four45</enum><value name="MajorSemiaxis" auto="true">19.8304866304385</value><value name="MinorSemiaxis" auto="true">19.8304866304385</value><value name="Angle">0</value></model>
</item><item name="Dataset"><enum name="TrendType">None</enum><model name="NeighbourSearch" options="CopyFromVariogram"><enum name="Type">Standard</enum><value name="NeighboursMax" auto="true">5</value><value name="NeighboursMin" auto="true">2</value><enum name="SectorType">Four45</enum><value name="MajorSemiaxis" auto="true">19.8304866304385</value><value name="MinorSemiaxis" auto="true">19.8304866304385</value><value name="Angle">0</value></model>
</item></items><model name="Variogram"><value name="DatalayerCount">2</value><value name="NumberOfLags" auto="true">12</value><value name="LagSize" auto="true">1.6525405525365418</value><enum name="PairsType" auto="true"><v>Semivariogram</v><v>Semivariogram</v></enum><bool name="NuggetOn">true</bool><value name="Nugget" auto="true"><v>1.1229850369211205</v><v>552.8246726873753</v></value><value name="MeasurementError100"><v>100</v><v>100</v></value><bool name="ShiftON">true</bool><value name="Shift" auto="true"><v>0</v><v>0</v><v>0</v><v>0</v></value><bool name="VariogramModelAuto">false</bool><model name="VariogramModel"><enum name="ModelType">Stable</enum><value name="Parameter" auto="true">2</value><value name="Range" auto="true">19.8304866304385</value><bool name="Anisotropy">false</bool><value name="Sill" auto="true"><v>1.6358684289646366e-006</v><v>0.06351430744842085</v><v>0.06351430744842085</v><v>2466.0096002988107</v></value></model>
</model>
</model>

Dear all,
I have a similar problem, I am trying to automatically process some shapefiles in order to use cokriging.
This is the python code:

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("GeoStats")
arcpy.env.overwriteOutput = True

# Local variables:
#CoKriging = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\CoKriging.lyr"
CoKriging = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\Krigingmodel2_2.xml"
TR1intermediate_point = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\TR1intermediate_point.shp"
GACreateGeostatisticalLayer1 = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\GACreateGeostatisticalLayer5"
#GALayerToPoints2 = "C:\\Users\\nkondapalli\\Documents\\ArcGIS\\Default.gdb\\GALayerToPoints5"
GALayerToGri2 = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\GALayerToGri2"

GALayerToPoints2_PointToRast = "C:\\Users\\nkondapalli\\Documents\\ArcGIS\\Default.gdb\\GALayerToPoints5_PointToRast"
test2ASCII_TXT = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\"#test2ASCII.TXT"

numFiles=5
# Process: Create Geostatistical Layer
for i in range(2,3):#numFiles+1):
    print "C:\\Users\\nkondapalli\\Desktop\NewFolder\\RG1intermediate_point_%d.shp X=Shape Y=Shape F1=DATA;C:\\Users\\nkondapalli\\Desktop\NewFolder\\TR1intermediate_point_%d.shp X=Shape Y=Shape F1=DATA"%(i,i)

    arcpy.GACreateGeostatisticalLayer_ga(CoKriging, "C:\\Users\\nkondapalli\\Desktop\NewFolder\\RG1intermediate_point_%d.shp X=Shape Y=Shape F1=DATA;C:\\Users\\nkondapalli\\Desktop\NewFolder\\TR1intermediate_point_%d.shp X=Shape Y=Shape F1=DATA"%(i,i), GACreateGeostatisticalLayer1)
    # Process: GA Layer To Points
    #arcpy.GALayerToPoints_ga(GACreateGeostatisticalLayer1, TR1intermediate_point, "", GALayerToPoints2)
# Process: GA Layer To Grid
    arcpy.GALayerToGrid_ga(CoKriging, GALayerToGri2, "0.063", "1", "1")
# Process: Point to Raster
    #arcpy.PointToRaster_conversion(GALayerToPoints2, "OBJECTID", GALayerToPoints2_PointToRast, "MOST_FREQUENT", "NONE", "0.05")

# Process: Raster to ASCII

    print 'print to file : '+'file%d.txt'%i
    arcpy.RasterToASCII_conversion(GALayerToGri2, test2ASCII_TXT+"file%d.txt"%i)

and this is the error:

C:\Users\nkondapalli\Desktop\NewFolder\RG1intermediate_point_2.shp X=Shape Y=Shape F1=DATA;C:\Users\nkondapalli\Desktop\NewFolder\TR1intermediate_point_2.shp X=Shape Y=Shape F1=DATA

Traceback (most recent call last):
File "C:\Users\nkondapalli\Desktop\NewFolder\testErik3.py", line 36, in <module>
    arcpy.GALayerToGrid_ga(CoKriging, GALayerToGri2, "0.063", "1", "1")
File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\ga.py", line 1341, in GALayerToGrid
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input geostatistical layer: Dataset C:\Users\nkondapalli\Desktop\NewFolder\Krigingmodel2_1.xml does not exist or is not supported
WARNING 000725: Output surface raster: Dataset C:\Users\nkondapalli\Desktop\NewFolder\GALayerToGri2 already exists.
Failed to execute (GALayerToGrid).

It looks like the error is in the xml file. We actually modified it, inserting optimize = "ByCrossvalidation" and changing the options auto to "true".
Any clue?
Thanks!

Michele

MicheleLazzarini

Dear all,
I have a similar problem, I am trying to automatically process some shapefiles in order to use cokriging.
This is the python code:

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("GeoStats")
arcpy.env.overwriteOutput = True

# Local variables:
#CoKriging = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\CoKriging.lyr"
CoKriging = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\Krigingmodel2_2.xml"
TR1intermediate_point = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\TR1intermediate_point.shp"
GACreateGeostatisticalLayer1 = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\GACreateGeostatisticalLayer5"
#GALayerToPoints2 = "C:\\Users\\nkondapalli\\Documents\\ArcGIS\\Default.gdb\\GALayerToPoints5"
GALayerToGri2 = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\GALayerToGri2"

GALayerToPoints2_PointToRast = "C:\\Users\\nkondapalli\\Documents\\ArcGIS\\Default.gdb\\GALayerToPoints5_PointToRast"
test2ASCII_TXT = "C:\\Users\\nkondapalli\\Desktop\NewFolder\\"#test2ASCII.TXT"

numFiles=5
# Process: Create Geostatistical Layer
for i in range(2,3):#numFiles+1):
    print "C:\\Users\\nkondapalli\\Desktop\NewFolder\\RG1intermediate_point_%d.shp X=Shape Y=Shape F1=DATA;C:\\Users\\nkondapalli\\Desktop\NewFolder\\TR1intermediate_point_%d.shp X=Shape Y=Shape F1=DATA"%(i,i)

    arcpy.GACreateGeostatisticalLayer_ga(CoKriging, "C:\\Users\\nkondapalli\\Desktop\NewFolder\\RG1intermediate_point_%d.shp X=Shape Y=Shape F1=DATA;C:\\Users\\nkondapalli\\Desktop\NewFolder\\TR1intermediate_point_%d.shp X=Shape Y=Shape F1=DATA"%(i,i), GACreateGeostatisticalLayer1)
    # Process: GA Layer To Points
    #arcpy.GALayerToPoints_ga(GACreateGeostatisticalLayer1, TR1intermediate_point, "", GALayerToPoints2)
# Process: GA Layer To Grid
    arcpy.GALayerToGrid_ga(CoKriging, GALayerToGri2, "0.063", "1", "1")
# Process: Point to Raster
    #arcpy.PointToRaster_conversion(GALayerToPoints2, "OBJECTID", GALayerToPoints2_PointToRast, "MOST_FREQUENT", "NONE", "0.05")

# Process: Raster to ASCII

    print 'print to file : '+'file%d.txt'%i
    arcpy.RasterToASCII_conversion(GALayerToGri2, test2ASCII_TXT+"file%d.txt"%i)

and this is the error:

C:\Users\nkondapalli\Desktop\NewFolder\RG1intermediate_point_2.shp X=Shape Y=Shape F1=DATA;C:\Users\nkondapalli\Desktop\NewFolder\TR1intermediate_point_2.shp X=Shape Y=Shape F1=DATA

Traceback (most recent call last):
File "C:\Users\nkondapalli\Desktop\NewFolder\testErik3.py", line 36, in <module>
    arcpy.GALayerToGrid_ga(CoKriging, GALayerToGri2, "0.063", "1", "1")
File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\ga.py", line 1341, in GALayerToGrid
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input geostatistical layer: Dataset C:\Users\nkondapalli\Desktop\NewFolder\Krigingmodel2_1.xml does not exist or is not supported
WARNING 000725: Output surface raster: Dataset C:\Users\nkondapalli\Desktop\NewFolder\GALayerToGri2 already exists.
Failed to execute (GALayerToGrid).

It looks like the error is in the xml file. We actually modified it, inserting optimize = "ByCrossvalidation" and changing the options auto to "true".
Any clue?
Thanks!

Michele

AnushuyaRamakrishnan

Hi Eric,

I followed your suggestions and ran the loop for 7 days and checked few days manually.

I got exactly matching results.

Thank you for your help and suggestions.

I appreciate it

Thanks

Anushuya

AnushuyaRamakrishnan

Thank you for your suggestion Eric

I will try the options that you suggest and let you know if it works.

Thanks

Anushuya

EricKrause

Sorry for taking so long. I've been busy and had to debug this problem.

You can recreate the defaults from the Geostatistical Wizard by adding one line to your XML file.

You should see a line at the top of the XML that says:
<model name="Kriging">

Change this line to:
<model name="Kriging" optimize = "BySill">

You actually don't need to change any "auto" flags. This "optimize" flag will override any auto flags, so it doesn't matter if they are false or true.

If you want to automate the Optimize Model button in the Geostatistical Wizard, you can change that same line to:
<model name="Kriging" optimize = "ByCrossvalidation">

AnushuyaRamakrishnan

Here is the XML input that I used.

Thanks

Anushuya

KrigingOrdMay1.zip

EricKrause

Your Python script looks fine. Can you post your XML file?

AnushuyaRamakrishnan

Hi Eric,

I generated my ordinary kriging output successfully using the python script for generating geostatistical layer and using it to generate raster and convert it to points and generate the table.

However, I see that the output generated by scripting does not exactly match with the output that I generate manually by running
ordinary kriging for the same day using the same values of range, sill and nugget.

In the script I only set auto=true for lagsize, number of lags, range, sill and nugget values.

I attach the python script for your review.

Please suggest.

Thank you

Anushuya

Ordinary Kriging_NoLoop.py.zip

AnushuyaRamakrishnan

Hi,

Do you input a raster or a layer or a data file like CSV.

Please go to add tools and add your specific input.

Then go to add xy tool and project it in arcmap.

Then perform specific geostatistical operation using the events layer and get the output layer

Next right click on layer and convert it to raster.

Use this raster to generate your output table.

In any module of arcmap the screen should not go blank.

Try these steps and check if it works.

Also try and reinstall arcmap if the problem persists.

Thanks

Anushuya

shabijunejo

hi, I am new to arcgis, I was working with georeferencing tool,,, i have followed all the process but in the and when I apply project raster the whole screen turns black and when I zoom to full extent it shows only a tiny part of it.... please help me if you may...

thanks

AnushuyaRamakrishnan

Thank you Eric... I will do that

Anushuya

EricKrause

The creation of a new geostatistical layer is the automation: that new geostatistical layer is the result of kriging on the new dataset. Once you have the geostatistical layer, you can export it to a raster with GA Layer to Grid, or you can predict directly to a point feature class with GA Layer to Points.

Now you just need to iterate all your datasets through the Create Geostatistical Layer tool.

AnushuyaRamakrishnan

Hi Eric,

I could succesfully generated new geostatistical layer from an existing xml input.

Can you guide me to the next step of kriging.

For kriging using SA I use the following script:

# Process: Kriging
arcpy.gp.Kriging_sa(merge3_csv_Events_lyr, "longitude", fieldname, "Exponential 0.003512 0.255986 1.092512 0.001093", "3.5120478474024E-03", "VARIABLE 12", Output_variance_of_prediction_raster)

1. I need the script for kriging using ga

2. What is the input for this script? Is it the same as SA.

Pls. let me know.

Thank you

Anushuya

EricKrause

You should change your code to something like:

X = arcpy.GeostatisticalDatasets("C:\\temp\\KrigingordinaryJan152.xml") #Make sure your xml is saved in C:\temp for this to work
X.dataset1 = "C:\\temp\\PM04Jan16.shp" #Make sure your shp file is in C:\temp
X.dataset1Field = "Y04PMJan16"
arcpy.GACreateGeostatisticalLayer_ga("C:\\temp\\KrigingordinaryJan16.xml", X, "PM04Jan16")

You can replace the X with any name you want. The first line creates a GeostatisticalDatasets object called "X". It populates the necessary properties by reading your XML. The second line sets the dataset1 property of X. The third line sets the field for X. And the last line uses the XML source and the GeostatisticalDatasets object ("X"), and it creates a geostatistical layer called "PM04Jan16".

AnushuyaRamakrishnan

Hi Eric,

I tried to run this in python:

# Import arcpy module
import arcpy

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

merge1.csv.events = arcpy.GeostatisticalDatasets("C:\\temp\\KrigingordinaryJan152.xml")
Harris_Lat_Long.csv.events = "C:\\temp\\PM04Jan16.shp"
merge1.csv,events.Field = "Y04PMJan16"
arcpy.GACreateGeostatisticalLayer_ga("C:\\temp\\KrigingordinaryJan16.xml", newDataset, "PM04Jan16")

It is throwing me an error:

Traceback (most recent call last):
File "C:\Users\aramakrishnan\Documents\ArcGIS\scratch\Ordinary Kriging_09252013_01.py", line 19, in <module>
merge1.csv.events = arcpy.GeostatisticalDatasets("C:\\temp\\KrigingordinaryJan152.xml")
File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\_ga.py", line 300, in GeostatisticalDatasets
return GADatasets(gapy.GALayerInfo(ga_model_source))
TypeError: bad argument type for built-in operation

Am I missing something?
>>>

AnushuyaRamakrishnan

Thank you Eric!!!

I will try this and get back to you

Anushuya

EricKrause

Ok, that makes things easier. We created the GeostatisticalDatasets arcpy class for 10.2, and it will be useful for you.

Here is the general workflow (notice that you don't need to use any XPATH commands):
1. Manually create an XML model source and change all the auto flags to "true". (You already did this)
2. Here is a quick shell of the Python code (replace the file paths and names with your data):

newDataset = arcpy.GeostatisticalDatasets("C:\\temp\\myXML.xml")
newDataset.dataset1 = "C:\\temp\\myFeatureClass1.shp"
newDataset.dataset1Field = "myField1"
arcpy.GACreateGeostatisticalLayer_ga("C:\\temp\\myXML.xml", newDataset, "outGALayer1")

In the above code, "myXML.xml" is the XML file that you created in step 1. "newDataset" is the GeostatisticalDatasets object that controls the datasets and fields that you want to interpolate. The second line sets the dataset1 property to one of your feature classes, and the third line sets the field. The fourth line creates the new geostatistical layer by using myXML as the model source and the newDataset object as the data. It then outputs a geostatistical layer called "outGALayer1". Again, you only need to make the XML file once, and you don't need any XPATH code.

3. You need to iterate this so that it cycles through the different datasets and fields. This will create a geostatistical layer for each dataset where the interpolation parameters have been recalculated for each new dataset.

AnushuyaRamakrishnan

I am using arcmap 10.2

Thanks

Anushuya

EricKrause

I see that you're using the Kriging tool from Spatial Analyst. This tool is not the same as the kriging methods implemented in the Geostatistical Wizard. If you want to use that tool, it's simple to automate, but it doesn't have advanced options like detrending, transformations, anisotropy, etc. The workflow I'm describing is about how to automate kriging from the Geostatistical Wizard.

I'll talk more about the workflow, but I need to know what version of ArcGIS you have.

AnushuyaRamakrishnan

HI Eric,

I am still not clear about the flow of python script for ordinary kriging.

# Process: Kriging
arcpy.gp.Kriging_sa(merge1_csv_Events_lyr, "Y04Mar02", Kriging_PM04Mar2, "Spherical 0.003512 0.339420 4.530740 0.517570", "3.5120478474024E-03", "VARIABLE 12",
Output_variance_of_prediction_raster)

(merge_1_csv_Events_lyr contains PM data for all the days of the study; Y04Mar02 is the Z field for one day)

I attach the XML file that I generated for one day of observation as ZIP file.

Do we built a model in model builder and give that model name as input for xmlpath.

What is happening in these steps:

xmlPath = "/model[@name = 'Kriging']/model[@name = 'Variogram']/value[@name = 'Range']"
newValue = 0
outModel = "C:/gapyexamples/output/outModel.xml"
model source="C:\\Users\\aramakrishnan\\Desktop\\Grid combo1\\merge1.csv Events"
xmlPath= "/model[@name = 'Kriging']/model[@name = 'Variogram']/model[@name = 'VariogramModel']/value[@name = 'Sill']"
newValue = 1
outModel = "C:/gapyexamples/output/outModel.xml"
model source="C:\\Users\\aramakrishnan\\Desktop\\Grid combo1\\merge1.csv Events"
xmlPath= "/model[@name = 'Kriging']/model[@name = 'Variogram']/value[@name = 'Nugget']/@auto;\
newValue = 2
outModel = "C:/gapyexamples/output/outModel.xml"

For this step what is the inlayer? Do we carry out this step after calling in range,sill and nugget values.

Execute CreateGeostatisticalLayer
arcpy.GACreateGeostatisticalLayer_ga(inLayer, inData, outLayer)

Please help me to understand this better.

Thank you

Anushuya

KrigingOrdinaryJan152.zip

AnushuyaRamakrishnan

Hi Eric,

Thank you for the suggestion. I will certainly try that.

EBK is relatively new and has not been reported in scientific literature for prediction of pollutant concentrations in adverse birth outcomes studies that we are involved in.

When we send our manuscript for review, there will be questions on the validity of the model and the applicability of it in related studies.

As ordinary kriging as an established model reported in a lot of manuscripts, use of it in the study will eliminate question on the validity of it.

Thanks

Anushuya

EricKrause

You might be able to convince them with crossvalidation and validation. EBK typically gives better results, and this is how you can demonstrate it. Take a couple of your datasets and interpolate them with Ordinary Kriging and EBK and compare the results.

AnushuyaRamakrishnan

Hi Eric,

Thank you for the detailed clarification.

I absolutely agree with you. EBK is a sophisticated model inbuilt in the geostatistical wizard and is very easy to automate.

Hopefully my advisors will be in favor of EBK, we are yet to reach a decision on EBK vs ordinary kriging.

Thanks Eric

Anushuya

EricKrause

The Create Geostatistical Layer tool takes a model source as input. This model source can be a geostatistical layer (either a layer in ArcMap or saved as a .lyr file on disk) or an XML model source. For your purposes, the XML will be easier to use. The tool reads all the interpolation parameters from the model source (type of kriging, nugget, range, sill, transformations, etc) and applies them to a new dataset. This is useful for something like temperature data taken daily. You can build the model for one day and easily apply that model to each subsequent day.

However, if you're using different kinds of data (temperature, elevation, pollution, etc), you don't want to keep using the same model over and over because the interpolation parameters will not fit different types of data. So, you need to tell the tool to recalculate all these parameters (explained below) for each new dataset.

You only need to manually create one XML file source. You can do this with the Geostatistical Wizard. Open the Wizard, choose Ordinary kriging, and give it a dataset. When you click Finish, the Method Report screen will pop up that shows all the parameters that you used. Click "Save..." and save the XML file in a convenient location. You then need to open the XML file with a text editor (Cooktop is a useful and free XML editor, but you can do it with Notepad too). Inside the XML, you'll see all of the parameters, and most of them will have auto = "false" after them. This auto flag tells the tool whether to recalculate that parameter or keep it fixed when it is used as a model source. For every parameter that you want to be recalculated, you need to change the flag to auto = "true" and save the XML. You only need to do this once, and you'll keep reusing this same model source.

You can then set up a loop in Python to iterate through your datasets. For each dataset, you'll use the XML as the model source in Create Geostatistical Layer. This will generate a geostatistical layer for each dataset, where the interpolation parameters have been recalculated, and you can use these layers to create the ASCII files that you described in your first post.

The difficulty in automating Ordinary/Simple kriging is one of the many reasons we made Empirical Bayesian Kriging as a geoprocessing tool.

AnushuyaRamakrishnan

Hi Eric,

I have few questions:

What is the source file called kriging.layer mentioned in the example? Does it refer to the layer file containing all the days of observation of data? or does it refer to one day of the Z field?

How do we get the Z field (containing PM and ozone data for different days of observation) into the model source= ..... script.

Can we generate the XML files you mentioned for a bunch of days and loop them for kriging process, or are we required to change the input and generate the xml files for Range, sill and nugget one day at a time.

Just wanted a clarification on this.

Thank you

Anushuya

AnushuyaRamakrishnan

Thank you Eric... I will try this and let you know if I run into trouble.

Thanks

Anushuya

EricKrause

Automation of Ordinary Kriging can be done with a template model and the Create Geostatistical Layer geoprocessing tool. This blog will take you through the basics, but you'll need to do some extra work:

http://blogs.esri.com/esri/arcgis/2010/06/18/automating-geostatistical-interpolation-using-template-layers/

However, you'll need to use a template XML file rather than a geostatistical layer. You can save this XML file at the end of the Geostatistical Wizard (after you hit Finish, click "Save..." on the Method Report screen). You'll need to change various "auto" flags in this XML file. You can change them with the "Set Model Parameter" gp tool, or you can just open the XML file in a text editor and manually change them (this way is probably easier).

Read the following topic (be sure to read the whole thing because I make a mistake in my first post):
http://forums.arcgis.com/threads/75743-Automatically-generate-anisotropic-semivariogram-parameters

The other poster's question was about anisotropic parameters and multiple semivariogram models, so the solution is more complicated that what you need. See if you can figure it out from those two resources, but I'm happy to help you if you run into any problems.

AnushuyaRamakrishnan

How do we automate ordinary kriging? Just was curious to know about it.

We need to generate range, sill and nugget values for each day of the study. Is there any way for me to generate those as a first step in python script and give the inputs in the subsequent step for kriging.

I am trying to figure out how this can be done?

Any suggestions/help would be highly appreciated.

Thank you

Anushuya

AnushuyaRamakrishnan

Hi Eric,

Thank you very much for the clarification.

This would help me a lot.

Thanks

Anushuya

Únete a la Conversación

Los miembros registrados pueden publicar, seguir actualizaciones y más. ¿Nuevo aquí? Regístrate gratis.

Enlaces Rápidos

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.