Steve Lynch: Is it possible to do Cokriging script?

[YumeiHu]

Hi Steven Llynch,

I need to interpolate more than 1000 precipitation maps with Cokriging method, but I don't know if it is possible to use python scripit programming, because I can't find this tool in ArcToolbox.
I saw your old reply about this topic as following:
       "Message - Create a cokriging layer using the Wizard
        - Use the Get and Set Model Parameter tool to modify parameters if required
        - use Create Geostatistical Layer tool to modify the input dataset or variable used in the cokriging layer "
However, I don't really understand, could you give me an example in python script?
I will be very grateful for this.

Thank you very much!

Yumei Hu

[BryceStath]

I'm trying to automate a bunch of GA processes and my initial script appears to be working.  But how can I view or add my output GA layer to ArcMap?  I know this sounds silly, but I cannot figure out how to view my output, tmpLyr after the process runs.  Either it's not producing it (which if I run the script a second time, Python tells me tmpLyr already exists) or I just don't know how to add it to ArcMap.
I don't really need to see it so much as I will need to export it to a raster in the next step. 

thx

(edited: This code works and I can view the output grid just fine but it might still be nice to be able to view the GA layer mid way though the process)

# Purpose: Create a new Geostatistical layer based on a previous one

# Create the Geoprocessor object.
import arcgisscripting
gp = arcgisscripting.create()

gp.workspace = "C:\\working\\Y2010\\ECSR\\Plant"

#gp.toolbox = "ga"

try:
    # Set the input GA layer.
    inputGALayer = "C:\\working\\Y2010\\ECSR\\xml\\DS\\DS01_Pre2010.xml"

    # Set the new input dataset
    inputDset = "C:\\working\\Y2010\\ECSR\\Plant\\ECSR_Plant2009_working.gdb\\DS_P120_2009\\DS01_P120_2009_Plant"

    # Set the field name
    inputDset = inputDset + " BV_4ft"

    # Set output layer name
    outLayer = "C:\\working\\Y2010\\ECSR\\rasters\\temp\\tempLyr"
              

    # Check out Geostatistical Analyst extension license.
    gp.CheckOutExtension("GeoStats")

    # Process: Create a Geostatistical layer...
    gp.GACreateGeostatisticalLayer( inputGALayer, inputDset, outLayer)

except:
    # If an error occurred while running a tool, then print the messages.
    print gp.GetMessages()


try:
    # Set the input GA layer.
    inputGALyr = "C:\\working\\Y2010\\ECSR\\rasters\\temp\\tempLyr"

    # Set the output grid name.
    outputGrid = "C:\\working\\Y2010\\ECSR\\rasters\\temp\\out_grid"

    # Set some of the parameters.
    cell_size = 1
    points_horiz = 1
    points_vert = 1

    # Check out GeoStatistical Analyst extension license.
    #gp.CheckOutExtension("GeoStats")

    # Process: GA Layer To Grid...
    gp.GALayerToGrid_ga (inputGALyr, outputGrid, cell_size, points_horiz, points_vert)

except:
    # If an error occurred while running a tool, then print the messages.
    print gp.GetMessages()

[SteveLynch]

The geostatistical layer is an in memory layer. To persist it, use the SaveToLayerFile tool, viz.,

Usage: SaveToLayerFile <in_layer> <out_layer>


-Steve

[BrennaForester]

I need to automate the production of my xml files (I have 120 xml files to produce, which will then be applied to cokrig more files using the code discussed above).

The only parameter I am changing in the cokriging process is to press the "optimize model" button.

Can I cokrig two shapefiles & "optimize" the model using Python, then alter the script to loop through my 120 files?

Thank you,
Brenna

[EricKrause]

I need to automate the production of my xml files (I have 120 xml files to produce, which will then be applied to cokrig more files using the code discussed above).

The only parameter I am changing in the cokriging process is to press the "optimize model" button.

Can I cokrig two shapefiles & "optimize" the model using Python, then alter the script to loop through my 120 files?

Thank you,
Brenna

You only need to create a single xml.  Run the wizard with cokriging, and save the xml at the end.  Open the xml with a text editor, and the second line of code should read:

<model name="Kriging">

Add the following optimize tag to that line (which is case-sensitive):

<model name="Kriging" optimize = "ByCrossvalidation">

Save the xml, and use it as the model source in the Create Geostatistical Layer tool.  Then provide two new cokriging datasets, and the output will be a geostatistical layer that will be the same as if you had pressed the optimize button on the new data.

Just keep looping the CGL tool and providing it the original altered xml, and all the outputs will behave as if you had pressed the optimize button in the Wizard with the new datasets.

More information can be found in these links:
http://forums.arcgis.com/threads/11981-Spatial-Analyst-Interpolation-versus-Geostatistical-Analyst-W....
http://blogs.esri.com/Dev/blogs/geoprocessing/archive/2010/06/18/Automating-geostatistical-interpola...

NOTE: The optimize button was added in ArcGIS 10, so this process won't work for any versions before 10.

[BrennaForester]

Eric - thank you!  You have saved me a huge amount of work.

One last question - I would like to export/save prediction errors for each model: mean, root-mean-square, mean standardization, RMS standardization and average standard error.

What script would I add to output these errors; perhaps save as a txt file.

Brenna

[SteveLynch]

Brenna

Have a look at the Cross validation tool that is new at 10
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Cross_Validation/00300000000z000000

in Python you'll do something like
import arcpy
cvResults = arcpy.CrossValidation_ga("CoKriging")
print cvResults.rootMeanSquare
print cvResults.the other ones you need

and so on (you will write these to a file)

Steve

[BrennaForester]

Thank you for the feedback, Steve & Eric - my script is working perfectly!
-Brenna

[RobertPazur]

Dear all,
i followed all your recommendation, but according to my result I think that change of auto="true" in XML is not working properly. The think is that if I do so, all my semivariograms parameters differ in comparison to the parameters which is possible to obtain through Geostatistical Wizard. I perform this "test" on modelling of an indicator ( of course i leave the auto="false" value for the indicator threshold) and my parameters differ (f.e. lag size in geostatistical wizard is 129.14 and in python script 1225.65)   Could somebody have some idea why? If somebody interests i also attaching example files. To run the script  just simply unzip the dictionary to c:\.
Kind regards,
Robert Pazur

[liutiejun]

Hi Steve,

    How can I save the outlayer into a folder where i specify in the script? I have set the outlayer
    Viz.
    outLayer = "D:/uk/out/m1"
    But after the execution of my script without any error messages, I find nothing in the folder "d:/uk/out/".

Bests,
Tiejun Liu
####################

    inputGALayer = "D:/uk/Universal Cokriging.xml"

    ids1 = "D:/2000_1995/samples1029/culd/1/pro/1.shp"
    ids2 = "D:/uk/culd2005.shp"

# Set the field name
    inputDset1 = ids1 + " culdArea"

    inputDset2 = ids2 + " GRID_CODE"

# Set output layer name
    outLayer = "D:/uk/out/m1"

# Check out Geostatistical Analyst extension license.
    gp.CheckOutExtension("GeoStats")

# Process: Create a Geostatistical layer...
    gp.GACreateGeostatisticalLayer(inputGALayer,inputDset1 + ";" + inputDset2, outLayer)
##############################

Hi Steve Lynch,

With your help, I have finally find the way to do  the programming, I think I can finish it soon.
Thank you very very....much!

Best regards,
Yumei

[SteveLynch]

Run the CreateGeostatisticalLayer geoprocessing tool and then look at the syntax reported in the processing dialog.

Steve