13.1 Map Algebra & Raster Overlay: Local, Focal, Zonal & Global Operations

13.1 Map Algebra & Raster Overlay: Local, Focal, Zonal & Global Operations

Quick Summary: Raster spatial analysis treats continuous and discrete geographic space as regular tessellations of grid cells. Formalized by Dr. C. Dana Tomlin, Map Algebra is the foundational mathematical language for raster modeling. It categorizes all raster operations into four fundamental spatial scopes based on the geographic extent of input cells required to compute an output cell: Local (cell-by-cell at identical coordinates), Focal (neighborhood moving windows), Zonal (irregular patches sharing categorical identifiers), and Global (the entire raster extent). Understanding how these operations manage mathematical transformations, handle NoData propagation, implement conditional statements, and construct least-cost cumulative friction paths is essential for advanced spatial modeling and multi-criteria decision analysis (MCDA).

---

1. Foundations of Map Algebra & Raster Cartographic Modeling

In raster analysis, geographic phenomena are modeled as georeferenced matrices of numeric values organized into rows ($r$) and columns ($c$). In the late 1970s and 1980s, Dr. C. Dana Tomlin established Map Algebra as an algebraic syntax that treats entire raster layers as variables within algebraic equations.

   Raster A (Input 1)         Raster B (Input 2)         Output Raster (A + B)
   +---+---+---+             +---+---+---+             +---+---+---+
   | 2 | 4 | 1 |             | 1 | 0 | 3 |             | 3 | 4 | 4 |
   +---+---+---+      +      +---+---+---+      =      +---+---+---+
   | 5 | 3 | 2 |             | 2 | 1 | 1 |             | 7 | 4 | 3 |
   +---+---+---+             +---+---+---+             +---+---+---+
   | 0 | 1 | 4 |             | 3 | 2 | 0 |             | 3 | 3 | 4 |
   +---+---+---+             +---+---+---+             +---+---+---+
   (Evaluated independently at every coincident cell coordinate)

Geometric Prerequisites for Multi-Layer Map Algebra

Before executing multi-layer Map Algebra, input rasters must satisfy rigorous geometric registration constraints. If rasters do not share identical geometric configurations, the spatial processing engine must resample and interpolate cell values on the fly, which introduces positional uncertainty and radiometric artifacts:

1. Identical Coordinate Reference System (CRS): Rasters must share the same projection and datum. On-the-fly reprojection during raster overlay can distort pixel boundaries and alter cell area.
2. Identical Cell Resolution: All participating rasters should share the same cell dimensions (e.g., $10,\text{m} \times 10,\text{m}$). If resolutions diverge, the software defaults to the coarsest cell resolution or the resolution specified in the geoprocessing environment settings.
3. Coincident Grid Registration (Snapping): The origin coordinates and bounding grid lines must match precisely. If cell corners are offset by even a fraction of a pixel, the software must resample one grid onto the other, blurring discrete categorical boundaries.
4. Consistent Spatial Extent: The spatial extent defines the minimum bounding rectangle enclosing the computation. By default, most engines set the output extent to the Intersection (common bounding box) or Union of all input layers.

---

2. Tomlin's Operational Taxonomy

Tomlin classified Map Algebra operations into four spatial scopes based on the spatial relationship between the input cells used for the calculation and the resulting output cell:

   LOCAL SCOPE                FOCAL SCOPE               ZONAL SCOPE               GLOBAL SCOPE
   +---+---+---+             +---+---+---+             +---+---+---+             +---+---+---+
   |   |   |   |             | x | x | x |             | Z1| Z1| Z2|             | x | x | x |
   +---+---+---+             +---+---+---+             +---+---+---+             +---+---+---+
   |   | X |   |             | x | X | x |             | Z1| Z1| Z2|             | x | X | x |
   +---+---+---+             +---+---+---+             +---+---+---+             +---+---+---+
   |   |   |   |             | x | x | x |             | Z2| Z2| Z2|             | x | x | x |
   +---+---+---+             +---+---+---+             +---+---+---+             +---+---+---+
   Input: Same Cell (x,y)    Input: Neighborhood       Input: Common Zone        Input: ALL Cells

