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

[AnushuyaRamakrishnan]

Here is the XML input that I used.

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]

Thank you for your suggestion Eric

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

Thanks

Anushuya

[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

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

[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

[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

[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

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

[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