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

[SteveLynch]

Create a Geostatistical Layer using cokriging. This is done once off using the Geostatistical Wizard.

Then in a Py script please do the following.
- use the Set Model Parameter to change any of the kriging parameters in the above layer (if required)
- use Create Geostatistical Layer tool to use the above and the data you want to interpolate

Repeat the above steps for each dataset.

PLease look at the examples in the Help.

Steve Lynch

[GerritLouw]

Hi  Steve,

I am trying to automate the cokriging process to test the different kriging parameters. I have created a geostatistical layer with the wizard and now I am trying to change the parameters with the Set Model Parameters Tool, however I am having trouble with the Xpath for some of the parameters.

Here is my code for setting the parameters:

params = "/model[@name = 'Kriging']/model[@name = 'Variogram']/value[@name = 'Anisotropy'];\

/model[@name = 'Kriging']/model[@name = 'Variogram']/model[@name = 'VariogramModel']/value[@name = 'ModelType'];\

/model[@name = 'Kriging']/model[@name = 'Variogram']/value[@name = 'LagSize'];\

/model[@name = 'Kriging']/model[@name = 'Variogram']/value[@name = 'NumberOfLags']"

The Xpath for the last two (LagSize and NumberOfLags) work fine, however the first two (Anisotropy and ModelType) does'nt work.

Can you show me where I went wrong?

Additionally, I am trying to print the RMSE results of the cross validation tool with some python script I found on the ESRI resource center, but it doesn't work:

results = arcpy.CrossValidation_ga(GACreateGeostatisticalLayer2)

print "RMSE: " + str(results.rootMeanSquare)

Any ideas?

Best regards,

Gerrit
(Stellenbosch)

[SteveLynch]

Gerrit

Please email me the XML file, I.e after you hit Finish in the Wizard then on the MethodReport dialog choose Save and send me (slynch@esri.com) this XML file.

For your second question,

GACreateGeostatisticalLayer2 is the name of the Geostatistical layer and not a variable I guess and therefore it should be in quotes, viz.

results = arcpy.CrossValidation_ga("GACreateGeostatisticalLayer2")

-Steve

[GerritLouw]

Hi Steve,

I've send you the XML.

With regard to my second question, I now get this error:

AttributeError: ResultObject: Get attribute rootMeanSquare not found

However, if I use the option to write the results out to the point file it does work.

Best regards,

Gerrit

[GerritLouw]

Thanks for the help!

Here is the final XML...

params = "/model/model[@name='Variogram']/model[@name='VariogramModel']/bool[@name='Anisotropy'];\

/model/model[@name='Variogram']/model[@name='VariogramModel']/enum[@name='ModelType'];\

/model[@name = 'Kriging']/model[@name = 'Variogram']/value[@name = 'LagSize'];\

/model[@name = 'Kriging']/model[@name = 'Variogram']/value[@name = 'NumberOfLags']"

And the assigned values could be 'true;Stable;355.1;10'

I could'nt get the cross validation tool to work, it kept saying that it can't find the attribute rootMeanSquare, but I worked around it by calculating the RMSE from the point output.

Best regards,

Gerrit

[YumeiHu]

Hi,

Thank you for your replying. I know how to do the programming for kriging.However,

Do you mean the step -'Create a Geostatistical Layer using cokriging. This is done once off using the Geostatistical Wizard.' can not be done with python script programming? Because I have large amount of data, it is reallly hard to click one by one.

Thank you very much!

Yumei

[YumeiHu]

I get your meaning now, but I also get new problem, when I try to use new shapefile to create Geostatistical layer, it always tell me
"Failed to execute. Parameters are not valid.
ERROR 045001: Input dataset(s) error. Table of inputs is not complete."
I don't know why.
When I create the layer in Geostatistical wizard, I used precipitation and elevation two fields, I tried to replace them with the new shapefile, but it doesn't work!

Thank you!
Yumei Hu

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

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

try:
    # Set the input GA layer.
    inputGALayer = "C:/Users/Public/Yumei/yumei/mytest0308/GA.lyr"

    # Set the new input dataset
    inputDset = "C:/Users/Public/Yumei/yumei/RE-compute/Interpp3ptV2/south/Interp_southV2_03.shp"

    # Set the field name
    inputDset1 = inputDset + " AveP"

    inputDset2 = inputDset + " Alt2k"

    # Set output layer name
    outLayer = "Interp_s03"

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

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

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

[SteveLynch]

Create a cokriging Geostatistical layer, then use this in the Create Geostatistical Layer GP tool and you'll notice from the results posted to results dialog that you are missing a semicolon between the 2 input datasets.

Viz.,
GACreateGeostatisticalLayer "Ordinary Cokriging" "C:\\Data\ca_ozone_pts.shp f1;C:\Data\ca_ozone_pts.shp f2" GAL1



Steve

[YumeiHu]

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

Create a cokriging Geostatistical layer, then use this in the Create Geostatistical Layer GP tool and you'll notice from the results posted to results dialog that you are missing a semicolon between the 2 input datasets.

Viz.,
GACreateGeostatisticalLayer "Ordinary Cokriging" "C:\\Data\ca_ozone_pts.shp f1;C:\Data\ca_ozone_pts.shp f2" GAL1



Steve