Spatial Interpolation

Spatial Interpolation

General

The Spatial Interpolation function offers the possibility to make spatial interpolations on a defined grid. Two spatial interpolation functions are available, the Kriging and the Inverse Distance methods.
The function interpolates point values to a user defined grid. The point series can be selected from the Series list box and the grid characteristics are entered in the Grid Characteristics frame. Interpolations can be made on a grid defined by geographic (latitude/longitude) co-ordinates or on a local meter co-ordinate system. Make sure that the station co-ordinates are entered in the same co-ordinate system as the selected system in this function.

Under normal use a grid is made for 1 timestep only, therefore the option "Make Interpolation for first timestep only" must be checked. In this case the user is asked for an Arc/Info ASCII filename to save the results to. In case the Kriging function is selected, two filenames must be entered, one for the estimation and one for the variance.
HYMOS also computes the average value of the grid, this value is stored in memory and can be viewed with the View button. If the option "Make Interpolation for first timestep only" is not checked, HYMOS computes the average value for all selected timesteps. This resulting series can be stored in the HYMOS database by pressing the <Save> button.
When an average value of a catchment must be computed, make sure the "Blank grid with catchment boundary data" checkbox is checked an a catchment is selected from the list. In this case HYMOS sets all interpolated values outside the catchment area missing and computes only the average value of the non-missing values.
Results of the interpolation function can be viewed and analyzed by the graph program, the viewer, the report, or by the Netter program.

Series Codes

The series to use in the spatial interpolation can be selected or de-selected by clicking the checkbox of the series. By clicking your right mouse button on top of the Series list box you can select or de-select all series.

Grid Characteristics

The characteristics of the grid on which the interpolation must be made has to be defined. First a selection must be made on the co-ordinate system to make the interpolation on. When a 'geographic' co-ordinate system is selected all the station co-ordinates of the selected series must have been entered in latitude/longitude co-ordinates in the 'station Characteristics' form. When a 'local co-ordinate' system is selected, the local X-co-ordinate and local Y-co-ordinate of the stations must have been entered in the 'station Characteristics' form. This local co-ordinate system is a co-ordinate system in meters.

After selecting a co-ordinate system, the centre of the lower left corner of the grid must be entered. This lower left co-ordinate can also be computed from the station co-ordinates by pressing the <Fill Grid> button. HYMOS will compute the minimum X and minimum Y co-ordinate of the selected stations and enters this co-ordinate in the 'Lower Left Corner' cells of the spreadsheet. To complete the 'Grid Characteristics' information, the number of cells in X and Y direction and the 'Cell Size' must be entered. The cell size is entered in kilometres.

The above Grid Image contains the following information:

• the location of the lower left cell, this is the centre of the cell,
• the Cell Size, equal in X and Y direction (entered in km),
• the number of grid cells in X and Y direction

Catchment Blanking

When an average value of a catchment must be computed, check the "Blank grid with catchment boundary data" checkbox to activate the catchment list box. The catchment list box contains the ID's and names of all catchments stored in the HYMOS database. Make sure that a closed polygon is entered as catchment boundary and that the co-ordinate system of the catchment boundary is equal to the selected co-ordinate system.

Extra Options

Two extra options are given:

• Make Interpolation for first timestep only
• Data values as % of long time averages

When the first option is checked, HYMOS interpolates only for the first timestep of the selected processing period.
When the second option is selected, HYMOS divides the point values by their long time average value before making the spatial interpolations. This value must be entered for each point series in the "Series Characteristics" function.

Make Interpolation for first timestep only
Compute a interpolated surface for the first timestep only. When this button is not checked then for all timesteps of the selected period an interpolated surface will be computed

Kriging

By 'Kriging Technique' we denote a class of interpolators based on a stochastic approach.
There is a large array of Kriging techniques, sometimes obscured by cumbersome notations and mathematics. All these different Kriging techniques have a very simple and common line: they are all regression techniques and differ only by the particular set of functions of the data which are combined into an estimator. A unique advantage of the Kriging interpolation method is its ability to quantify the reliability of prediction, to provide an estimate plus a confidence interval.

Basic Idea of Kriging

The basic idea of the different Kriging methods is the same. If x denotes a point in space and Z(x) is a function of x which is known in the n observation points x₁, x₂, .., x_n, we look for an estimate Z^*(x₀) of Z(x₀) in an non-measured location of the form:

a linear combination of the observed values, supplied with a weight λ_i for each observation point.

The difference between estimate value Z^*(x₀) and the true value Z(x₀) for any set of weights is called the estimation error ε_est. If there is no trend (i.e. the stationary hypothesis holds) then Z^*(x₀) is an unbiased estimate of Z(x₀), the expected value of the difference between the estimated value and the true value will be zero:

The error variance, also called the estimation variance, of this difference, ε_est , is not zero: it is the mean of the squared difference. In the general case with the estimated point being a linear combination of values 1, 2, .., n, the error variance is as follows:

Each observation will contribute a different proportion to the total estimation variance of Z(x₀). There are an infinite number of ways in which the weights can be allocated, and each will produce a different estimation variance. Among these at least one combination of weights must produce a minimum estimation variance. It is this combination which kriging seeks to find. If the covariance function; the semivariance or the covariance is known, the weights λ_i can be calculated.

Ordinary Kriging

Ordinary Kriging is the most utilised type of Kriging, it can be used when the variable is stationary and the covariance is known, in contrary the mean is unknown. Ordinary Kriging refers to the following model assumption:

where m is the unknown stationary mean. The expected value of the error at any particular location if often referred to as the bias. Implementing the unbiasedness condition, equation on the linear estimation function, we can write:

This yields to the first constraint:

 