Overview of Operational Scopes

---

3. Local Operations: Mathematical Functions, Reclassification & Conditional Logic

Local operations compute an output value at location $(x,y)$ solely from the value(s) residing at $(x,y)$ in the input raster(s). They do not interrogate neighboring cells.

Mathematical and Trigonometric Operations

• Unary Operations (Single Layer): Evaluates a mathematical function on individual cells: $\text{Output} = f(\text{Input})$. Examples include $\log(x)$, $\exp(x)$, $\text{abs}(x)$, $\sqrt{x}$, and trigonometric functions ($\sin$, $\cos$, $\tan$) commonly used in solar radiation and slope aspect decomposition.
• Binary / Multi-Layer Operations: Combines multiple overlapping rasters: $\text{Output} = \text{Raster}_1 + \text{Raster}_2$ or $\text{Output} = \frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red}}$ (Normalized Difference Vegetation Index).

Reclassification

Reclassification reassigns cell values to new output values using lookup tables or mathematical intervals. It is used to:

• Group continuous values into discrete classes (e.g., converting continuous slope degrees into 5 slope-hazard categories).
• Standardize disparate thematic scales into a uniform suitability ranking (e.g., reclassifying land cover classes into scores from 1 to 10).
• Replace unwanted or erroneous values with NoData.

Conditional Evaluation: The Con Statement

Conditional modeling performs branch logic on a cell-by-cell basis, mirroring the standard IF-THEN-ELSE programming construct:

$\text{Output} = \text{Con}(\text{Conditional\_Raster},\; \text{True\_Expression},\; [\text{False\_Expression}])$

If the condition evaluates to TRUE (non-zero numeric value), the cell is assigned the True_Expression. If FALSE (zero), it receives the False_Expression (or NoData if omitted).

# Example Map Algebra conditional syntax in Python (ArcPy / Rasterio / NumPy)
# Assign value 10 to steep slopes (>25%), otherwise retain original elevation
OutRaster = Con(Slope_Grid > 25, 10, Elev_Grid)

# Nested Conditional: Multi-tiered ecological classification
OutHabitat = Con(Landuse == 4, 1, Con((Landuse == 5) & (Slope < 10), 2, 0))

NoData Propagation and Handling

In raster architectures, NoData represents missing, uncollected, or masked spatial data. In standard local arithmetic operations, NoData is infectious:

$\text{NoData} + 15 = \text{NoData}$ $\text{NoData} \times 0 = \text{NoData}$

If an analyst adds two elevation rasters and a single cell in one raster is NoData, the output cell becomes NoData. To prevent unintended data loss, analysts use conditional statements paired with IsNull() testing functions:

# Replace NoData cells with a baseline value of 0 prior to addition
Clean_Raster1 = Con(IsNull(Raster1), 0, Raster1)
Final_Sum = Clean_Raster1 + Raster2

---

4. Focal Operations: Neighborhood Moving Windows & Filter Kernels

Focal operations calculate the output value at location $(x,y)$ as a function of the input cells falling within a specified neighborhood (or kernel) surrounding $(x,y)$. The neighborhood shifts sequentially across every cell in the grid.

   Input 3x3 Neighborhood               Kernel Weights              Output Focal Mean Cell
   +----+----+----+                     +----+----+----+            
   | 10 | 12 | 14 |                     | 1  | 1  | 1  |            Sum = (10+12+14+12+20+16+11+13+15)
   +----+----+----+                     +----+----+----+                = 123
   | 12 | 20 | 16 |          *          | 1  | 1  | 1  |    ===>    Count = 9
   +----+----+----+                     +----+----+----+            Output Center Cell = 123 / 9
   | 11 | 13 | 15 |                     | 1  | 1  | 1  |                               = 13.67
   +----+----+----+                     +----+----+----+            (Replaces original value 20)

