CE 397 Statistics=
in
Water Resources
Exercise 9
Geostatistical An=
alysis
by:
John Shaw, Yao You, Ruth Haberm=
an and
David Maidment
University of Texas at Austin
April 2009
Contents
Introduction
Basics
of Kriging
Case
Study
Geostatistical
Analysis
Simple
Kriging and Inversed Distance Weight Methods.
To
be turned in
In this exe=
rcise 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 interpolati= ng a spatially continuous variable from point samples. For example, many spatial= ly 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 a= nd ordinary kriging.
The IDW is simple and intuitive deterministic interpola= tion 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 (hig= her 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 su= fficiently 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 spati= al autocorrelation. It uses a semivariogram.
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 dista= nce 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 spa= tial 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 implic= itly used in IDW. In kriging, one must = model the spatial autocorrelation using a semivariogram instead of assuming a dir= ect, linear relationship with separation distance.
Semivariance equals one-half the squared difference be= tween points separated by a distance d±∆d (assuming no direction preference). As the distance betwe= en samples increase, we expect the semivariance to also increase (again, becau= se near samples are more similar than distant samples). This is true, however, only up to some = given separation distance. For this dist= ance and up, points are unrelated. Stat= ed another way, if 50m is this critical separation distance, two points separa= ted 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.
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 semivariog= ram 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 mo= del first flattens out is known as the range Sample locations separated by distances closer than the range are spatially autocorrelated, whereas locat= ions 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 part= ial 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 ye= llow 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 interest= ed in the derivation of the weighting parameters (or any of the other topics discussed here), Applied Geostatistics by Edward H. Isaaks and R. Mo= han Srivastava is an excellent resource. Or for the more mathematical folks, try Statistics for Spatial Data by = Noel A.C. Cressie.
Goals of this Exe=
rcise
Now we know the basics of spatial interpolation. Let’s= use our knowledge to develop a digital elevation model of the Wax Lake Outlet A= rea using the 1963 bathymetric data.
Computer Requirem=
ents
This exercise is to be performed in ESRI ArcMAP 9.x wi= th 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 fro= m a 1963 bathymetry chart of the Wax Lake Outlet Area which is located in the d= elta of the Mississippi River in Louisiana, as shown in Figure 3.
Figure 3 := <= /span>from Wellner et al. 2005
The Wax Lake Outlet is a ma=
n-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 rout=
ed
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 Ou=
tlet
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 a= l. (2005)
We have uncovered a bathyme=
try
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 bathymetr=
ic
chart. This exercise will turn this chart into a digital elevation model.
Figure 4 : from USACE
Goals of this Exercise
Now we know the basics of spatial interpolation. Let’s use our knowledge to develop a di= gital 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 wi= th 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 fro= m the 1963 bathymetry chart of the Wax Lake Outlet Area. Additional data came from NHDPlus Regio= n 8
http://www.ho= rizon-systems.com/nhdplus/HSC-wthMS.php
Getting Started <= o:p>
In this exercise we need to use spatial analyst and
geospatial analyst toolbar. Open =
ArcGIS. We n=
eed to
enable this extension in the first place. Open Tools à
Extension.
Check the ‘Sp= atial analyst’ and ‘Geospatial analy= st’ box. Now we can use all the functions in both toolbars.
The digitized water depths from the bathymetric chart a= re contained in the spreadsheet Depth.xls and c= onsist 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, t= he latitudes and longitudes would have been referenced to the North American D= atum of 1927 (NAD27), the standard for earth datums = at that time.
These data were extracted from the Excel file and conve= rted to ArcGIS as point events, to become the feature class DepthEvents. This feature class was coupled with another for NHDWaterbody representatio= n to provide spatial context for the point set. The NHD dataset is presented in geographic coordinates with the North American Datum of 1983 (NAD83). Be= cause 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 Sta= te Plane Coordinate system for South Louisiana, using the NAD83 datum. The DepthEvents were first converted from NAD27 to NAD83 datum before this projection to St= ate Plane Coordinates was done. The re= sulting dataset is stored in the WaxLakeOutlet feature dataset in the Ex9 geoddatabase, as shown below= :
Open ArcGIS, and add the
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 Geospati= al 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 la= yer as DepthEvent= s 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 elev= ations 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 transfo= rm 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= =3DX, then the data follow a normal distribution.
Open Geospati=
al analyst à
Explore data à
The QQ plot shows the linear relationship between log(Depths_Events) and the standard normal distributio= n. 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 samp= les 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 Geospati= al analyst à Explore data àTrend analysis.=
The screen will look similar to the figure below. We can choose the different graph optio= ns 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 samp= les 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 Geospati=
al analyst à
Geospatial Wizard
Select the Kr=
iging
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.
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 Spheri= cal model for the interpolation
Other commonly used models are exponential and Gaussia= n. 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 statio= ns. Since we have the total 315 da= ta points, 12 and minimum 10 stations is an appropriate number for interpolati= on. Click on the third circle for sect= or type.
Click Next to proceed.
These are the error stat= istics for the interpolation. This is based on one-leave-out <= span style=3D'mso-spacerun:yes'> 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 l=
ight
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 th= e DEM, and the shallow areas that have been caused by sediment deposition.
Lets repeat the process and= get the Prediction Standard Error Map<= /p>
You can click on the map at any point in the interpola= tion 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.
Now go back to Geospatial analyst à Geospatial Wizard. Still choose Kriging and click next. This time choo= se Simple Kriging as the interpolati= on method. The rest is the same as Or= dinary Kriging procedures. You will get a result similar to the following (aft= er overlay the original data points on top):
Notice that this map agrees perfectly with the data po= ints 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 value= s at the sample sites. You can learn more about this method by reading the docum= ent provided by the software. Discuss in your homework whether this method is better than other two methods for this exercise.