Unbiasedness is guaranteed because the coefficients sum to 1.
The second constraint is that the weights λ_i have to be solved with a minimal error variance. To find the set of weights that will give the minimum mean square error, the technique of Lagrangian multipliers is used. The optimal parameters satisfy:

where:
γ(x_i, x_j) = the semivariance of Z between the sampling points,
γ(x₀, x_i) = the semivariance between the sampling points x_i and the estimated point x₀ .

The reliability of interpolation (e.g. error variance) is found with the solutions of the weights:

The estimation error variance σ_k² can be regarded as depending exclusively on the number and the locations of the measurement locations. Therefore σ_k² is an efficient tool for solving network optimisation problems such as the optimal choice of measurement locations. It must be emphasised that σ_k² is not the variance of the actual real spatial estimation error but a modelled error. σ_k² Provides a theoretical measure of the relative accuracy's of the various estimates.

Kriging in HYMOS

In HYMOS the Kriging interpolations are performed according to ordinary Kriging. Different spherical models can be selected among which the Spherical, Exponential, Gaussian and Power model. For more information on these models see the Spatial correlation function.

While doing the interpolations the following default values are set, these can not be changed by the user:

• Search radius = 1E10
• Angle of Continuity = 0
• Anisotropy = 1
• Discretization of Blocks = 1*1
• Minimum number of samples = 4
• Maximum number of samples = Number of selected Stations

The figure shows a typical Kriging Variance map. The variance is low around the stations, this indicates that the estimation accuracy is high.

After the computation is done HYMOS asks the user to save the computed Grid files. For the Kriging interpolations the user must first enter a name for the interpolated values, the second file to save is the grid file with Kriging variances.

Inverse Distance

The basis of the inverse distance method is that the weights are inversely proportional to the distances. The weights may be raised to a power, not necessarily an integer value, to increase the effect of the weighting function. inverse distance squared is most commonly used:

where D_i, (D_j) is the distance between point x_i (x_j) and the point to be estimated, and p is the power.

The weights assigned to the measured data points depend only on a distance function. Delfiner and Delhomme made the following comment with regard to distance weighting schemes: 'Clearly, no general rule can be derived from experiment on particular data and point configurations. Consequently, the choice of a distance weighting function is more or less a matter of personal belief, of tradition or of confidence in the device of influential authorities'.
Distance weighting schemes suffer from a certain arbitrariness in the selection of interpolation parameters, also they do not cope well with clustered data. They are however much used for having a first impression of the spatial distribution of a measured parameter.

Notes:
From a scientifically point of view one must exercise extreme caution about likely accuracy of any interpolation scheme, including those of great mathematical complexity. This is because the validity of any interpolation scheme has to be seen against the extreme variability caused by non linearity of the phenomenon (i.e. the variable) under study.

Before using the Kriging function the phenomenon under study must first be analysed thoroughly. The semi-variogram must be determined ("Spatial correlation" function) and the parameter values of the semi-variogram entered.

The Kriging function is much more complex than the inverse distance function. However, when the phenomenon is studied and the semi-variogram is determined, this function is much more powerful than the inverse distance function. The advantages of Kriging are amongst others the declustering of measurement stations and a variance map with the accuracy's of the interpolation..

Application

The Spatial Interpolation function can be used with any kind of time series data, the only restriction is that no missing values are allowed in the time series. Computation procedure:

1. Select some series from the Select series list box.
2. Select the interpolation function to use and enter the function parameters (i.e. Kriging or inverse distance).
3. Enter the appropriate Grid Characteristics, explained in the Grid Characteristics section.
4. Press the <Execute> button to start the computation.
5. Enter a filename for the ASCII Grid file where the computed data must be stored to.
6. To show the results in Netter, select the "Options" - "Map Options" menu and load the ASCII file as a grid.

ASCII Raster File Format

The interpolated data are stored in an ASCII raster file format, also called the Arc/Info ASCII Grid file format. This is a simple and common format that can be used by many GIS, like ArcInfo and ArcView.
The format is very simple, basically a header followed by a list of cell values. The header that includes the following keywords and values:

The first row of the data is the top of the data set, moving from left to right. Cell values should be limited by spaces. The number of cell values must be equal to the number of rows times the number of columns.
Example:
377 372 366 362 357 353 348 344 340 336 ....
333 330 328 326 326 327 330 334 339 345 ....
......

Presentation of Grids in Netter

The Spatial Interpolation function produces two types of grid files, static and dynamic grid files. To produce static grid files, the option "Make Interpolation for first timestep only" must be checked. In this case Arc/Info ASCII files are saved containing interpolation results for only 1 timestep. To produce dynamic grid files the above mentioned option must be unchecked. In this case dynamic BIL files are saved, containing interpolated grids for more than 1 timestep.

To add a static grid file as a map layer in Netter follow these steps:
1. Select "Options" - "Map Options" from the Netter menu.
2. Press the <Add Layer>, top left button and select a ASCII grid file

3. Press the <OK> button

From the Map Options menu the grid map can be customised to the needs of the user. Some of the possible customisation options are:

• Smooth the contours with the "Data-Smooth contours" checkbox.
  Show Isolines with the "Data-Show Isolines" checkbox. The isolines can be customised (change line colours, with and frequency, or change label colours, font and frequency) by pressing the <Isolines> button.
• Show the datarange in the map legend by checking the "Show in legend" checkbox of the "Properties" tab. Then select "Option" - "Legend Options" and make the map legend visible including grid data.

• If you want to remove the grid colours, but want to keep the isolines select: "Options" - "Map Options" - "Data-tab" - Options" and let the data range (Lower bound and Upper bound" of the "Missing Value" option cover your whole data range.

• You can also use elevation maps as background map layers and show these maps using the shading option. The GTOPO30 maps, distributed by the USGS can be loaded directly in Netter.