13.2 Surface Analysis & Spatial Interpolation: IDW, Kriging, Spline & Viewshed Analysis

13.2 Surface Analysis & Spatial Interpolation: IDW, Kriging, Spline & Viewshed Analysis

Quick Summary: Continuous geographic fields—such as elevation, barometric pressure, precipitation, and soil contaminant concentrations—cannot be exhaustively measured at every point across a landscape. Geospatial analysts employ spatial interpolation to estimate unknown values at unmeasured locations from discrete sample point observations. Grounded in Tobler's First Law of Geography, interpolation algorithms span deterministic techniques (Inverse Distance Weighting, Spline) that rely purely on geometric proximity and surface smoothness, and geostatistical methods (Kriging) that model spatial autocorrelation using semivariograms to generate optimal estimates alongside quantified error variance surfaces. In addition, terrain visibility modeling through Viewshed Analysis evaluates line-of-sight intervisibility while accounting for observer offsets, target heights, Earth curvature, and atmospheric refraction.

---

1. Theoretical Foundations: Tobler's First Law & Surface Modeling

Spatial interpolation rests on the foundational axiom articulated by Waldo Tobler in 1970:

"Everything is related to everything else, but near things are more related than distant things." — Waldo Tobler

This principle, known as Tobler's First Law of Geography, implies that geographic space possesses spatial continuity and spatial dependence (spatial autocorrelation). If near things were not more related than distant things, spatial interpolation would be impossible, and unmeasured locations could only be estimated using a global mean.

Classification of Spatial Interpolators

Interpolation techniques are formally classified across three primary structural dimensions:

1. Deterministic vs. Geostatistical:
   • Deterministic: Uses mathematical functions based on Euclidean distance or geometric smoothness (e.g., IDW, Spline, Trend Surface). They do not evaluate statistical autocorrelation and cannot quantify prediction standard errors.
   • Geostatistical: Leverages statistical and probabilistic models of spatial autocorrelation (Kriging). They provide both a predicted value surface and a statistical measure of uncertainty (prediction error surface).
2. Exact vs. Inexact:
   • Exact: The interpolated surface passes precisely through every sampled observation point (e.g., IDW, Thin-Plate Spline, Kriging with zero nugget). At sample point $s_i$, $\hat{Z}(s_i) = Z(s_i)$.
   • Inexact: The surface smooths through the data without necessarily matching sample values exactly (e.g., Global Trend Surface, Kriging with a measurement error nugget).
3. Local vs. Global:
   • Local: Computes estimates using only neighboring points within a specified search radius or neighbor count (IDW, Spline, Ordinary Kriging).
   • Global: Uses every sample point across the entire dataset simultaneously to fit a single overarching mathematical equation (Trend Surface Polynomials).

---

2. Inverse Distance Weighting (IDW)

Inverse Distance Weighting (IDW) is the most common deterministic interpolation method. It operates under the premise that the attribute value at an unmeasured location $\hat{Z}(s_0)$ is a distance-weighted average of neighboring sample values.

$\hat{Z}(s_0) = \frac{\sum_{i=1}^{n} w_i Z(s_i)}{\sum_{i=1}^{n} w_i}, \quad \text{where } w_i = \frac{1}{d_i^p}$

Here, $Z(s_i)$ is the measured value at point $i$, $d_i$ is the Euclidean distance from sample point $i$ to the prediction location $s_0$, and $p$ is a positive real exponent termed the power parameter.

   Sample P1 (Z=100)               Sample P2 (Z=20)          IDW PREDICTION AT P0:
          *                               *                  Distance to P1: d1 = 20 m
           \                             /                   Distance to P2: d2 = 80 m
            \           P0 (?)          /                    With Power p = 2:
             \            *            /                     w1 = 1 / (20^2) = 0.0025
              \                       /                      w2 = 1 / (80^2) = 0.000156
               \                     /                       
                +-------------------+                        P0 is heavily weighted toward P1!

The Power Parameter ($p$)

The power parameter dictates the rate at which sample influence decays with increasing distance:

• Standard Default ($p = 2$): Known as inverse-distance-squared weighting. Weights decrease rapidly with distance, giving local samples strong influence while maintaining a smooth transition.
• Low Power ($p \to 0$ or $p = 1$): Weight diminishes slowly. Distant sample points exert substantial influence, pulling estimates toward the regional average and producing a flattened, over-smoothed surface.
• High Power ($p \ge 3$ or $p = 8$): Weight drops precipitously with distance. The nearest single observation completely dominates the prediction, causing the surface to resemble stepped polygonal plateaus similar to Voronoi/Thiessen polygons.

Critical Characteristics and Limitations of IDW

