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

Vind meer berichten getagd met

Accepted answers

All comments

SteveLynch

1) to automate the kriging process, have a look at http://blogs.esri.com/esri/arcgis/2010/06/18/automating-geostatistical-interpolation-using-template-layers/
2) for steps 10-14 do the following,
- you have a point feature class where the homes are, use this in the Extract values to points
- run the AddXY tool on the output feature class created above
- select all the rows in the output table and export, or if the output was sent to a shapefile you can use Excel to open the .dbf file

Hope this helps
Steve

AnushuyaRamakrishnan

Hi Steve,

Thank you very much for your prompt response.

Thank you for guiding me with steps 10-14.

I will try them.

Also, I will try model builder or python scripting and get in touch with you if I have any difficulty.

Thank you

AnushuyaRamakrishnan

Hi Steve,

Going back to Steps 10-14, I tried to use the extract values to points tool and used the kriged prediction raster as input raster.

Extraction was successful.

I am trying to add xy layer as suggested by you by right clicking the extracted points shp file.

I am not finding the option for it.

Also I am trying to use the sample tool to integrate all the extracted points shp file to generate a single file. A table is generated.

However, I am not able to convert it into ascii or dbf. When i try to export it, I am getting a 0KB file that I am not able to access.

Can you suggest suitably?

Thank you

Anushuya

Eric6346

Rather than converting the geostatistical layer to a raster and extracting the points, you should just use the "GA Layer to Points" geoprocessing tool. It will bypass the creation of the raster, and merging all of the outputs will be easy.

Also, you might want to consider using Empirical Bayesian Kriging (EBK) instead of Ordinary Kriging. The results are generally more accurate, and EBK is implemented as a geoprocessing tool, so it is much simpler to automate.

AnushuyaRamakrishnan

Thank you Eric!!! I will try running EBK and will keep you updated.

Thanks

Anushuya

AnushuyaRamakrishnan

Hi Eric/Steve,

I have another question:

I use a 99 gridded (12 km x 12 km grid) Ozone PM input file to predict concentration at each mother's residence.

The "Extract values to Points" output table generates raster values at each mother's location.

Do I multiply or divide this value by 144 (12 x 12) to adjust for the gridded area.

Please let me know on this.

Thanks

Anushuya

VHolubec_0-1713883009885.png

Eric6346

You don't need to do any adjustment. Ozone is a continuous phenomenon, so its value changes everywhere. When you extract the points from a raster, all you're doing is asking for the ozone value at the closest cell center. You don't need to do any corrections for the size of the grid cell.

Again, you can eliminate any imprecision from using the closest cell center by using GA Layer to Points. It will calculate the ozone value at the exact point rather than looking to the closest cell center.

AnushuyaRamakrishnan

Thank you for the clarification Eric

Anushuya

SteveLynch

I am trying to add xy layer as suggested by you by right clicking the extracted points shp file.

I am not finding the option for it.

It's a geoprocessing tool, http://resources.arcgis.com/en/help/main/10.1/index.html#//001700000032000000

AnushuyaRamakrishnan

Hi Steve/Eric,

I have created a model for my project using a model builder and incorporated the tools for Emprical bayesian kriging. I could run the model successfully and get the output that I expected.

Now, I need to loop the model to change the Z field of my merge1.CSV input point file (containing PM values for each day from Jan 1 2004 to Dec 31 2006). to run all the tools continuously to generate outputs for the entire dataset.

How I can loop through the merge1.CSV input point file in my model builder?

Thanks

Anushuya

AnushuyaRamakrishnan

Hi Steve/Eric,

I attach my python script for Emprical Bayesian Kriging.

I need to inlcude a last step for looping through the entire merge1.csv data (Jan 1 2004 to Dec 31 2006).

How do I do this looping?

Pls. suggest

Thanks

Anushuya

Emprical Bayesian Kriging_09162013.py.zip

Eric6346

This topic explains what you need to know:

http://resources.arcgis.com/en/help/main/10.2/index.html#/Examples_of_using_iterators_in_ModelBuilder/00400000001n000000/ESRI_SECTION1_33EDB1C8829F4CF4BBD5DABAA3ED34DD/

AnushuyaRamakrishnan

Hi Eric/Steve,

I am trying to follow the steps below for iteration of EBK:

# Create a list variable to store the names of all the fields in the air quality data table
fieldList = arcpy.ListFields(File path to your table)

# Loop through each field in the air quality table
for zField in fieldList:

    # Set output raster name to match the field name with an added EBK suffix
    outRaster = zField + "_EBK"

    # Run the EBK interpolation for each field.
    arcpy.EmpiricalBayesianKriging_ga(...)

When I try to run my model in python, arcmap 10.2 crashes... I do not know the reason.

Have you experienced this problem?

Can you suggest me suitably.

AnushuyaRamakrishnan

I have written python scripts for Empirical bayesian kriging and I am using it to generate predicted values for my sample mothers.
However, there is a concern that EBK is not well established in the literature.

I need certain clarifications before I can proceed to apply EBK for my study:

1. What is the underlying model to estimate parameters in semivariogram in EBK?

2. What is the bayesian Prior /bayesian rule that is applied in EBK?

3. What is the basis of calling it more accurate?

Any clarifications for these questions would be highly appreciated.

Thank you

Anushuya

Eric6346

1. If no transformation is applied, the underlying model is an intrinsic random function, specifically a Power semivariogram model. If a transformation is applied, it is a simple kriging model with an exponential semivariogram.

2. The prior distribution is estimated from the data using restricted maximum likelihood (REML). This is the difference between "Bayesian" and "empirical Bayesian." In empirical Bayesian models, the prior is estimated from the data, and the posterior distribution is generated through simulations.

3. There are several reasons to think EBK will, in general, be more accurate. First, REML is known to be a better estimator of semivariogram parameters than weighted least-squares (which is what is used in other kriging methods). Second, EBK works on subsets, so it eases the stationarity assumption of kriging. Other kriging methods assume global stationarity. EBK, however, only assumes stationarity within subsets, so it can work effectively even in the presence of global nonstationarity, as long as the subsets are close to being stationary. Third, EBK does a better job of estimating kriging variances because it can account for error in estimating the underlying semivariogram (it does this with simulations); other kriging models assume the semivariogram is modeled absolutely correctly, which often results in standard errors that are too small.

Here is a reference for EBK (a longer paper with all the details is in the works):

Krivoruchko K. and Gribov A. (2014) Pragmatic Bayesian kriging for non-stationary and moderately non-Gaussian data. Submitted. In Mathematics of Planet Earth. Proceedings of the 15th Annual Conference of the International Association for Mathematical Geosciences. Eds: Pardo-Igúzquiza, E.; Guardiola-Albert, C.; Heredia, J.; Moreno-Merino, L.; Durán, J.J.; Vargas-Guzmán, J.A. Springer 2014, pp. 61-64.

(We're not sure why the date is 2014 because it's already been published)

AnushuyaRamakrishnan

Hi Eric,

Thank you very much for the clarification.

This would help me a lot.

Thanks

Anushuya

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

Eric6346

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

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

Thanks

Anushuya

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

Eric6346

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,

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

Eric6346

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 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

VHolubec_1-1713883021042.png

GeoJosh_0-1713883318678.png

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

Eric6346

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

I am using arcmap 10.2

Thanks

Anushuya

Eric6346

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

Thank you Eric!!!

I will try this and get back to you

Anushuya

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?
>>>

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!!!!

Doe mee aan het gesprek

Aangemelde leden kunnen berichten plaatsen, updates volgen en meer. Nieuw hier? Registreer een gratis account.

Snelle Links

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