CE 397 Statistics in Water Resources

Exercise 9

Geostatistical Analysis

by:

John Shaw, Yao You, Ruth Haberman and David Maidment

University of Texas at Austin

April 2009

Contents

Introduction. 1

Basics of Kriging. 2

Case Study. 4

Geostatistical Analysis. 6

Simple Kriging and Inversed Distance Weight Methods. 28

To be turned in. 30

Introduction

In this exercise we will explore Spatial Interpolation Methods

Figure 1 : The interpolated value at the unmeasured yellow point is a function of the neighboring red points (From ArcGIS Help Menu).

A very basic problem in spatial analysis is interpolating a spatially continuous variable from point samples. For example, many spatially explicit hydrologic/watershed models require continuous surfaces of temperature. Three commonly used interpolation methods to model spatially distribution from point data are Inverse Distance Weighting (IDW), spline and ordinary kriging.

The IDW is simple and intuitive deterministic interpolation method based on the principle that sample values closer to the prediction location have more influence on prediction value than sample values farther apart. Using higher power assigns more weight to closer points resulting in less smoother surface. On the other hand, lower power assigns low weight to closer points resulting in smoother surface. We optimize the power parameter using ArcGIS. The major disadvantage of IDW is the “bull's eye” effect (higher values near observed location) and edgy surface. Spline is deterministic interpolation method which fits a mathematical function through input data to create smooth surface. A spline can generate sufficiently accurate surfaces from only a few sampled points and retain small features (Anderson, 2008). Spline works best for gently varying surfaces like temperature. In ArcGIS Spline is Radial Basis Function.

Unlike IDW and spline, kriging is method based on spatial autocorrelation. It uses a semivariogram.

Basics of Kriging

Kriging was developed in the 1960s by the French mathematician Georges Matheron.  The motivating application was to estimate gold deposited in a rock from a few random core samples.   Kriging has since found its way into the earth sciences and other disciplines.  It is an improvement over inverse distance weighting because prediction estimates tend to be less bias and because predictions are accompanied by prediction standard errors (quantification of the uncertainty in the predicted value).

The basic tool of geostatistics and kriging is the semivariogram.  The semivariogram captures the spatial dependence between samples by plotting semivariance against separation distance (semivariance will be explained in the next paragraph).  The premise of any spatial interpolation is that close samples tend to be more similar than distant samples (this is also called spatial autocorrelation).  This property of spatial data is implicitly used in IDW.  In kriging, one must model the spatial autocorrelation using a semivariogram instead of assuming a direct, linear relationship with separation distance.

Semivariance equals one-half the squared difference between points separated by a distance d±∆d (assuming no direction preference).  As the distance between samples increase, we expect the semivariance to also increase (again, because near samples are more similar than distant samples).  This is true, however, only up to some given separation distance.  For this distance and up, points are unrelated.  Stated another way, if 50m is this critical separation distance, two points separated by 50m are likely to be just as similar (or dependent on one another) as samples separated by 100, 200, 300, or any distance greater than 50m.

Suppose we have the semivariogram shown in Figure 2.  What information does the plot provide?  Well, the semivariance between samples separated by no distance is about 1.5E-4.  This is called the nugget.  What it says is that if you measure the variable at locations very, very close to one another, the values measured might be quite different.  Why would this happen?  Suppose you had a gold nugget in the middle of an otherwise gold-free rock.  If you sample just on the edge of the nugget you get a high gold estimate.  If you sample just outside of the edge, you get no gold in your estimate.  The presence of a nugget in the semivariogram therefore tells you that, assuming no measurement error, the variable is not spatially continuous.

The semivariogram also tells us that points separated by 60,000 m are likely to have the same average difference as points separated by 100,000, 150,000, 200,000 m or any distance above 60,000m.  60,000 m is the range of the semivariogram and suggests the area of influence for any given point.  An unmeasured location can be predicted based on its neighboring samples closer than 60,000m.  A sample collected 61,000 m away from the sample will likely have no influence on the actual value at the unmeasured location.