1. Exact Interpolator: When predicting at a coordinate that coincides with an existing sample point ($d_i = 0$), the weight approaches infinity, returning the exact measured sample value.
2. Bounded Range (No Extrapolation): An IDW surface cannot predict values lower than the minimum sample value or higher than the maximum sample value. It cannot model natural mountain peaks or river trenches that were not directly captured in the sample dataset.
3. The "Bullseye" Effect: Around isolated sample points that diverge from local averages, IDW produces concentric circular contour rings resembling bullseyes or donuts. This artifact occurs because distance weights decay symmetrically in all directions from the isolated point.
4. No Error Estimation: IDW provides no mathematical mechanism to calculate confidence intervals, standard errors, or prediction uncertainty.

---

3. Spline & Radial Basis Functions

Spline interpolation fits a continuous mathematical function that passes through sample points while minimizing total surface curvature. The mathematical analogy is bending a thin, flexible sheet of metal or rubber through fixed elevation pins.

$\min \int \left( \frac{\partial^2 f}{\partial x^2} + 2\frac{\partial^2 f}{\partial x \partial y} + \frac{\partial^2 f}{\partial y^2} \right)^2 dx\, dy$

Spline Formulations: Regularized vs. Tension

• Regularized Spline: Minimizes higher-order derivatives, producing a very smooth, gently curving surface. It is well suited for phenomena with gradual, continuous spatial transitions such as barometric pressure or regional groundwater tables.
• Tension Spline (Spline with Tension): Adds a tension parameter that increases the stiffness or rigidity of the simulated sheet. Higher tension forces the surface to transition more abruptly between sample points, suppressing wild oscillations.

   REGULARIZED SPLINE                          SPLINE WITH TENSION
   (Smooth, allows gentle curves)             (Stiffer, restrains extreme peaks/overshoots)
            ^                                            ^ 
           / \                                          / \
   *------/---\------*                          *------*---*------*
  P1     /     \    P2                         P1     P_peak    P2
        v       v
   (Risk of extreme overshoot/undershoot!)     (Constrained closer to sample bounds)

Analytical Trade-offs of Splines

• Extrapolation Capability: Unlike IDW, Splines can predict values above the sample maximum and below the sample minimum, allowing them to model realistic ridge crests and valley bottoms.
• Overshoot / Undershoot Hazard: In areas with steep localized gradients (e.g., a high-elevation point immediately adjacent to a low-elevation point), Regularized Spline surfaces often overshoot wildly, predicting impossible negative elevations or extreme artificial spires.

---

4. Geostatistical Interpolation: Kriging & Semivariogram Modeling

Developed by South African mining engineer Danie Krige and formalized by French mathematician Georges Matheron, Kriging is an advanced geostatistical interpolation method. Unlike deterministic approaches, Kriging models the statistical spatial correlation among sample points to generate Best Linear Unbiased Estimates (BLUE) alongside a continuous prediction standard error surface.

The Empirical Semivariogram

The foundation of Kriging is the semivariogram, which quantifies spatial autocorrelation by calculating the average dissimilarity (semivariance $\gamma$) between pairs of sample points separated by a specific geographic distance (lag distance $h$):

$\gamma(h) = \frac{1}{2 N(h)} \sum_{i=1}^{N(h)} \left[ Z(s_i) - Z(s_i + h) \right]^2$

Where $N(h)$ is the number of sample point pairs separated by lag distance $h$.

   Semivariance γ(h)
      ^
      |                             SILL (Total Variance)
      |                       -----------------------------
      |                     / |
      |                   /   |
      |                 /     | <--- PARTIAL SILL (Sill - Nugget)
      |               /       |
      |              /        |
      |      +------/         |
      |      |                |
      | NUGGET                |
      +------+----------------+-----------------------------> Lag Distance (h)
             0              RANGE
      (Microscale error)   (Distance beyond which points are uncorrelated)

Anatomy of the Semivariogram

1. Nugget: The semivariance value where the curve intercepts the vertical axis at lag $h = 0$. In theory, points at distance zero should have zero variance. A positive nugget reflects measurement error or microscale spatial variation occurring at distances smaller than the sample spacing.
2. Range: The lag distance at which the semivariogram curve levels off into a horizontal plateau. Sample points separated by distances less than the range are spatially autocorrelated; points separated by distances greater than the range are spatially independent and share no spatial autocorrelation.
3. Sill: The plateau value of semivariance reached at the range. The sill represents the total sample variance of the dataset.
4. Partial Sill: The difference between the Sill and the Nugget ($\text{Partial Sill} = \text{Sill} - \text{Nugget}$). It represents the proportion of total variance explained by spatial autocorrelation.

Fitting Theoretical Models

To compute Kriging weights, an analyst must fit a continuous theoretical mathematical model to the discrete empirical semivariogram points:

• Spherical Model: Rises progressively and flattens out completely at the range. Widely used for geological and environmental modeling.
• Exponential Model: Rises steeply and approaches the sill asymptotically (effective range is defined as 95% of the sill).
• Gaussian Model: Exhibits an S-shaped curve with gentle initial slope, representing extremely continuous and smooth physical processes.