Neighborhood Morphologies

• Rectangle: Specified by dimensions (e.g., $3 \times 3$, $5 \times 5$, $7 \times 7$ cells). Must have odd integer dimensions to preserve an unambiguous center cell.
• Circle: Defined by a radial distance from the center cell. Cells whose centroids fall within the radius are included.
• Annulus (Doughnut): Defined by an inner and outer radius, excluding immediate adjacent cells while capturing surrounding contextual rings.
• Wedge: Defined by a directional angle and radius, oriented along wind or hydrological flow vectors.

Focal Statistical Functions and Filter Types

Boundary and Margin Effects in Focal Processing

When a moving window reaches the edge of a raster, the kernel extends beyond the valid data boundary. Raster engines handle boundary conditions using three common methods:

1. Truncated Kernel (Default): Calculates the statistic using only valid cells falling inside the raster extent (e.g., a $3 \times 3$ corner cell calculates the mean across 4 available cells instead of 9).
2. NoData Padding: If any cell in the kernel falls outside the raster, the output cell is assigned NoData, shrinking the output raster boundary by $(k - 1) / 2$ cells.
3. Mirror / Constant Padding: Extrapolates edge cells or reflects adjacent values into the outer buffer.

---

5. Zonal Operations: Aggregating Across Discrete Spatial Units

Zonal operations calculate output statistics for a Value Raster based on spatial zones delineated by a Zone Layer (which may be an integer raster or a vector polygon layer).

   Zone Layer (Integer Raster)         Value Layer (Continuous DEM)       Zonal Mean Calculation
   +----+----+----+----+               +----+----+----+----+              
   | 1  | 1  | 2  | 2  |               | 10 | 20 | 50 | 60 |              Zone 1 Cells: {10, 20, 15, 25}
   +----+----+----+----+               +----+----+----+----+              Zone 1 Mean = 70 / 4 = 17.5
   | 1  | 1  | 2  | 2  |               | 15 | 25 | 70 | 80 |              
   +----+----+----+----+               +----+----+----+----+              Zone 2 Cells: {50, 60, 70, 80}
   | 3  | 3  | 3  | 3  |               | 10 | 12 | 14 | 16 |              Zone 2 Mean = 260 / 4 = 65.0
   +----+----+----+----+               +----+----+----+----+              
   | 3  | 3  | 3  | 3  |               | 18 | 20 | 22 | 24 |              Zone 3 Cells: {10,12,14,16,18,20,22,24}
   +----+----+----+----+               +----+----+----+----+              Zone 3 Mean = 136 / 8 = 17.0

Core Characteristics of Zones

• Definition of a Zone: A zone comprises all cells possessing the exact same zone identifier, regardless of whether those cells are physically contiguous. For example, all fragmented patches coded Zone = 4 (e.g., Deciduous Forest) constitute a single zone.
• Zonal Statistics: Summarizes the value raster per zone, producing either an output summary table (one row per zone ID) or a new raster where every cell in a zone is stamped with the summary statistic (e.g., mean, maximum, minimum, sum, range, standard deviation, variety, majority, or minority).
• Zonal Geometry: Measures the spatial and geometric properties of each zone itself (without requiring an external value raster), such as total area, perimeter, centroid coordinates, and circularity/thickness.

[!IMPORTANT] Focal vs. Zonal Distinction:

• In a Focal operation, the neighborhood shifts with every individual cell; each cell has its own unique, overlapping neighborhood.
• In a Zonal operation, zones are fixed geographic parcels or categories. Every cell belonging to the same zone receives the identical aggregated statistical value.

---

6. Global Operations: Distance Surfaces & Least-Cost Paths

Global operations compute output cell values as a function of the entire raster grid. Distance modeling and least-cost routing are the most widely used global operations.