When you look at the model of a semivariogram, you'll notice that at a certain distance, the model levels out. The distance where the model first flattens out is known as the range Sample locations separated by distances closer than the range are spatially autocorrelated, whereas locations farther apart than the range are not. The value that the semivariogram model attains at the range (the value on the y-axis) is called the sill. The partial sill is the sill minus the nugget

Figure 2 : The semivariogram is used to model the spatial relationships between samples separated by some distance, d

For kriging estimation, the semivarogram model (the yellow line in figure 2) is used to obtain estimates for the weighting parameters of Equation 1.  This process is done automatically by the geostatistical analyst once the user is satisfied with the semivariogram.  If you are interested in the derivation of the weighting parameters (or any of the other topics discussed here), Applied Geostatistics by Edward H. Isaaks and R. Mohan Srivastava is an excellent resource.  Or for the more mathematical folks, try Statistics for Spatial Data by Noel A.C. Cressie.

Goals of this Exercise

Now we know the basics of spatial interpolation. Let’s use our knowledge to develop a digital elevation model of the Wax Lake Outlet Area using the 1963 bathymetric data.

Computer Requirements

This exercise is to be performed in ESRI ArcMAP 9.x with the Geostatistical Analyst extension.  The data for this exercise is at: http://www.ce.utexas.edu/prof/maidment/StatWR2009/Ex9/Ex9.zip (11MB).

The data used in this exercise were taken directly from a 1963 bathymetry chart of the Wax Lake Outlet Area which is located in the delta of the Mississippi River in Louisiana, as shown in Figure 3.

Case Study

Figure 3 :   from Wellner et al. 2005

The Wax Lake Outlet is a man-made channel that connects the Atchafalaya River to Atchafalaya Bay in southern Louisiana. The Atchafalaya River is connected to the Mississippi river in northern Louisiana by the Old River Control Structure.  At the Old River Control Structure, approximately 30% of the Mississippi’s water and suspended sediment is routed into the Atchafalaya basin. As shown in Figure 4, this sediment is building sub-areal deltas at the mouths of the Atchafalaya River and the Wax Lake Outlet at rates of roughly 10km2 sub-areal land per year. Understanding the mechanisms of this delta growth is seen as a key to better remediation of wetland loss in coastal Louisiana.

2002

1974

Figure 4 : Wellner et al. (2005)

We have uncovered a bathymetry chart of Atchafalaya Bay from 1963. We can make an estimate of how much sediment has accumulated in the delta since then by subtracting a current digital elevation model from an elevation model created from this bathymetric chart. This exercise will turn this chart into a digital elevation model.

Figure 4 : from USACE

Geostatistical Analysis

Goals of this Exercise

Now we know the basics of spatial interpolation.   Let’s use our knowledge to develop a digital elevation model (DEM) of the Wax Lake Outlet Area using the 1963 data.

Computer Requirements

This exercise is to be performed in ESRI ArcMAP 9.3 with the Geostatistical Analyst extension.  This is available in the CE-Learning Resource Center in Room 3.301.  The data for this exercise is at: http://www.ce.utexas.edu/prof/maidment/StatWR2009/Ex9/Ex9.zip

The data used in this exercise were taken directly from the 1963 bathymetry chart of the Wax Lake Outlet Area.  Additional data came from NHDPlus Region 8

http://www.horizon-systems.com/nhdplus/HSC-wthMS.php

Getting Started

In this exercise we need to use spatial analyst and geospatial analyst toolbar.   Open ArcGIS.  We need to enable this extension in the first place. Open Tools à Extension.

Check the ‘Spatial analyst’ and ‘Geospatial analyst’ box. Now we can use all the functions in both toolbars.

The digitized water depths from the bathymetric chart are contained in the spreadsheet Depth.xls  and consist of latitude and longitude points associated with a depth of water at that point, measured in feet from the water surface to the bed.  Since these data were mapped in 1963, the latitudes and longitudes would have been referenced to the North American Datum of 1927 (NAD27), the standard for earth datums at that time.