Kriging Methodologies

---

5. Summary Comparison of Spatial Interpolation Methods

---

6. Viewshed & Line-of-Sight (LOS) Analysis

Viewshed analysis calculates the intervisibility of geographic terrain between one or more observer locations and all cells on a Digital Elevation Model (DEM). It returns a binary or cumulative raster identifying which cells are visible from the observer points.

   Observer (Z + observer-height offset)
      o
     / \            Line of Sight (LOS)
    /   \--------------------------------- - - - - - - - (Obstructed Target)
   /     \                                     X
  /  Hill 1 (Blocks LOS)                      / \
 /         \                                 /   \
+-----------+-------------------------------+-----+------------------ DEM Surface

Critical Viewshed Parameters

In GIS viewshed engines, visibility is parameterized using specific attribute fields assigned to observer points:

• Observer-height offset: Vertical height added to the observer point's surface elevation (e.g., a 1.75-meter human eye level or a 30-meter communications tower).
• Target-height offset: Vertical height added above the terrain surface at each target cell being evaluated (e.g., evaluating visibility of a 100-meter wind turbine blade tip versus ground-level smoke).
• AZIMUTH1 & AZIMUTH2: Defines horizontal scan limits (e.g., restricting field of view between 90° East and 180° South).
• VERT1 & VERT2: Defines upper and lower vertical viewing angles in degrees.
• RADIUS1 & RADIUS2: Establishes minimum inner and maximum outer visibility search distances.

Earth Curvature & Atmospheric Refraction Corrections

Over long distances ($> 3\text{--}5,\text{km}$), two physical factors alter line-of-sight visibility:

1. Earth Curvature: The Earth's surface curves downward away from the horizontal line of sight at a rate proportional to the square of the distance. An object 10 km away drops approximately 7.85 meters below the flat tangent plane.
2. Atmospheric Refraction: Variation in air density bends light rays downward toward the Earth, partially counteracting the curvature effect (restoring approximately 13% of apparent elevation).

The combined geometric correction formula is:

$\Delta Z = \frac{D^2}{2R_E} (1 - k)$

Where $D$ is the distance to the target, $R_E$ is the mean radius of the Earth ($\approx 6,371,000,\text{m}$), and $k$ is the atmospheric refraction coefficient (typically taken as $0.13$). Failing to enable curvature and refraction corrections during long-range viewshed modeling severely overestimates visibility across open plains and maritime environments.

---

7. Practical Geospatial Scenario: Industrial Contamination & Communication Tower Siting

Scenario A: Regulatory Groundwater Monitoring

An environmental consultancy must map trichloroethylene (TCE) concentrations in an aquifer based on 32 monitoring wells. The state environmental agency requires not only an estimated concentration map but also an explicit statistical confidence interval surface to justify the placement of remediation extraction wells.

• Decision: The analyst selects Ordinary Kriging rather than IDW. Because Kriging evaluates the semivariogram, it delivers both the predicted TCE concentration grid and the prediction standard error surface, demonstrating to regulators where data uncertainty is high and where additional monitoring wells are required.

Scenario B: Emergency Dispatch Tower Viewshed

A county sheriff's department must place a new 50-meter radio repeater tower to ensure line-of-sight communications with emergency patrol vehicles across a rugged valley:

• Parameter Configuration: The observer feature is set with observer-height offset = 50.0 (tower height). The target vehicle antennas are modeled with target-height offset = 2.0 (two meters above terrain). The analyst enables Earth curvature and atmospheric refraction corrections because line-of-sight paths extend 28 kilometers to the county perimeter.

---

8. Common Exam Traps & Pitfalls

[!CAUTION] Exam Trap 13.2.1: IDW Cannot Extrapolate Beyond Sample Min/Max. A standard exam question asks: "A GIS technician interpolates high-elevation mountain peaks using Inverse Distance Weighting from a sample of valley weather stations. Why does the resulting model fail to show the true mountain summit elevations?" The answer is that IDW is a distance-weighted average; it can never predict a value greater than the maximum sample value or lower than the minimum sample value.

[!CAUTION] Exam Trap 13.2.2: Kriging vs. Deterministic Error Surfaces. Questions frequently ask which interpolation algorithm provides an accompanying measure of prediction uncertainty or standard error. The correct answer is Kriging. Deterministic methods (IDW, Spline, Natural Neighbor) are purely geometric and cannot produce a prediction error variance surface.

[!CAUTION] Exam Trap 13.2.3: observer-height offset vs. target-height offset in Viewshed Analysis. Memorize the distinction: observer-height offset is the vertical height of the observer added to the surface DEM. target-height offset is the vertical height of the target being viewed. Setting observer-height offset = 100 and target-height offset = 0 models a 100-meter observer looking at ground-level targets, which produces a completely different visibility pattern than an observer at ground level looking at a 100-meter structure.