Euclidean Distance & Allocation

• Euclidean Distance: Calculates the straight-line distance from the center of each cell to the nearest source feature (points, lines, or polygons): $d = \sqrt{(x_2 - x_1)^2 + (y_2 - y_1)^2}$. Produces a continuous floating-point surface.
• Euclidean Direction: Calculates the azimuth angle (in degrees from 0 to 360, measuring clockwise from North) to the nearest source.
• Euclidean Allocation: Assigns each cell the unique identity or attribute of the closest source feature, generating a raster equivalent of Voronoi/Thiessen polygons.

Cost Distance & Least-Cost Path Modeling

When travel across a landscape is impeded by terrain, vegetation, water bodies, or infrastructure, straight-line Euclidean distance fails to reflect real-world travel effort. Cost Distance modeling calculates the minimum cumulative cost of traveling from any cell back to the nearest source across a continuous Cost Raster (friction surface).

   +------------------+       +------------------+       +--------------------------+
   | Source Features  |       | Cost/Friction    |       | Cost Distance Surface    |
   | (Points / Lines) |   +   | Surface          |  ==>  | (Cumulative minimum cost |
   |                  |       | (Slope, Landuse) |       |  from each cell to source|
   +------------------+       +------------------+       +--------------------------+
                                                                      |
                                                                      v
                                                         +--------------------------+
                                                         | Cost Backlink Direction  |
                                                         | Raster (Codes 1 to 8)    |
                                                         +--------------------------+
                                                                      |
                                                                      v
                                                         +--------------------------+
                                                         | Least-Cost Path Tool     |
                                                         | (Traces optimal path     |
                                                         |  from target to source)  |
                                                         +--------------------------+

Cost Distance Mechanics

The algorithm spreads outward from source cells using Dijkstra-like graph searching across the cell lattice. The cost to travel between two adjacent cells ($A$ and $B$) is the average of their cost values multiplied by the distance:

• Orthogonal Movement: $\text{Cost} = \frac{\text{Cost}_A + \text{Cost}_B}{2} \times \text{Cell_Size}$
• Diagonal Movement: $\text{Cost} = \frac{\text{Cost}_A + \text{Cost}_B}{2} \times \sqrt{2} \times \text{Cell_Size} \approx \frac{\text{Cost}_A + \text{Cost}_B}{2} \times 1.4142 \times \text{Cell_Size}$

The Cost Backlink (Direction) Raster

To compute a Least-Cost Path, the system requires two outputs from the Cost Distance tool:

1. The Accumulated Cost Surface (floating-point cumulative cost).
2. The Cost Backlink Raster (integer grid coding direction from 1 to 8, identifying the exact neighboring cell to traverse to return to the source along the cheapest path).

The Least-Cost Path tool begins at the designated destination cell and follows the backlink directions downhill across the cumulative cost surface until it reaches the source, guaranteeing the mathematical minimum-cost route.

---

7. Multi-Criteria Decision Analysis (MCDA) & Weighted Overlay

Multi-Criteria Decision Analysis (MCDA) combines multiple spatial factors to evaluate land suitability or vulnerability. While Boolean overlay evaluates strict pass/fail criteria using binary logic (AND/OR), Weighted Overlay accommodates trade-offs by evaluating continuous degrees of suitability.

   CRITERION 1: Slope            CRITERION 2: Land Cover        CRITERION 3: Dist. to Road
   Reclassified (1 to 9)         Reclassified (1 to 9)          Reclassified (1 to 9)
   Weight: 40%                   Weight: 35%                    Weight: 25%
        \                              |                             /
         \                             |                            /
          +----------------------------+---------------------------+
                                       |
                                       v
               WEIGHTED OVERLAY = (Slope * 0.40) + (LandCover * 0.35) 
                                 + (Distance * 0.25)
                                       |
                                       v
               FINAL SUITABILITY RASTER (Continuous Scores 1 to 9)

Standard MCDA Workflow Steps

1. Define Objective & Criteria: Identify environmental, legal, and engineering requirements (e.g., siting a solar energy facility).
2. Derive Continuous Factor Surfaces: Calculate slope from a DEM, Euclidean distance from transmission lines, and solar radiation surfaces.
3. Reclassify to a Standardized Scale: Convert incompatible physical units (degrees, meters, dollars) into a common suitability scale (e.g., 1 to 5, 1 to 9, or 0 to 100), where higher values indicate greater suitability.
4. Assign Criterion Weights: Allocate percentage weights to criteria reflecting their relative importance. All weights must sum to exactly 100% (or 1.0).
5. Apply Weighted Sum / Overlay: Multiply each reclassified raster by its assigned weight and sum the results: $\text{Suitability} = \sum_{i=1}^{n} w_i \times S_i$
6. Apply Absolute Constraints (Restricted Masking): Multiply the final suitability surface by a binary constraint raster ($1 = \text{Permitted}, 0 = \text{Prohibited}$) to exclude legally or physically unsuitable areas (e.g., wetlands, public parks, water bodies).

---

8. Practical Geospatial Scenario: Emergency Wildfire Evacuation Corridor

Analytical Problem

A regional emergency management authority must identify an optimal, secondary overland evacuation corridor to connect an isolated mountain community to an interstate highway during active wildfire conditions.

Workflow Implementation

   Step 1: Build Friction/Cost Surface (Local Operations)
           - Slope Raster derived from DEM (steep slopes increase travel time exponentially)
           - Vegetation Fuel Load (dense brush = high transit risk)
           - Road Class Grid (paved = low cost, dirt trail = moderate cost, off-road = high cost)
           Cost_Surface = (Slope_Cost * 0.40) + (Fuel_Cost * 0.35) + (Road_Cost * 0.25)
           
   Step 2: Mask Absolute Hazards (Local Con / IsNull)
           Cost_Surface_Masked = Con(Active_Fire_Buffer == 1, NoData, Cost_Surface)
           
   Step 3: Calculate Cumulative Cost & Backlink (Global Operation)
           Accum_Cost, Backlink = CostDistance(Source=Highway, CostRaster=Cost_Surface_Masked)
           
   Step 4: Delineate Optimal Corridor (Least-Cost Path)
           Evacuation_Path = CostPath(Destination=Community, CostDistance=Accum_Cost, 
                                      Backlink=Backlink, PathType="EACH_CELL")

By executing this workflow, the agency identifies the mathematically optimal path that minimizes travel time and hazard exposure, avoiding impassable terrain while steering clear of active fire fronts.

---

9. Common Exam Traps & Pitfalls

[!CAUTION] Exam Trap 13.1.1: Conflating Focal and Zonal Neighborhoods. Questions frequently ask which operation computes statistics across areas sharing an attribute versus areas defined by a spatial filter. Remember: Focal operations move cell-by-cell across the grid, evaluating a sliding geometric window centered over each target cell. Zonal operations aggregate across predefined thematic parcels or classes sharing the same integer ID; zones are static across the dataset and do not require cells to be spatially contiguous.

[!CAUTION] Exam Trap 13.1.2: Infectious NoData in Local Arithmetic. If an exam question asks what happens when a raster containing valid numbers is multiplied or added to a raster containing NoData cells, the answer is that the result is NoData at those locations. The only way to preserve underlying values is to use a conditional statement (e.g., Con(IsNull(A), B, A + B)).

[!CAUTION] Exam Trap 13.1.3: Least-Cost Path Input Requirements. A common question asks: "What inputs are required to compute a least-cost path from a destination to a source?" The incorrect distractor is "The Source Layer and the Friction Surface." The correct answer requires the Accumulated Cost Surface and the Cost Backlink (Direction) Raster generated by the Cost Distance tool. You cannot trace a least-cost path directly from a raw friction surface alone.