Spatial Analysis

Part 1: Concept, Purpose, and DEM Generation

1. What is Spatial Interpolation?

Spatial interpolation is the mathematical and statistical process of estimating the values of a continuous geographic phenomenon at unknown, unsampled locations based on the known values of surrounding sampled points.

2. What is the purpose of spatial interpolation?

In reality, it is impossible, physically and financially, to measure data (like elevation, rainfall, or soil strength) at every single square meter of a study area. The primary purpose of interpolation is to take isolated sample points (like weather stations, borehole logs, or GPS spot heights) and create a continuous surface (a raster grid) that models the entire region.

Comparison: IDW vs. Minimum Curvature Spline

3. Why is interpolation needed for producing a DEM?

A Digital Elevation Model (DEM) is a continuous raster grid where every single pixel contains an elevation value. However, field surveyors and LiDAR point clouds only capture elevation at discrete, scattered points. Interpolation is absolutely strictly required to calculate the elevation of the "empty" spaces between the surveyed points, effectively creating a solid, continuous topographic surface from isolated data.

Part 2: Classification and Types of Interpolation

If asked to classify or explain the different types break them down into two main families:

A. Deterministic Methods

These use predefined mathematical equations based on the distance or geometric layout of the known points. They do not account for statistical probability.

• Inverse Distance Weighted (IDW): Assigns values to unknown points by using a weighted average of the surrounding known points. The closer a known point is the more heavily it influences the unknown value.
• Spline (Minimum Curvature): Fits a mathematically smooth flexible "rubber sheet" through all the known data points minimizing the overall curvature of the surface.
• Thiessen Polygons (Nearest Neighbor): Creates sharp boundaries. Any unknown location simply takes the exact value of the single nearest known point.

B. Geostatistical Methods

These rely on both the distance and the statistical relationship (spatial autocorrelation) among the measured points.

• Kriging: A highly advanced method that not only estimates the unknown values but also provides a map of the estimated error/uncertainty. It looks at the overall spatial trend and variance in the dataset.

Part 3: Detailed Method Comparison (IDW vs. Spline)

The examiners frequently ask you to describe one method in detail discuss its pros and cons or compare IDW and Spline directly.

Inverse Distance Weighted (IDW)

• Concept: Assumes that "things that are closer are more alike."
• Advantages: Very simple to understand fast to compute and respects the exact values of the sampled points.
• Disadvantages: It cannot predict a value higher than the highest known point or lower than the lowest known point (no extrapolation). It often creates artificial "bullseye" or circular patterns around isolated data points.

Spline (Minimum Curvature)

• Concept: Fits a smooth curve through the points much like bending a flexible ruler.
• Advantages: Creates highly visually appealing smooth surfaces. Unlike IDW a Spline can predict values outside the range of the input data (e.g. predicting a mountain peak even if surveyors only measured points on the slopes).
• Disadvantages: Highly sensitive to outliers. If one data point is incorrectly recorded as extremely high the Spline will wildly distort the surrounding surface to force the "rubber sheet" through that error.

Part 4: Engineering Applications

To answer "application in different engineering fields" effectively provide these specific examples:

1. Hydrology (mHM Modeling): Interpolating isolated daily rainfall data from a few regional meteorological stations to create a continuous precipitation raster across an entire river basin.
2. Geotechnical Engineering (Site Response Analysis): Interpolating shear wave velocity (V_s30) or bedrock depth from scattered borehole logs to create a continuous subsurface vulnerability map for the Kathmandu Valley.
3. Water Resources (HEC-RAS): Interpolating riverbed bathymetry from surveyed cross-sections to generate a continuous 3D channel terrain for 2D flood inundation modeling.
4. Environmental Engineering: Modeling the dispersion plume of airborne pollutants or groundwater contamination spreading from an industrial site based on a few scattered air/water quality monitoring sensors.

Part 5: The Numerical Question

Question: Interpolate a value using any one method if the distance from four known values 160 170 180 and 90 are 1000 400 600 and 120 m respectively.