These data were extracted from the Excel file and converted to ArcGIS as point events, to become the feature class DepthEvents.   This feature class was coupled with another for NHDWaterbody representation to provide spatial context for the point set.  The NHD dataset is presented in geographic coordinates with the North American Datum of 1983 (NAD83).  Because spatial interpolation has to be done in a Cartesian (x,y) coordinate system rather than a (latitude, longitude) coordinate system, all the data have been projected into the State Plane Coordinate system for South Louisiana, using the NAD83 datum.  The DepthEvents were first converted from NAD27 to NAD83 datum before this projection to State Plane Coordinates was done.  The resulting dataset is stored in the WaxLakeOutlet feature dataset in the Ex9 geoddatabase, as shown below:

Open ArcGIS, and add the DepthEvents and NHDWaterbody layers to the map display.  The following view will appear on your screen after some manipulation of the legends and feature class labels:

Notice the point labels, which identify the depth in feet.  Most of the points are between 2 and 7, with an area at the very bottom of the dataset with much more extreme depths (62 feet, 29 feet, 14 feet).  These are not an anomaly, but reflect a deep trench off-shore.

Exploratory Spatial Data Analysis

Exploratory Spatial Data Analysis (ESDA) is a process of understanding the properties of a spatial dataset in order to best model the data using geostatistics.  The word “Explore” should tell you that ESDA is more of an adventure than a strict – you must follow the path at all times – procedure.  In this exercise we show how the Geostatistical Analyst tools can be used to understand the population distribution of the attribute of interest and how to understand the large-scale patterns in the dataset through 3D visualization.  This is not a comprehensive list of ESDA procedures, but a good start.  Through this process, keep in mind that the better one understands the spatial characteristics of the data, the better kriging model one can build to interpolate the data, and, consequently, the better estimates one will produce.

Histogram

Open Geospatial analyst à Explore Data à Histogram.

Select the 40 bars and check the statistics option to view the statistics of the surface elevations. Keep the transformation as None. Select layer as DepthEvents and Attribute as Depth.  A similar window to below will appear on screen if you increase the number of bars to 40.

We will now analyze the histogram for the surface elevations at the Wax Lake Outlet.  The upper right corner shows the statistics of the surface elevations.  The histogram shows that our data is not perfectly normally distributed but the skewness is not large. We should remember that we did not fully follow this assumption when analyzing these data.  One of the crosschecks of normal distribution of data is that mean should be close to the median. In our case mean is 5.7302 feet and median is 6 feet, so we aren’t too far off. We can consider our data as normally distributed. User can transform the data into log or box-cox distribution if it is not normally distributed.

QQ (Quantile-Quantile) plot

Another way to understand the data’s distribution is by using the Normal QQ Plot tool.  In a QQ (Quantile-Quantile) plot we test whether data are normally distributed by plotting it against a dataset with a known normal distribution.  If the plot is linear along the line Y=X, then the data follow a normal distribution.

Open Geospatial analyst à Explore data à Normal QQ plot.

The QQ plot shows the linear relationship between log(Depths_Events) and the standard normal distribution. The data is not exactly linear on the upper because of points from the large trench.

So, lets log-transform the data to allow for that:

Ok, this looks better.

Trend analysis

The trend analysis tool provides a 3D plot of the samples and a regression on the attribute in the XZ and YZ planes.  The purpose of the tool is to visualize the data and to observe any large-scale trends that the modeler might want to remove prior to estimation.  It is best to keep the kriging model as simple as possible and to only remove a trend if it significantly improves prediction errors.

Open Geospatial analyst à Explore data àTrend analysis.

The screen will look similar to the figure below.  We can choose the different graph options and rotate the location to see the trend in the data.   It looks like these data have some significant trends in them (which we would expect because water flows downhill!).

