topic Re: Universal Kriging with External Drift in ArcGIS GeoStatistical Analyst Questions

Universal Kriging with External Drift

Tae-Jung_JonathanKwon — Mon, 24 Dec 2012 01:59:01 GMT

Hi

I am looking for a way to use UK with external drift and the external drift. Or in others words, is there any way to perform a regression kriging on ArcGIS? (i.e., kriging with underlying trend..).
Thanks for your input in advance.

Re: Universal Kriging with External Drift

EricKrause — Wed, 26 Dec 2012 16:17:51 GMT

Universal kriging is available in the Geostatistical Wizard. It's one of the six kriging types. The trend is calculated from polynomials of the (x,y) coordinates; it does not support covariates other than the spatial location. However, covariates can be used as cokriging variables.

Re: Universal Kriging with External Drift

Tae-Jung_JonathanKwon — Wed, 26 Dec 2012 17:13:06 GMT

Thank you very much for your answer.
Could you please be little more specific as to what I should do?
You are right on that I do not need to use the trend derived from polynomials of the x and y coordinates.
This is being taken into account when I run the multiple linear regression where I consider 1. lat, 2 Long. 3. Dist to water. 4 Altitude, and 5. relative topography (differences in elevations between the target cell and nearby cells (within 1km). Here the dependent variable is mean air temperature. So I calibrated the model and and by using map algebra, I used the coefficients of the significant parameters (in my case, they were lat, alt, distance to water) to calculate the result of my regression for the whole area (thus producing a continuous surface).

Referring back to what you said, how could I use the continuous map generated from the regression as cokringing variables? It will be greatly appreciated if you could further explain or elaborate..

One more thing, is there any tool that allows me to calculate the correlation coefficient between the target cell and the nearby cells? For example, from MLR or Kriging, I can generate the interpolated map in which each cell contains a unique value (in my case it would be mean air temperature). I want to know if there is any tool available to see how much the value is each cell is related with averaged values in other neighbouring cells.

Thank you very much for your answers in advance and happy new year!

Universal kriging is available in the Geostatistical Wizard. It's one of the six kriging types. The trend is calculated from polynomials of the (x,y) coordinates; it does not support covariates other than the spatial location. However, covariates can be used as cokriging variables.

Re: Universal Kriging with External Drift

EricKrause — Wed, 26 Dec 2012 18:22:47 GMT

If you have a good MLR model already, I wouldn't try to use the covariates as cokriging variables. If you want to try it anyway, in the Geostatistical Wizard, when you choose kriging on the first page, you can enter up to four datasets. The first one you enter is the variable you will interpolate, and the three additional datasets will be used as cokriging variables.

"Kriging" is often called "residual kriging," and there's a reason for this: you always perform kriging on the residuals of some model. This model can be almost anything, but "regression kriging," "kriging with external drift," "universal kriging," and "linear mixed model kriging" all generally refer to the simultaneous estimation of the covariate coefficients and the kriging parameters. However, you may find success with doing a sequential estimation: first calculate the coefficients using your MLR model. Then calculate the residuals, and perform Simple kriging on these residuals (you should use Simple kriging instead of Ordinary kriging because you know that the mean of the residuals is 0). Then add these interpolated residuals back into the MLR predictions. You'll lose some power because you are sequentially calculating parameters (rather than simultaneously estimating them), but you should still get defensible results.

As for comparing the values in one point to the average of the neighboring points, the Semivariogram/Covariance cloud is probably the best way to visualize this, but the result is a graph rather than a single correlation coefficient. If you really need to calculate the correlation coefficient (and you're ok with ignoring spatial correlation in the analysis), we have a tool called Neighborhood Selection that selects the neighbors of an input (x,y) location (use the same neighborhood parameters that you used in kriging). It will probably take a lot of work, but I'm sure you can write a Python script that will do what you're trying to do. I've never personally done this, so I don't want to try to outline an algorithm, but I'm sure all the tools are there to accomplish this task.

Re: Universal Kriging with External Drift

EricKrause — Wed, 26 Dec 2012 18:25:26 GMT

Also, if you're going to do kriging on the residuals of your MLR model, recalculate the model without using Lat as a covariate. Otherwise you'll be "double-counting" (for lack of a better phrase) the spatial location.

Re: Universal Kriging with External Drift

EricKrause — Wed, 26 Dec 2012 18:33:57 GMT

Also, we heavily suggest using a projected coordinate system rather than a Lat-Long GCS. Distance calculations get badly distorted when using Lat-Long, and this distortion will get propagated through all the kriging calculations.

Re: Universal Kriging with External Drift

Tae-Jung_JonathanKwon — Fri, 28 Dec 2012 14:22:51 GMT

[ATTACH=CONFIG]20260[/ATTACH]
Dear Eric6346,

Thanks for your kind and detailed explanations:D

I would just like to confirm with you the method that you suggested; a sequential estimation by performing simple-kriging residuals
and combining this with MLR predictions. Here I have listed a few technical issues that I am facing;
1. I have created a grid surface using fishnet with 1km by1km covering my entire study area (500 cells by 500 cells)
2. I have 100 weather stations from which I calculated the mean air temperature (using 5 years of temperature data) for each station
3. I then performed MLR in SPSS; dependent variable being mean air temperature and independent variables being lat, altitude (using DEM layer), distance to water (using "near" function; each station to the nearest large water body).
4. From MLR, I now have an equation and with this, I calculate the residuals. The residuals were calculated by taking the difference between the observed and estimated (using the MLR equation). So I have 100 residuals.
5. With the MLR equation, I also generate an interpolated gridded surface; For each and all cells (1km^2) created in step1, and I extracted lat, alt, and distance to water thereby producing mean surface temperature estimated map on entire 500 cells by 500 cells. (please see the attached sample showing in a bigger grid size)
6. Now with the 100 residuals, I perform Simple Kriging thereby generating the continuous residual estimated map.
7. Then I assign the continuous surface found in step 6 onto each and every cell (i.e., 500 cells by 500 cells). This can be done by taking the average of the estimated residuals that fall into each cell.
8. Now simply add the products calculated in Steps 5 and 7.

Do these steps seem make sense to you?

Later on, I intend to perform a sensitivity analysis by using different grid size (5km by 5km, 10km by 10km, and 20km by 20km).
For this reason, it would be better if I could generate a MLR continuous surface (same resolution as continuous surface generated by Kriging) so that I can simply add both maps and assign them to whichever gird size that I want to test on. But I do not know the way to calculate the continuous MLR map (probably it is because I performed MLR outside of ArcGIS but in other software like SPSS).

As for the second question I had about calculating the correlation coefficient. Let me explain what I intend to do with this.
With the map generated in above steps (i.e., estimated mean air temperature map), I intend to analyse the current weather station setting and to recommend other potential sites. Let's assume that I want to install a weather station in "cold spots" meaning that a cell with the coldest estimated air temperature gets priority for weather station installation (this is why I intend to generate the interpolated air temperature map in the first place). In other words, I can have a rank for all cells by assigning the coldest cell with 1 and the warmest cell with say 100 (thus the rank ranges 1 to 100). But from the interpolated map, it is possible that a group of highly ranked cells are nearby each other. In this case, I then have to make a decision as to which cell gets more weight/priority; a cell with a high correlation coefficient gets priority as it would provide more "benefit" (since it is more representative) than a cell with a low correlation coefficient.
I will go ahead and explore the feasibility of what you suggested. But please suggest if there is any other way in ArcGIS to accomplish this task..

I am sorry for this long question, but your assistance would be greatly appreciated. Thanks again.

If you have a good MLR model already, I wouldn't try to use the covariates as cokriging variables. If you want to try it anyway, in the Geostatistical Wizard, when you choose kriging on the first page, you can enter up to four datasets. The first one you enter is the variable you will interpolate, and the three additional datasets will be used as cokriging variables.

"Kriging" is often called "residual kriging," and there's a reason for this: you always perform kriging on the residuals of some model. This model can be almost anything, but "regression kriging," "kriging with external drift," "universal kriging," and "linear mixed model kriging" all generally refer to the simultaneous estimation of the covariate coefficients and the kriging parameters. However, you may find success with doing a sequential estimation: first calculate the coefficients using your MLR model. Then calculate the residuals, and perform Simple kriging on these residuals (you should use Simple kriging instead of Ordinary kriging because you know that the mean of the residuals is 0). Then add these interpolated residuals back into the MLR predictions. You'll lose some power because you are sequentially calculating parameters (rather than simultaneously estimating them), but you should still get defensible results.

As for comparing the values in one point to the average of the neighboring points, the Semivariogram/Covariance cloud is probably the best way to visualize this, but the result is a graph rather than a single correlation coefficient. If you really need to calculate the correlation coefficient (and you're ok with ignoring spatial correlation in the analysis), we have a tool called Neighborhood Selection that selects the neighbors of an input (x,y) location (use the same neighborhood parameters that you used in kriging). It will probably take a lot of work, but I'm sure you can write a Python script that will do what you're trying to do. I've never personally done this, so I don't want to try to outline an algorithm, but I'm sure all the tools are there to accomplish this task.

Re: Universal Kriging with External Drift

EricKrause — Fri, 28 Dec 2012 16:01:54 GMT

The steps you've outlined are correct.

Unfortunately, there's no way to make a continuous OLS surface. You'll need to make separate rasters for each resolution you want to test, but you only need to krig on the residuals once. You can export the interpolated residual surface to any cell size and extent that you want.

I now understand what you're trying to do with the correlation coefficient. However, I don't think a correlation coefficient will work here because if you want to correlate a single point to the mean of its neighbors, you'll only be able to calculate a single coefficient for the entire surface (since you need repeated samples), so it won't help you in deciding which particular locations should be given preference.

The first thing that comes to mind is the Voronoi Map tool. It's an interactive graphical tool, and if you use Standard Deviation, Entropy, or Interquartile Range, you'll get an estimate of the local variability. A small local variability indicates that the predictions are more constant in that area, so they might be good candidates for new sites because the area can be better represented by a single value. Note that you'll need to convert your rasters to points to run the tool.

Re: Universal Kriging with External Drift

Tae-Jung_JonathanKwon — Sun, 30 Dec 2012 00:13:37 GMT

Dear Eric6346,

Thanks for all your explanations.

I at first intended to calculate correlation coefficient by including its direct neighbours (just like how Entropy is calculated) for each cell and for entire map.
Yes, your suggestion of using the Voronoi Map does really do what I wanted to do! 😄

So I generated the Voronoi Map which also contained Standard Deviation, Entropy, and etc. I looked up the descriptions for each stat generated from the Voronoi Map, and like you well mentioned, standard deviation and/or entropy would best describe the local variability. Since both measures explain the variability of the target cell and its neighbours, I would expect both measures be strongly correlated. However their degree of correlation wasn't so strong. (please see the attached diagram)[ATTACH=CONFIG]20319[/ATTACH]. I do understand that they are calculated in a different manner but shouldn't both measures have a strong positive correlation?

This leads me to wonder how both measures are calculated and which measure I should use for deciding the next potential weather station site in my study. I attached a sample Excel sheet in which I listed values for the target cell (i.e., center) and its neighbours. I also included the calculated values for Standard deviation and Entropy for the target cell. Would you be able to explain what was done in ArcGIS to calculate these two values? I looked up the description for Entropy, but it seemed like ArcGIS uses "Smart Quantiles" to categorize the values into five different classes. Unless I know the exact boundary of these five classes, it is difficult to apply the equation listed in this page..

Again, thank you very much for your continued support and helpful comments!

Jonathan

The steps you've outlined are correct.

Unfortunately, there's no way to make a continuous OLS surface. You'll need to make separate rasters for each resolution you want to test, but you only need to krig on the residuals once. You can export the interpolated residual surface to any cell size and extent that you want.

I now understand what you're trying to do with the correlation coefficient. However, I don't think a correlation coefficient will work here because if you want to correlate a single point to the mean of its neighbors, you'll only be able to calculate a single coefficient for the entire surface (since you need repeated samples), so it won't help you in deciding which particular locations should be given preference.

The first thing that comes to mind is the Voronoi Map tool. It's an interactive graphical tool, and if you use Standard Deviation, Entropy, or Interquartile Range, you'll get an estimate of the local variability. A small local variability indicates that the predictions are more constant in that area, so they might be good candidates for new sites because the area can be better represented by a single value. Note that you'll need to convert your rasters to points to run the tool.

Re: Universal Kriging with External Drift

Tae-Jung_JonathanKwon — Mon, 31 Dec 2012 03:20:36 GMT

I figured out how the software would calculate the statistics (at least in my case).. Taking an example of numbers in the attached Excel file (or below);

1       2        3
4       0        5
6       7        8

And the values assigned to each number are
0:      1.610996 (target cell)
1: 2.082779
2: 1.677525
3: 1.553863
4: 1.876214
5: 1.437138
6: 1.81845
7: 1.548435
8: 1.425129

If I were to calculate mean, it should be sum(0:8)/9= 1.670059. But the tool (Voronoi) ALWAYS ignores the values in 3 and 6 (top right corner and bottom left corner) when calculating mean thereby giving 1.665459 instead of 1.670059... This affects calculations for all other stats as well (e.g., stdev, entropy, and etc).. Any idea as to why this is happening and/or suggestion to fix this problem?
Many thanks!

Jonathan

Re: Universal Kriging with External Drift

EricKrause — Wed, 02 Jan 2013 14:35:37 GMT

I'll look into why the top-right and lower-left cells are being ignored, but after asking around, I think the tool you want to use is Focal Statistics. It gives lots of options for defining cell neighbors, and you can calculate the standard deviation of these neighbors.

And I'm not surprised at the fairly weak R^2 between standard deviation and entropy. Because entropy works with classified values (rather than raw values), the entropy map tends to be smoother. Entropy also has a maximum, but standard deviation has no maximum. So, two cells can have the same entropy value but still have very different standard deviations. If you look at your scatterplot, you can even see this; there are clear vertical columns that all have the same entropy, but the standard deviations vary a lot. This variance is what is pulling down the R^2.

Re: Universal Kriging with External Drift

EricKrause — Thu, 03 Jan 2013 17:02:09 GMT

I've looked into why the Voronoi map is ignoring the top-right and lower-left cells. The problem is that when you have gridded points, the Delaunay triangulation is not unique. Since we define polygon neighbors at the triangulation step (the first step of creating the Voronoi polygons), our algorithm drops the top-right and lower-left neighbors. A different but analogous implementation would drop the upper-left and lower-right polygons. We could fix this by defining neighbors after the polygons are created, but this would slow down the tool. We'll have to think whether this hit in performance is worth it, especially considering that Focal Statistics is specifically built to deal with gridded data.

Re: Universal Kriging with External Drift

Tae-Jung_JonathanKwon — Fri, 04 Jan 2013 16:51:04 GMT

Thank you very much for your kind explanations. Yes I understood why those cells are dropped.
Yes use of Focal Statistics would do the job since it gives standard deviation (which really indicates the local variability with its neighbours. I will go ahead and try it out!

Again, thanks a lot for all your support! Your support truly helped very much!!!
I will come back with more question if I have any 😄

Have a good weekend.

Jonathan

I've looked into why the Voronoi map is ignoring the top-right and lower-left cells. The problem is that when you have gridded points, the Delaunay triangulation is not unique. Since we define polygon neighbors at the triangulation step (the first step of creating the Voronoi polygons), our algorithm drops the top-right and lower-left neighbors. A different but analogous implementation would drop the upper-left and lower-right polygons. We could fix this by defining neighbors after the polygons are created, but this would slow down the tool. We'll have to think whether this hit in performance is worth it, especially considering that Focal Statistics is specifically built to deal with gridded data.

Re: Universal Kriging with External Drift

EricKrause — Mon, 07 Jan 2013 20:26:51 GMT

Glad I could help. Feel free to ask more questions if anything else comes up.

Re: Universal Kriging with External Drift

Tae-Jung_JonathanKwon — Wed, 09 Jan 2013 18:24:26 GMT

[ATTACH=CONFIG]20608[/ATTACH]
Hi Eric,
One quick question; following the steps mentioned above, I have krigged the residuals using 40 points (appeared as black pin). But as you can see in the attached image, the krigged region does not cover the entire study area. I know that it only allows to krig to the extent where the known points are.. Is there any way to extend this and to cover the entire study area? So far I have not found any option to do so.. Thanks!

Jonathan

The steps you've outlined are correct.

Unfortunately, there's no way to make a continuous OLS surface. You'll need to make separate rasters for each resolution you want to test, but you only need to krig on the residuals once. You can export the interpolated residual surface to any cell size and extent that you want.

I now understand what you're trying to do with the correlation coefficient. However, I don't think a correlation coefficient will work here because if you want to correlate a single point to the mean of its neighbors, you'll only be able to calculate a single coefficient for the entire surface (since you need repeated samples), so it won't help you in deciding which particular locations should be given preference.

The first thing that comes to mind is the Voronoi Map tool. It's an interactive graphical tool, and if you use Standard Deviation, Entropy, or Interquartile Range, you'll get an estimate of the local variability. A small local variability indicates that the predictions are more constant in that area, so they might be good candidates for new sites because the area can be better represented by a single value. Note that you'll need to convert your rasters to points to run the tool.

Re: Universal Kriging with External Drift

EricKrause — Wed, 09 Jan 2013 18:28:18 GMT

Yes, you can change the extent of a geostatistical layer. Right-click the layer in ArcMap's Table of Contents and choose "Properties." Go to the Extent tab and specify the new extent.

Re: Universal Kriging with External Drift

Tae-Jung_JonathanKwon — Wed, 09 Jan 2013 18:35:19 GMT

Sorry I didn't know it was that simple.. I should have looked into more closely.. I thought such feature would be available when running the kriging..

Thank you!

Yes, you can change the extent of a geostatistical layer. Right-click the layer in ArcMap's Table of Contents and choose "Properties." Go to the Extent tab and specify the new extent.

Re: Universal Kriging with External Drift

SteveLynch — Wed, 09 Jan 2013 21:05:00 GMT

Such features are available through the geoprocessing environment. I.e. when using GP tools. The Wizard, however, does not honor these environment settings.

You could, however, use the layer that you created in the GALayerToGrid tool and specify an output extent in the environment and your output raster will then have this new extent.

Steve

Re: Universal Kriging with External Drift

Tae-Jung_JonathanKwon — Fri, 11 Jan 2013 00:30:17 GMT

Steve,

Thanks for providing additional info. Yes it is definitely worthwhile using it since I have equal-sized grids (MLR surface), and I will assign a single value to each and every grid by taking the average of the converted raster data to get the final output of "regression-Kriging"...
Have a good weekend.

Jonathan

Such features are available through the geoprocessing environment. I.e. when using GP tools. The Wizard, however, does not honor these environment settings.

You could, however, use the layer that you created in the GALayerToGrid tool and specify an output extent in the environment and your output raster will then have this new extent.

Steve

Re: Universal Kriging with External Drift

MERLINDAVIS — Fri, 11 Jan 2013 15:08:16 GMT

Hello

I would like to know how exactly we can do the external drift kriging using ArcGIS. I couldnot figure it out(got confused) from the above discussion.
Also can anyone tell whether it is possible to do it for a time series of data. I have the hourly precipitation data for some staions. I would like to interpolate it for the entire area.

Thanks in advance! 🙂