Spatial Analyst/Interpolation versus Geostatistical Analyst Wizard.

[CedricWannaz]

Dear all,

     I have interpolations to perform of several hundred gridded (global 2x2.5 degrees) datasets that I need as 4000x2000 pixels rasters. I am currently testing a script that loops over all the datasets. I started with the Geostatistical Analyst's Wizard (Ordinary Kriging, default/automatic parameters) on one particular dataset and I obtained a smooth interpolation as shown on the first page of the attached document. The parameters (range, sill, ..) automatically determined by the wizard seemed to work. I then used the Spatial Analyst/Interpolation/Kriging tool, hoping that  parameters (that seem to be optional) would be determined the same way, and I obtained what is displayed on the second page of the attached document. There is clearly a structure appearing. I remember having tuned interpolation parameters for getting rid of such structures years ago (and under MATLAB), but I had to detect them visually before acting, and I cannot do this for all the datasets that have to be interpolated.

So my question: as my script is based on the Spatial Analyst tool (I don't think that the wizard's internal tool(s) is/are available as methods of the Geoprocessor (?)), is there a method(ology) for me to evaluate interpolation parameters the way the wizard does (so I can setup gp.Kriging_sa() the same way)?

The issue for me with any "interpolation artifact" is that I have to perform then a zonal statistics with zones that can be as small as 10km by 10km "squares". I wanted to have at least one pixel center per zone, and for these small zones I can clearly not count on any averaging effect that would smoothen to some extend the artifacts.


Thank you and best regards,

Cedric

[EricKrause]

Can you elaborate on what you mean by "structure"?  Your kriging surfaces look quite smooth, so I'm not sure what artifact you're talking about.

Also, the semivariogram appears to have three distinct strata... any idea where that might be coming from?

You might also be interested in this recent blog post:
http://blogs.esri.com/Dev/blogs/geoprocessing/archive/2010/06/18/Automating-geostatistical-interpola...

If you want to use the same kriging parameters (range, sill, lag, etc) for all of your interpolations, that post will help.

[CedricWannaz]

Thank you for your answer. I attach to this post a zoom that shows the structure that I mentioned. It appears when I use the Spatial Analyst/Interpolation/Kriging tool with no parameter specified except the cell size,  but not when I use the Geostatistical Analyst Wizard and then export to a raster.

Datasets are components along different directions of a vector field of wind speeds. Each dataset represents a layer (there are 55 of them at different altitudes), and I am performing an interpolation of each layer to have horizontal wind speed components at any point of them. I am not a specialist of atmospheric dynamic and I cannot explain the starta .. a cheap guess would be that we have three "modes/types" of flow (at the scale of the grid we can not really talk about laminar/turbulent I guess), each one having its own amplitude of variation of the speed component that I plotted.

Finally, I cannot use a template, because interpolation parameters will be different for each layer and each direction. In any case, as I am programming the loop, it wouldn't be an issue to leave constant parameters for each interpolation.

Best regards,

Cedric

[EricKrause]

At the Search Neighborhood screen inside the Wizard, try changing "Neighborhood type" to Smooth.  If that doesn't work (or if your version does not have the Smooth option), try increasing the Major and Minor semiaxes.  If that still doesn't fix it, let me know.

For reference, which version of ArcGIS are you using?

[CedricWannaz]

Actually, the Wizard is working perfectly well (its output, after exporting to raster, is smooth with default parameters). Problems come when I use the tool Spatial Analyst/Interpolation/Kriging instead, because it seems that it does not compute/provide default parameters the same way the wizard does.

As I can only use the aforementioned tool in my script but not the wizard, my question is: is there a way to find the set of default parameters that is generated by the wizard, so I can use them to setup the tool for it to work as well as the wizard?

My version of the software is ArcGIS Desktop 9.3.1 (built 3000), with an ArcInfo license.

Best regards,

Cedric

[EricKrause]

Here is my recommendation, and I wish there was a simpler way to do this. 

Use the Wizard to create a kriging layer.  On the Method Properties screen in the Wizard, save the xml file on your harddrive.  Then open the xml with a text editor, and wherever you see auto="false", change it to auto="true", and save the xml file.  The "auto" flag tells the software whether to update with new default values.