Always use the Inverse Distance Weighted (IDW) method for manual exams assuming a power of p=1 (or p=2 but p=1 is faster to calculate manually unless specified). Let's use p=1.

This is the linear weighting method we did earlier. It assumes influence drops off linearly with distance.

Inverse Distance Weighted (IDW) with Power = 2

This is the default setting in most GIS software (like QGIS or ArcGIS). Squaring the distance means that points further away lose their influence exponentially faster.

Notice how changing the power to 2 pulled the final answer much closer to 90? That is because squaring the distance made the far-away points (like the one 1000m away) almost mathematically irrelevant.

7.2 Raster Analysis Including Topological Overlay

Part 1: Theoretical Foundations of Raster Analysis

1. What is Raster Analysis?

Raster analysis is the process of performing spatial operations on raster (grid-based) data to derive new information. Unlike vector analysis which works with discrete features, raster analysis processes continuous surfaces where each cell holds a numeric value.

Advantages and Disadvantages of Raster-Based Analysis

Map Overlay vs. Map Calculations

2. Types of Raster Operations (Local, Focal, Zonal, Global)

Raster operations are classified by the spatial extent of cells involved in each calculation:

Why NDVI is Local and Slope is Focal: NDVI is a local function because it calculates each pixel's value using only the corresponding NIR and Red band values at that same pixel location. Slope is a focal (neighborhood) function because calculating the slope at any cell requires reading the elevation values of its 8 surrounding neighbors in a 3×3 moving window.