The trend analysis tool provides a 3D plot of the samples and a regression on the attribute in the XZ and YZ planes.  The purpose of the tool is to visualize the data and to observe any large-scale trends that the modeler might want to remove prior to estimation.  It is best to keep the kriging model as simple as possible and to only remove a trend if it significantly improves prediction errors.

Perform Kriging interpolation

Open Geospatial analyst à Geospatial Wizard

Select the Kriging as method and DepthEvents as input data. Select attribute as Depth.

Click Next to proceed.

Lets use a Log transformation of the data, and a First Order trend removal to allow for the trends that our earlier exploratory data analysis suggested.

We will choose the most widely used ordinary kriging method and prediction map. As discussed in

Click Next to proceed.

You’ll see a map that shows the fitted linear trend in the data

Hit Next, and lets choose the default Spherical model for the interpolation

Other commonly used models are exponential and Gaussian. Lag size and number of lags are the important parameters to be selected. These two parameters are used to group the number of pairs of data. For our analysis we will use the default values. For given lag size and number of lags, ArcGIS automatically calculates the nugget, range and sill.  Semivariogram for particular directions can be investigated using show search direction tool.

Click Next to proceed.

Change Neighbours to include to 12 and include at least 10 stations. Since we have the total 315 data points, 12 and minimum 10 stations is an appropriate number for interpolation. Click on the third circle for sector type.

Click Next to proceed.

These are the error statistics for the interpolation. This is based on one-leave-out  at a time cross validation method. We can save the cross validation table to perform further statistical analysis.

Click Finish.

You get a report of the analysis.

This is the summary of interpolation parameters chosen. Click Ok.

Now you will see the following digital elevation map (DEM) of the Wax Lake Outlet.

We will save this as a raster dataset. Right click on the Ordinary Kriging layer. Select Data à Export to raster.

Select output raster as folder and file name as Depths_Events

Add this layer to the map.

You can tweak the symbology to show shallow water as light and deep water as a darker blue by right clicking on the raster layer, and go to Properties à Symbology, and choosing an appropriate color ramp.

This is what you will see.

You can actually see the deep area at the bottom of the DEM, and the shallow areas that have been caused by sediment deposition.

Lets repeat the process and get the Prediction Standard Error Map

You can click on the map at any point in the interpolation and use the Identify tool to get the predicted depth and its standard error of estimate: in this case 6.8 ft +/- 3.0 ft.

Simple Kriging and Inversed Distance Weight Methods

Now go back to Geospatial analyst à Geospatial Wizard. Still choose Kriging and click next. This time choose Simple Kriging as the interpolation method. The rest is the same as Ordinary Kriging procedures. You will get a result similar to the following (after overlay the original data points on top):

Notice that this map agrees perfectly with the data points where the Ordinary Kriging method doesn’t. This is because in Simple Kriging the values at those original data points are kept the same while in Ordinary Kriging every point on is map is reassigned with an interpolated values. Discuss in the homework which one is more reasonable here.

We go back to Geospatial analyst à Geospatial Wizard again, and this time let’s try Inverse Distance Weight method. This method does not assume any distribution of the data. Try with different parameters and create one you like the most (record the parameters you used for your final result). The following is one example of the output.

Obviously this method does not keep the original values at the sample sites. You can learn more about this method by reading the document provided by the software. Discuss in your homework whether this method is better than other two methods for this exercise.

To be turned in

A histogram, normal QQ plot, trend analysis, for the data.

A semivariogram, method summary, a DEM, and a Prediction Standard Error Map developed using Ordinary Kriging, Simple Kriging, and Inverse Distance Method.  Discuss the difference between these techniques.  Note that Simple Kriging assumes we have 100% confidence in data, which we do not have for the 1963 bathymetry.  Ordinary Kriging does not assume that we trust the data 100%, and does a better job of “smoothing out" the output.

Try the same process, but see what happens if you don’t log transform the data or remove its trends. Do you think the result is significantly different from what you got earlier?