Run the Create Geostatistical Layer tool using the altered xml as the model source, and point to a new set of data.  This will create a kriging layer with the default parameters for the new data, and you can convert this layer to raster.

You should be able to automate this within Python, just keep using the original altered xml file in your loop.

[CedricWannaz]

Thank you for your answer. It seems that there is a detail (or more) that I don't get; I develop more what I am doing below, with illustrations in the attached file, hoping that it will help to determine what I don't understand.

I have two options to perform an Ordinary Kriging interpolation:
  A. The Geostatistical wizard: figures A.1 to A.9.
  B. A tool from the toolbox Spatial Analyst: figures B.1 to B.2

When performing A.1 to A.9, I obtain a smooth output (A.9).
When performing B.1 to B.2, I obtain an output that shows some structure (B.2).

My conclusion, at first, was that option A (the wizard) generates more appropriate default interpolation parameters (sill, range, lags, ..) than option B (the tool from the Spatial Analyst toolbox).

As it seems that I can only use the tool from option B in my script (whose prototype is shown in B.3), my question was: is there a way to have the tool from B computing default parameters the way A (the wizard) does. Or is there a method of the Geoprocessor that would provide me with parameters similar to what the wizard generates (but it has to be the output of a method (a function), because it is not an operation that I can repeat by hand), so I can setup the tool form B with them.

Reading your answers, I thought at first that B (the tool from Spatial Analyst) could be using an XML file outputed by the wizard (with auto="true") and would then be computing default parameters the way A (the wizard) does. But I can see no such input in B (on screenshot B.1 or in B.3 for the method), and I am unable to find in A.1 to A.9 where it is possible for me to output an XML file with interpolation parameters.

I hope that this explanation shows you clearly what I don't understand.

There is an additional point that comes to my mind: is it possible that the apparent smoothness in A comes from the conversion to raster that is done in a way "as smooth as possible" while still being within the estimation of the error (?), whereas in B we just have the best estimate? With Kriging I can expect the interpolant to have a symmetrical pattern; the surface is moving towards the mean when getting far form data points, and the symmetry comes from the regularity of the grid. The reason I used earlier the word "artifacts" is that I was expecting smoother "ups and downs" I guess and I was surprised to see a more geometrical pattern. If this is correct, then I have no more question because it means that B (the one that I can use in my script) is a better interpolant for me than A (that is no more a best estimate).


Thank you and best regards,

Cedric

[EricKrause]

The default parameters for the SA tool and the GA wizard are not computed in the same way, and I don't know of any way to export the parameters from the wizard to the SA tool.  You also need to understand that Kriging is our bread-and-butter here in Geostatistical Analyst; the wizard offers many more parameters and options than the Spatial Analyst tool.  So even if you could get the wizard parameters, many of them could not be used in the SA tool.

Again, I think the best (and simplest) thing to do here is to use the Create Geostatistical Layer tool in conjunction with an altered xml.  When you click Finish in the wizard, a "Method Summary" window will pop up.  On that window, click "Save", and you will be able to save the xml on your harddrive.  Go in with a text editor, and change all "auto" flags to "true".  Feed this altered xml into the Create Geostatistical Layer tool as the model source, and give it a new set of data.  The "auto=true" flag will tell the CGL tool to update default parameters for the new dataset.  Then convert this geostatistical layer to raster.  When looping through your datasets in python, just keep running the CGL tool with the original altered xml and converting the output to raster. 

As for why the SA tool is giving this geometrical pattern, I'm not sure, but I doubt it has much to do with the interpolation parameters like range, nugget, and sill.  My intuition is that it has to do with the search radius and the spatial orientation of the dataset.

[CedricWannaz]

I understand now; part of my confusion came from the fact that I was trying to apply your answers using partly the Spatial Analyst toolbox and I should have sticked to Geostatistical Analyst.

Thank you very much for your answers! This is now perfectly working as a script and I have no problem using the trick of the XML file (I am printing this post as a PDF for the record, and will update my script later if your toolbox allows once to have auto set to true through methods' parameters).

Best regards,

Cedric