Map Algebra (Tomlin's Framework)

Map algebra is a set of rules and operators for performing cell-by-cell mathematical operations on raster data, first formalized by Dana Tomlin. It is the raster equivalent of SQL for vector data.

• Arithmetic operations: Addition, subtraction, multiplication, division of raster layers (e.g., Weighted Overlay = 0.5 × Slope + 0.3 × Soil + 0.2 × LandUse)
• Boolean operations: AND, OR, NOT — used for suitability queries (e.g., Slope < 10 AND LandUse = "Grassland")
• Relational operations: Greater than, less than, equal to — create binary rasters (1/0)
• Conditional operations: IF-THEN-ELSE — reclassify based on conditions

Reclassification

Reclassify is the process of reassigning cell values in a raster to new values based on defined criteria. It converts continuous data into discrete categories for analysis.

• Example: Reclassify a slope raster: 0–5° → 1 (Flat), 5–15° → 2 (Moderate), 15–30° → 3 (Steep), >30° → 4 (Very Steep)
• Purpose: Simplifies complex data for overlay analysis. Essential in suitability modeling where continuous data must be converted to a common ranking scale (e.g., 1–10).
• Types:
  • Range-based: Group values into ranges (0–100 → Low, 100–200 → Medium)
  • Individual value: Replace specific values (NoData → 0)

Hillshade

Hillshade creates a 3D shaded relief representation of a terrain surface by simulating the illumination of a surface from a hypothetical light source (the sun). It is a focal function that uses a 3×3 neighborhood window.

• Input: DEM (Digital Elevation Model)
• Parameters:
  • Altitude (Sun Elevation Angle): Height of the sun above the horizon (default: 45°)
  • Azimuth (Sun Direction): Compass direction of the sun (default: 315° = northwest)
• Formula: Hillshade = 255 × [(cos(Zenith) × cos(Slope)) + (sin(Zenith) × sin(Slope) × cos(Azimuth − Aspect))]
• Output: Grayscale raster (0–255) where bright areas face the sun and dark areas are in shadow
• Application: Visual enhancement of terrain maps; used as a background layer beneath semi-transparent thematic layers for better 3D perception

Viewshed (View Shade)

Viewshed analysis determines which areas on a terrain surface are visible from one or more observation points. It uses line-of-sight calculations across a DEM.

• Input: DEM + one or more observer points (with optional height offset)
• Output: Binary raster where 1 = visible from the observer, 0 = not visible (hidden by terrain)
• How it works: For each cell in the DEM, a straight line is drawn from the observer to that cell. If any intermediate cell's elevation blocks the line-of-sight, the target cell is marked as not visible.
• Applications:
  • Locating telecom towers to maximize signal coverage
  • Road alignment to identify scenic routes
  • Military: Finding hidden positions not visible from enemy locations
  • Environmental: Visual impact assessment of proposed buildings or wind turbines

Euclidean Distance

Euclidean Distance calculates the straight-line distance from each raster cell to the nearest source feature. It is a global function.

• Input: Source raster or feature (e.g., river, road, hospital locations)
• Output: Distance raster where each cell value = its distance to the nearest source
• Application: Finding areas within a certain distance of a river (for flood analysis) or determining accessibility to hospitals

Cost Path Analysis

Cost Path finds the least-cost route between a source and destination across a surface where movement cost varies by location.

• Input:
  • Cost Surface: Raster where each cell value represents the cost of traversing that cell (based on slope, land cover, etc.)
  • Source point: Starting location
  • Destination point: Ending location
• Process: (1) Create a Cost Distance raster showing accumulated cost from source, (2) Create a Back Link raster recording direction of cheapest path, (3) Trace back from destination to source using the back link
• Application: Finding optimal road alignment that minimizes construction cost across varying terrain; pipeline routing; wildlife corridor planning

• Concept: This function calculates summary statistics (like Area, Sum, Mean, Maximum, Minimum) for a continuous "Value Raster" strictly within the boundaries defined by a separate "Zone Dataset". The Zone Dataset can be either a vector polygon shapefile or an integer raster.
• Mechanism: The software identifies all the pixels of the Value Raster that physically fall inside Zone A, lumps them together, calculates the requested statistic (e.g., the average), and then moves on to Zone B.
• Exam Example: For your final year project on Site Response Analysis, you have a highly detailed, continuous raster grid showing Peak Ground Acceleration (PGA) across the Kathmandu Valley. You also have a polygon shapefile of the municipal wards. Using the Summarize Zone tool, the software will calculate and output a table showing the exact Mean and Maximum PGA specifically for each individual ward, allowing you to rank the wards by seismic vulnerability.

3. Proximity Map Calculation & Distance

Distance and Proximity are often used interchangeably in raster analysis, but they represent slightly different ways of modeling "how close" something is.

Distance (Euclidean Distance)

• Concept: Calculates the straight-line distance from every single cell center to the closest source feature.
• Mechanism: It generates a continuous raster where the numerical value of every single pixel represents its distance (in meters or feet) to the target.
• Exam Example: When mapping the "Smart Village" initiative, you need a general overview of water accessibility. You run a Distance function on the "Water Tap" point layer. The output is a gradient raster where pixels right next to the taps have values near 0m, and pixels on the outskirts of the village have values of 500m or more.

Proximity (Raster Buffering / Allocation)

• Concept: Instead of a continuous distance gradient, a Proximity calculation often groups cells into distinct bands or allocates them to the nearest source.
• Mechanism: It classifies pixels based on a threshold. For example, assigning a value of "1" to all pixels within 30 meters of a source, and "0" to everything else (acting like a vector buffer, but in grid format).
• Exam Example: To enforce sanitation regulations for a new rural toilet, it must be located at least 30 meters away from any drinking water well. You run a Proximity map calculation on the wells, instructing the tool to flag all raster cells within a 30m radius. This instantly generates a restricted "no-build" pixel zone across your study area.

4. Comprehensive Raster Surface & Spatial Functions

(A) Surface

A surface in raster GIS is a continuous geographic phenomenon represented as a grid where every cell stores a measured or modelled value of a variable that varies smoothly across space (elevation, temperature, rainfall, pollution concentration, slope, etc.). The most common surface is a Digital Elevation Model (DEM), where every pixel stores a Z value (elevation in metres).

• Why "surface"? If you imagine the Z value as the height above the X–Y plane, the raster looks like a 3D bumpy sheet — hence the term surface.
• Derivatives of a surface: Many other rasters are mathematically derived from a surface DEM — Slope, Aspect, Hillshade, Curvature, Contour, Flow Direction, Flow Accumulation, Viewshed.
• Engineering Example: A DEM of the Kathmandu Valley generated from contour interpolation is the master surface from which a civil engineer derives slope (for landslide hazard), aspect (for solar exposure of buildings), hillshade (for visualisation in a project report) and viewshed (to select a microwave tower location).

(B) Slope

Slope is a focal (neighbourhood) function that measures the maximum rate of change in elevation between a cell and its eight surrounding neighbours in a 3×3 moving window. It is reported in degrees (0°–90°) or percent rise.

Formula (Horn’s method, used by ArcGIS):

where a…i are the elevations of the 3×3 neighbourhood (a = NW, b = N, c = NE, d = W, e = centre, f = E, g = SW, h = S, i = SE).

• Engineering Example: A road alignment in Sindhupalchok must avoid slopes > 30° to limit landslide risk. The Slope tool is run on the SRTM DEM, the output is reclassified to a binary mask (slope ≤ 30° = 1, slope > 30° = 0), and the road centreline is restricted to the “1” cells.

(C) Aspect

Aspect is also a focal function derived from the same 3×3 elevation neighbourhood. It identifies the compass direction (in degrees clockwise from north, 0°–360°) that the downhill slope of each cell faces. Flat areas (slope = 0) are coded as −1.

Formula:

The result is then converted to compass bearing: 0° = North, 90° = East, 180° = South, 270° = West.

• Engineering Example: For a solar PV farm site selection in Mustang, the design team wants south-facing slopes (135°–225°) that receive maximum daily insolation. The Aspect raster is reclassified into 9 classes (N, NE, E, SE, S, SW, W, NW, Flat) and only cells classified as “S” are kept for further suitability scoring.

(D) Raster Buffer

A raster buffer creates a zone of cells within a specified distance of one or more source cells. Unlike the smooth polygon produced by a vector buffer, a raster buffer is a stepped block of pixels whose resolution is dictated by the cell size. It is implemented as a Boolean Euclidean-distance threshold: buffer = (EuclideanDistance ≤ D).

• Mechanism: (1) Identify source cells (e.g., river pixels = 1, others = NoData). (2) Run Euclidean Distance → produces a continuous distance raster. (3) Reclassify: cells with distance ≤ D get value 1 (inside buffer), all others get 0.
• Engineering Example: The Department of Roads requires that no quarry be permitted within 200 m of the Trishuli River. The river pixels are buffered using a Euclidean Distance ≤ 200 m rule, producing a binary raster where the “1” cells form a no-quarry corridor along both banks.

(E) Raster Intersection & Overlay

In raster GIS, Intersection is the cell-by-cell logical-AND of two binary rasters: an output cell is “1” only when both input cells are “1”. Overlay is the more general process of combining two or more rasters using map-algebra (logical, arithmetic, or weighted-sum) to derive a single output raster.

• Boolean Intersection: Out = (A == 1) AND (B == 1) — classic suitability filter.
• Arithmetic Overlay: Out = A + B or Out = w1·A + w2·B + … — produces a continuous suitability score.
• Engineering Example: A landfill must satisfy both “slope < 5°” AND “land-use = barren”. Each criterion is reclassified to a binary raster, and the cell-by-cell AND (intersection) yields the candidate cells. A weighted overlay (0.5×Slope + 0.3×Soil + 0.2×LandUse) then ranks them.

(F) Local Operations — Mean, Sum, Min, Max

A local operation works on each cell independently, using the corresponding cell from one or more input rasters to compute the output cell value. No neighbour cells are consulted — each pixel is processed in isolation.

• Single-raster local: e.g., Out = ln(A), Out = A × 2.
• Multi-raster local: e.g., Out = (A + B + C) / 3 — the mean of three rasters cell by cell.
• Statistics available: Sum, Mean, Min, Max, Range, Standard Deviation, Majority, Variety.
• Engineering Example: Three monthly rainfall rasters (Jun, Jul, Aug) of the Bagmati basin are stacked. Cell-Statistics “Sum” gives total monsoon rainfall per pixel; “Mean” gives the average monthly rainfall — both used for runoff modelling.

(G) Focal Statistics

A focal (neighbourhood) operation computes each output cell from the values inside a moving window (typically 3×3, 5×5, or a circular kernel) centred on the corresponding input cell. The window walks across the raster pixel by pixel.

• Statistics available: Mean, Sum, Min, Max, Range, Std-Dev, Median, Majority, Variety.
• Common uses: smoothing/filtering noisy DEMs (3×3 mean), edge detection (max−min), removing speckle from satellite images.
• Engineering Example: A LiDAR DEM of a building site contains noisy spikes from vehicles. A 3×3 focal-mean filter is applied to produce a smoothed DEM that is then used for cut-and-fill volume estimation.

(H) Zonal Statistics

Zonal statistics summarises the values of a continuous “value raster” within the boundaries of a discrete “zone raster” (or polygon). The zone raster groups cells that share the same integer ID; the tool reports Sum, Mean, Min, Max, Range, Std-Dev, Variety, Majority, Minority for each zone.

• Inputs: Zone dataset (e.g., ward boundaries, watershed IDs) + Value raster (e.g., elevation, rainfall, NDVI, PGA).
• Output: A statistics table (one row per zone) or a raster where every cell of a zone is replaced by the zonal statistic.
• Engineering Example: The mean Peak Ground Acceleration (PGA) is computed for each ward of Kathmandu Metropolitan City. Wards 16, 19 and 32 emerge with mean PGA > 0.45 g and are flagged as priority for seismic retrofitting.

Consolidated Summary — All Raster Functions at a Glance

Scenario-Based Questions

1. The Environmental Impact Scenario

Question: You are working on a project to analyze the impact of deforestation on local water resources. What types of spatial data would you use, and how would you represent them in GIS? Describe with sketches.

This question aligns perfectly with soil conservation and watershed management studies. You need to list the exact data types and how the GIS software handles them.

Data Types & Representation:

1. Land Cover / Land Use (Raster): Used to represent the changing extent of the forest over time. A raster grid where pixels are classified as "Forest" or "Deforested/Barren".
2. Digital Elevation Model / DEM (Raster): A continuous elevation grid used to delineate the watershed boundaries and calculate terrain slopes, which directly affect surface runoff and soil erosion.
3. River Network (Vector – Line): A polyline shapefile representing the streams and rivers to analyze where the increased runoff and sediment load will ultimately flow.
4. Meteorological Data / Rainfall (Raster): An interpolated continuous surface (using IDW) of precipitation data to model the water input into the deforested basin.
5. Soil Hydrologic Groups (Vector – Polygon): Polygons representing different soil types to estimate infiltration rates vs. surface runoff.

The Sketch: Draw a simple flowchart showing: DEM + Land Cover + Soil Type pointing into a Hydrologic Model (e.g., mHM or HEC-HMS), which then points to a River Network Line showing "Increased Peak Discharge."

2. The Reservoir Site Selection Scenario

Question: You are tasked with identifying suitable locations for building a new reservoir in a region. The ideal location must meet the following criteria:

1. The area must have a slope of less than 10 degrees to ensure stability.
2. It must be within 1 kilometer of a river to allow for water supply.
3. The land cover must be non-forested (e.g., grassland or bare soil) to minimize environmental impact.
4. The area must have a high water-retention capacity, indicated by soil types classified as "clay" or "silt".

The data provided are DEM, Soil type, River, and Land cover data. Using raster data analysis, describe the steps you would take to identify suitable locations. Include at least four raster analysis functions in your explanation.

Approach: This is a classic Boolean Overlay (True/False) suitability model. Each dataset is independently processed, reclassified into binary (1=Suitable, 0=Unsuitable), and then multiplied cell-by-cell. The multiplication acts as a logical AND — a cell must pass all four tests to be marked suitable.

The Four Raster Analysis Functions Used

1. Slope (from DEM): A Focal/Neighborhood function that computes the rate of elevation change using a 3×3 moving window. Input: DEM raster. Output: Continuous slope raster in degrees (0°–90°). This is a neighborhood function because each pixel's slope depends on its 8 surrounding neighbors.
2. Euclidean Distance (from River): A Distance function that calculates straight-line distance from every cell center to the nearest river feature. Input: River vector (or raster). Output: Continuous distance raster in meters. Cells near rivers have low values; distant cells have high values.
3. Reclassify: A Local function that reassigns ranges of continuous values into discrete binary classes. Applied independently to all four datasets:
   • Slope: 0°–10° → 1, >10° → 0
   • Distance: 0–1000 m → 1, >1000 m → 0
   • Land Cover: Grassland/Bare soil → 1, Forest → 0
   • Soil Type: Clay/Silt → 1, Sand/Gravel → 0
4. Raster Calculator / Map Algebra: A Local function that performs cell-by-cell arithmetic across multiple aligned rasters. Here we multiply all four binary grids: Final = Slope_Rec × Distance_Rec × LandCover_Rec × Soil_Rec. Any cell with a 0 in any layer becomes 0 in the output (Boolean AND logic).

Workflow Overview:

Step 1: Calculate Slope (Tool: Slope)

Input the DEM raster and run the Slope function to generate a continuous raster showing the steepness in degrees.

Step 2: Calculate Proximity (Tool: Euclidean Distance)

Input the River vector data and run the Euclidean Distance tool to create a raster where every pixel value represents its distance from the river.

Step 3: Standardize the Criteria (Tool: Reclassify)

You must convert all datasets into a binary scale where 1 = Suitable and 0 = Unsuitable.

• Reclassify Slope: 0° to 10° = 1. Greater than 10° = 0.
• Reclassify Distance: 0 to 1000m = 1. Greater than 1000m = 0.
• Reclassify Land Cover: Grassland/Bare soil = 1. Forested = 0.
• Reclassify Soil: Clay/Silt = 1. Sand/Gravel = 0.

Step 4: Combine the Constraints (Tool: Raster Calculator / Map Algebra)

Multiply all four reclassified grids together using the equation:

Because you are multiplying, if any single criteria fails (has a 0), the final pixel value becomes 0. The output raster will only have a value of 1 where all criteria are perfectly met.

3. The Raster Condition Numerical

Question: Generate an output raster of the suitable area for landfill site from the given raster below. The values in both rasters given below denote the distances in kilometers from settlement buildings and river network. The area located farther from the buildings and closer from the river network is the suitable for it. First of all set the conditions for both raster.

Step 1: Set the Conditions (Engineering Judgment)

The question leaves the thresholds to your judgment. We must define explicit rules:

Step 2: Reclassify Both Grids to Binary (1/0)

4. Rice Crop Suitability — Weighted Overlay (14 Marks Numerical)

This is a Weighted Overlay problem. Unlike the Boolean (AND) approach above, here each input layer has a different importance (weight), and we assign suitability scores on a 1–10 scale, then compute a weighted sum. The final value determines suitability.

The Formula

Step-by-Step Solution Method

The question provides three input grids (slope grid, soil grid, landuse grid) and a reclassification table. Follow these steps:

Step 1: Reclassify each input grid using the given Table 1 (which maps grid codes to suitability scores 1–10).

Step 2: Apply weights — multiply each reclassified cell by its weight factor (0.50, 0.30, 0.20).

Step 3: Add — sum the three weighted grids cell-by-cell to get the final suitability score per cell.

Step 4: Threshold — any cell with suitability > 6 is suitable; others are not suitable.

Step 5: Calculate area — count suitable cells × cell area (25m × 25m = 625 m² per cell).

Worked Example (Generic 3×3 Grid)

Since the actual grid data varies per exam paper, here is a complete worked example showing the method:

Suppose after reclassification, the 1–10 suitability scores for a 3×3 grid are:

Final Suitability Map & Area Calculation

Key exam tip: In the actual exam, you will be given specific grid codes and a reclassification table. Reclassify each cell first, then apply the weighted formula to every cell. Show all arithmetic clearly for full marks.

End of Chapter 7 – Spatial Analysis