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

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

Vind meer berichten getagd met

Accepted answers

All comments

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

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

N_G_

I found the problem, and how to answer to my question above.
I put my code below if that can help ...

To use a list in the GACreateGeostatisticalLayer:

  #Set environment
    arcpy.env.workspace = "C:/INTERP/PJ.gdb"
    # Check out Geostatistical Analyst extension license.
    gp.CheckOutExtension("GeoStats")
    
    #Variables and loop
    inputGA = "C:/INTERP/COK.xml"
    
    # fc = list component/ fcs = list
    for fc in fcs:
    
        #1st variable /"PJ" = targetfield
        inputDset1 = fc + " "+ "X=Shape Y=Shape "+ "F1=" +"PJ"
        
        #2nd variable
        inputDset2 = "C:/INTERP/data.gdb/dem"
        VARCOVAR = inputDset1 + ";" + inputDset2
        
        #outLayer
        out = "cok"  
        outLayer = out + fc
        
        # Process : GACreateGeostatisticalLayer
        for fc in fcs:      
            gp.GACreateGeostatisticalLayer(inputGA,VARCOVAR,outLayer)
        print gp.GetMessages()

note : in this case it's not useful to make field list and value list for the target field

N_G_

Dear Steve and all,

I have read with attention posts about cokriging. I also have to produce cokriging with more than 2500 variables (Elevation+rainfall).

I'm not sure to be in the good thematic but I try, because some of you have probably solved (I think) my problem in their script.

I'm new in python and I don't manage lists. I would like to use them to create geostatical layer. I used model builder (models work) but I must write now cokriging script without model builder (export from model builder to python script is not enough...).

So, I made list of shapefiles and list of values for the field rainfall. But I have this answer in the python shell:

Traceback (most recent call last):
  File "C:\INTERP\sCRIPT PYTHON\FC4.py", line 75, in <module>
    inputDset1 = fc + NewList
TypeError: coercing to Unicode: need string or buffer, list found

I think that my code is not good and my lists must be transformed before to use them to create geost. layers. But I don't know how and where.......

If someone have an idea, a solution? It would be great.

My code:

import sys, string, os, arcgisscripting
# Create the Geoprocessor object
gp = arcgisscripting.create()

def listFcsInGDB():
    # list all Feature Classes in a geodatabase 
    gp.workspace = "C:\\INTERP\\PJ_class.gdb"

    fcs = []
    for fds in gp.ListDatasets('','feature') + ['']:
        for fc in gp.ListFeatureClasses('','',fds):
            #yield os.path.join(fds, fc)
            fcs.append(os.path.join(fds, fc))
    return fcs

fcs = listFcsInGDB()
for fc in fcs:
    print fc

# List of value in fields
input_dataset =  fc   
Atts = 'PJ_RacPJ'    #field with rainfall values       

rows = gp.searchcursor(fc)
row = rows.next()

NewList = []
for fc in fcs:
    for row in gp.SearchCursor(fc):
        fcValue = row.getvalue(Atts)
        NewList.append(fcValue)
#print NewList

arcpy.CheckOutExtension("GeoStats")

# Load required toolboxes
arcpy.ImportToolbox("C:/Documents and Settings/ArcGIS/Toolbox.tbx")

# Local variables:
Cokriging_xml = "C:\\INTERP\\COK_.xml"
CK = "CK2"
ids1 = "C:\\INTERP\\data.gdb\\dem1000"


for fc in fcs:
    for fcValue in NewList:
        inputDset1 = fc + NewList #Variables
        inputDset2 = ids1            #Covariable
        InputDset = "inputDset1;inputDset2"

# Create cokriging layer (code from model builder)

    tempEnvironment0 = gp.autoCommit
    arcpy.env.autoCommit = "1000"
    tempEnvironment1 = gp.spatialGrid1
    arcpy.env.spatialGrid1 = "0"
    arcpy.GACreateGeostatisticalLayer_ga(Cokriging_xml,InputDset,CK)
    arcpy.env.autoCommit = tempEnvironment0
    arcpy.env.spatialGrid1 = tempEnvironment1

Thanks

Naho

SteveLynch

Please use the help and search for terms like delete, remove.

liutiejun

Dear Steve,

Could you please tell me how to release or remove the outlayer from memory.

Bests,
Tiejun Liu

Yes it is in memory. To persist it use the SaveToLayerFile geoprocessing tool

Steve

liutiejun

Dear Steve,

Thank you very much.

Bests,
Tiejun Liu

Yes it is in memory. To persist it use the SaveToLayerFile geoprocessing tool

Steve

SteveLynch

Yes it is in memory. To persist it use the SaveToLayerFile geoprocessing tool

Steve

liutiejun

Hi Steve,

I have found the mistake in my script, I have misspelled the folder name of input layer.
However, another question puzzles me that where the outlayer is. Does it only exist in memory as an object?

Bests,
Tiejun Liu

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

Steve

SteveLynch

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

Steve

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

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

semivariogram_issue.zip

BrennaForester

Thank you for the feedback, Steve & Eric - my script is working perfectly!
-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

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

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-Wizard.
http://blogs.esri.com/Dev/blogs/geoprocessing/archive/2010/06/18/Automating-geostatistical-interpolation-using-template-layers.aspx

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

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

SteveLynch

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

Usage: SaveToLayerFile <in_layer> <out_layer>

-Steve

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

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

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

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

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

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

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.