Geographic Information Systems and Remote Sensing
Principles, Models, and Analytical Methods for Advanced Spatial Science
Preface
This text is written for advanced undergraduate and graduate students, and for practitioners who wish to place their spatial work on firmer theoretical ground. It treats Geographic Information Systems (GIS) and remote sensing (RS) as a single, integrated discipline: remote sensing supplies calibrated, repeatable measurements of the Earth’s surface, and GIS supplies the data models, reference frameworks, and analytical operators that turn those measurements into information, inference, and decisions.
The emphasis throughout is on understanding why a method works. Each chapter therefore develops the physical, geometric, or statistical principle first, states it mathematically, and then examines its assumptions, limits, and typical failure modes. The book intentionally contains no program code; the ideas are software-independent and can be implemented in any modern geospatial environment. Instructors may pair each chapter with laboratory sessions in the software of their choice.
How the book is organised
- Chapters 1–3 establish the foundations of geospatial science: what geographic information is, how position is defined on a curved Earth, and how spatial phenomena are represented in digital form.
- Chapters 4–6 develop the physics of remote sensing, the design of sensors and platforms, and the preprocessing needed to convert raw digital numbers into physically meaningful, geometrically correct data.
- Chapters 7–9 cover image enhancement, classification and accuracy assessment, and the specialised domains of radar, LiDAR, and imaging spectroscopy.
- Chapters 10–12 address spatial analysis and geostatistics, time-series and change detection, and the questions of uncertainty, ethics, and reproducibility that separate rigorous work from map-making.
Prerequisites. Introductory calculus, linear algebra (vectors, matrices, eigen-decomposition), and elementary statistics (distributions, regression, hypothesis testing).
Notation
| Symbol | Meaning |
|---|---|
| \(\lambda\) | wavelength |
| \(\nu\) (nu), \(f\) | frequency (Hz); in Chapter 2, \(\nu\) also denotes the prime-vertical radius of curvature |
| \(\varphi\), \(\ell\) | geodetic latitude and longitude (\(\ell\) avoids a clash with wavelength) |
| \(h\), \(H\), \(N\) | ellipsoidal height, orthometric height, geoid undulation (Chapter 2) |
| \(L\) | radiance (\(\mathrm{W\,m^{-2}\,sr^{-1}\,\mu m^{-1}}\) when spectral) |
| \(E\) | irradiance (\(\mathrm{W\,m^{-2}}\)) |
| \(\rho\) | reflectance (unitless, 0–1) |
| \(\tau\) | atmospheric transmittance |
| \(\theta_s\), \(\theta_v\) | solar zenith angle, view zenith angle |
| \(\mathbf{x}\) | pixel vector of spectral values, \(\mathbf{x}=(x_1,\dots,x_B)^{\top}\) |
| \(\boldsymbol{\mu}\), \(\boldsymbol{\Sigma}\) | mean vector and covariance matrix |
| \(\gamma(h)\) | semivariance at lag \(h\) |
| \(\mathbf{W}\), \(w_{ij}\) | spatial weights matrix and its elements |
1 Foundations of Geospatial Science
1.1 What is geographic information?
Geographic information describes where something is, what it is like there, and when it was so. Every geographic datum therefore has, at minimum, a spatial component (position, extent, shape), a thematic component (attributes), and a temporal component (validity in time). Collectively these are sometimes called the where–what–when triple of a spatial observation.
Definition 1.1 (Geographic phenomenon) A geographic phenomenon is something that can be observed at, or is associated with, locations on or near the Earth’s surface, and that varies in a way that can be described by measurement or classification.
Geographic phenomena are conceptualised in two complementary ways:
- Discrete objects have well-defined boundaries and identity: a building, a river channel, a health facility, an administrative district. They are naturally represented by points, lines, and polygons.
- Continuous fields have a value at every location: temperature, elevation, soil moisture, vegetation greenness. They are naturally represented by rasters, triangulated surfaces, or functions.
The distinction is one of conceptual model, not of physical reality. A forest may be treated as an object (a polygon with attributes) or as a field (a canopy-cover surface). The choice determines which operations are meaningful and which errors are likely.
1.1.1 Tobler’s first law and spatial dependence
The organising principle of spatial analysis is Tobler’s first law of geography: everything is related to everything else, but near things are more related than distant things. Statistically, this is spatial autocorrelation: observations close together tend to be more similar (positive autocorrelation) than observations far apart. Autocorrelation is a blessing when we wish to interpolate or classify, and a curse when we wish to apply methods that assume independent observations. It appears repeatedly throughout this book (Chapters 8, 10, and 12).
A second property, spatial heterogeneity (non-stationarity), means that relationships themselves vary over space: the association between rainfall and vegetation greenness in a humid highland need not equal that in a semi-arid plain.
1.2 Components and functions of a GIS
A GIS is commonly described as an organised collection of five components:
- Data — spatial and attribute data, together with metadata.
- Software — tools for capture, storage, analysis, and display.
- Hardware and infrastructure — from field devices to cloud platforms.
- Methods — the analytical procedures and conceptual models applied.
- People — the analysts, developers, and decision-makers whose judgement completes the system.
Functionally, a GIS supports a cycle of capture (surveying, GNSS, digitising, remote sensing), management (databases, indexing, versioning), analysis (queries, overlays, modelling, statistics), and communication (cartography, dashboards, web services). Geographic information science (GIScience) is the scholarly field that studies the concepts and methods underlying these functions, including representation, uncertainty, and scale.
1.3 Remote sensing within the geospatial workflow
Definition 1.2 (Remote sensing) Remote sensing is the science and practice of obtaining information about an object, area, or phenomenon through the analysis of data acquired by a device that is not in contact with it, typically by measuring electromagnetic radiation reflected, emitted, or backscattered by the surface.
Remote sensing systems are passive when they measure radiation from a natural source (reflected sunlight, or thermal emission from the Earth) and active when they supply their own illumination (radar, LiDAR). The general chain from observation to decision can be summarised as:
\[ \text{energy source} \rightarrow \text{atmosphere} \rightarrow \text{target} \rightarrow \text{atmosphere} \rightarrow \text{sensor} \rightarrow \text{processing} \rightarrow \text{information} \rightarrow \text{decision}. \tag{1.1}\]
Each link in Equation 1.1 modifies the signal and therefore contributes to the uncertainty in the final product. Much of Chapters 4 to 6 is devoted to understanding and removing those modifications.
GIS and remote sensing are complementary. Remote sensing produces rasters rich in radiometric information but poor in semantics: a pixel value is not a land-cover class. GIS provides the reference frame (coordinate systems), the ancillary layers (elevation, soils, administrative boundaries, field data) and the analytical framework needed to interpret, validate, and apply the imagery.
1.4 Scale, resolution, extent, and support
Four related concepts are frequently confused. Cartographic scale is the ratio of a distance on a map to the corresponding distance on the ground, expressed as a representative fraction (RF):
\[ \mathrm{RF} = \frac{d_{\text{map}}}{d_{\text{ground}}}, \qquad \text{e.g. } \mathrm{RF}=\frac{1}{50\,000}. \tag{1.2}\]
Large-scale maps (large RF, such as 1:5,000) show small areas in detail; small-scale maps (1:1,000,000) show large areas generalised. The other three concepts concern data rather than maps:
- Extent is the total area or time span covered by a dataset.
- Resolution (or grain) is the smallest unit distinguished: the pixel size of an image, the minimum mapping unit of a polygon layer.
- Support is the size, shape, and orientation of the spatial unit on which a value is defined (a point measurement, a \(30\times30\) m pixel, a district). Values defined on different supports cannot be compared or combined without accounting for the change of support.
Statistical results computed from areal aggregates depend on the scale (how many zones) and the zoning (how the boundaries are drawn) of the aggregation. Correlations, regression coefficients, and cluster patterns can change, sometimes reversing sign, when the same underlying data are aggregated differently. MAUP is not an error to be eliminated but a property to be acknowledged and tested (see Section 10.7).
1.5 Application domains
Integrated GIS and remote sensing underpin work in: agriculture and food security (crop condition, yield forecasting, irrigation mapping); hydrology and water resources (surface-water dynamics, snow cover, soil moisture, flood mapping); forestry and carbon accounting (biomass, deforestation, degradation); urban planning (impervious surface, land-use change, heat islands); disaster management (burned area, landslides, earthquake deformation); public and environmental health (environmental covariates of vector-borne disease, exposure surfaces, health-facility accessibility); and climate science (albedo, land-surface temperature, sea-ice extent).
1.6 A brief history
Systematic spatial analysis has old roots: the 1854 cholera map of London that linked deaths to a water pump is a classic early example of thematic mapping used for inference. Aerial photography matured in the two World Wars; the launch of the first Earth Resources Technology Satellite (later Landsat 1) in 1972 began the continuous civil satellite record. The Canada Geographic Information System of the 1960s pioneered computer-based spatial analysis, while the availability of civil GPS, the opening of the Landsat archive (2008), the Copernicus programme’s free and open Sentinel data, and cloud-based computing have since made global-scale analysis routine.
Exercises
Exercise 1.1 For each of the following, state whether the most natural conceptual model is an object or a field, and justify your choice in one sentence. Then describe a situation in which the other model would be preferable: (a) a lake; (b) sea-surface temperature; (c) a road network; (d) population density; (e) a household.
Exercise 1.2 A topographic map has a representative fraction of 1:50,000. (a) What ground distance corresponds to 3.2 cm on the map? (b) A field measures 4 cm × 2.5 cm on this map; what is its ground area in hectares? (c) On a 1:250,000 map, what map distance represents the same 1.6 km?
Exercise 1.3 Explain, with a small numerical thought-experiment involving four spatial units, how the correlation between two variables can differ when the units are merged in two different ways. Why does this matter for studies that relate disease incidence to environmental covariates aggregated by administrative district?
Exercise 1.4 Design problem. A ministry wishes to monitor the extent of irrigated agriculture in a river basin every season. List the spatial, thematic, and temporal requirements of the product, and state which components of a GIS and which remote sensing data characteristics (Chapter 5) you would need to satisfy them.
2 Geodesy, Coordinate Systems, and Map Projections
Every spatial analysis presupposes an answer to the question: where, exactly? Because the Earth is neither flat nor a perfect sphere, position must be defined relative to a mathematical model of its shape and a coordinate framework tied to that model.
2.1 The shape of the Earth
2.1.1 Sphere, ellipsoid, and geoid
Three surfaces are relevant, as illustrated in Figure 2.1:
- The physical surface of the land and ocean, on which we measure.
- The geoid, an equipotential surface of the Earth’s gravity field that coincides, in the oceans, with mean sea level extended under the continents. It is irregular because the Earth’s mass is unevenly distributed.
- The reference ellipsoid, a smooth, mathematically defined oblate spheroid chosen to approximate the geoid, on which geodetic coordinates are computed.
An ellipsoid of revolution is defined by its semi-major axis \(a\) and semi-minor axis \(b\), or equivalently by \(a\) and the flattening \(f\):
\[ f=\frac{a-b}{a}, \qquad e^{2}=\frac{a^{2}-b^{2}}{a^{2}}=2f-f^{2}, \tag{2.1}\]
where \(e\) is the first eccentricity. The WGS 84 ellipsoid, used by GPS and most global datasets, has \(a=6\,378\,137\) m and \(1/f=298.257223563\), so that \(b\approx 6\,356\,752\) m: the Earth’s equatorial radius exceeds its polar radius by about 21 km.
The three heights are linked by
\[ h = H + N, \tag{2.2}\]
where \(h\) is the ellipsoidal height (what GNSS receivers measure), \(H\) the orthometric height (elevation above mean sea level, what most maps and engineers use), and \(N\) the geoid undulation (separation of the geoid from the ellipsoid, ranging globally from about \(-106\) m to \(+85\) m). Mixing \(h\) and \(H\) is a common source of vertical errors of tens of metres.
2.1.2 Radii of curvature
Distances and areas on the ellipsoid depend on latitude because the curvature varies. The meridian radius of curvature \(\rho_m\) and the prime-vertical radius of curvature \(\nu\) at geodetic latitude \(\varphi\) are
\[ \rho_m=\frac{a\,(1-e^{2})}{\left(1-e^{2}\sin^{2}\varphi\right)^{3/2}},\qquad \nu=\frac{a}{\sqrt{1-e^{2}\sin^{2}\varphi}}. \tag{2.3}\]
Consequently, one degree of latitude is about \(\rho_m\,\pi/180\) (≈ 110.6 km at the equator, ≈ 111.7 km at the poles), while one degree of longitude spans \(\nu\cos\varphi\,\pi/180\) and shrinks to zero at the poles.
2.2 Geodetic datums
Definition 2.1 (Geodetic datum) A geodetic datum specifies the reference ellipsoid and the position and orientation of that ellipsoid relative to the Earth (and, for a modern datum, the realisation of the reference frame by a network of control points).
Local datums are fitted to a region so that the ellipsoid closely follows the local geoid; they give good agreement with local surveys but are offset by tens to hundreds of metres from a global frame. Global (geocentric) datums such as WGS 84 and the International Terrestrial Reference Frame (ITRF) place the ellipsoid’s centre at the Earth’s centre of mass.
Geodetic coordinates \((\varphi,\ell,h)\) can be converted to Earth-centred, Earth-fixed (ECEF) Cartesian coordinates by
\[ \begin{aligned} X&=(\nu+h)\cos\varphi\cos\ell,\\ Y&=(\nu+h)\cos\varphi\sin\ell,\\ Z&=\left[\nu\,(1-e^{2})+h\right]\sin\varphi. \end{aligned} \tag{2.4}\]
Converting between datums is usually done in this Cartesian space with a seven-parameter Helmert (similarity) transformation:
\[ \begin{pmatrix}X'\\Y'\\Z'\end{pmatrix} = \begin{pmatrix}T_X\\T_Y\\T_Z\end{pmatrix} +(1+s)\, \begin{pmatrix} 1 & -r_Z & r_Y\\ r_Z & 1 & -r_X\\ -r_Y & r_X & 1 \end{pmatrix} \begin{pmatrix}X\\Y\\Z\end{pmatrix}, \tag{2.5}\]
with three translations \((T_X,T_Y,T_Z)\), three small rotations \((r_X,r_Y,r_Z)\) in radians, and one scale factor \(s\). When the datum shift is a grid-based correction (as in several national frameworks), interpolation in a correction grid replaces Equation 2.5.
2.3 Map projections
A map projection is a systematic transformation of positions on the curved surface to positions on a plane: \[ (x,y)=\bigl(f(\varphi,\ell),\;g(\varphi,\ell)\bigr). \] No projection can preserve all geometric properties: a sphere is not a developable surface (Gauss’s Theorema Egregium), so every flat map distorts area, shape, distance, or direction.
2.3.1 Classification
Projections are classified by:
- the developable surface used: cylindrical, conical, or azimuthal (planar);
- the aspect: normal, transverse, or oblique;
- the contact: tangent or secant to the reference surface;
- the property preserved: conformal (local shape and angles), equal-area (equivalent), equidistant (distances from one or two points, or along meridians), or compromise.
2.3.2 Distortion and scale factors
Let \(h_s\) and \(k_s\) be the local linear scale factors along the meridian and the parallel respectively (ratio of projected to true length). Then:
\[ \text{conformal: } h_s=k_s;\qquad \text{equal-area: } h_s\,k_s=1 . \tag{2.6}\]
Tissot’s indicatrix, the image of an infinitesimal circle on the sphere, visualises distortion: a circle of unchanged shape indicates conformality, an ellipse of unchanged area indicates equivalence.
2.3.3 Selected projections
For a sphere of radius \(R\) and central meridian \(\ell_0\):
Mercator (conformal cylindrical; straight rhumb lines): \[ x=R(\ell-\ell_0),\qquad y=R\ln\tan\!\left(\frac{\pi}{4}+\frac{\varphi}{2}\right),\qquad k_s=h_s=\sec\varphi . \tag{2.7}\]
Because the scale grows as \(\sec\varphi\), areas are inflated by \(\sec^{2}\varphi\): at \(60^{\circ}\) latitude linear scale doubles and area scale quadruples. The Web Mercator variant used by web maps is therefore unsuitable for measuring areas or distances.
Lambert cylindrical equal-area: \[ x=R(\ell-\ell_0),\qquad y=R\sin\varphi . \tag{2.8}\]
Transverse Mercator and UTM. The Universal Transverse Mercator system divides the world into 60 zones each \(6^{\circ}\) wide and applies a transverse Mercator projection with a central-meridian scale factor \(k_0=0.9996\) (a secant-like reduction that limits scale error to about \(\pm 0.04\%\) within a zone), a false easting of 500 000 m, and, for the southern hemisphere, a false northing of 10 000 000 m. The zone number and central meridian for longitude \(\ell\) (degrees east, from \(-180\) to \(+180\)) are
\[ Z=\left\lfloor\frac{\ell+180}{6}\right\rfloor+1,\qquad \ell_0=6Z-183 . \tag{2.9}\]
Conic projections (Lambert conformal conic, Albers equal-area conic) suit mid-latitude regions extended in the east–west direction; azimuthal projections suit polar regions or circular extents.
| Analytical need | Prefer | Reason |
|---|---|---|
| Area statistics, land-cover totals | Equal-area (Albers, Lambert CEA, Mollweide globally) | Preserves area |
| Navigation, large-scale mapping, angles | Conformal (UTM, Lambert conformal conic) | Preserves local shape and bearing |
| Distances from a facility | Azimuthal equidistant centred on the facility | True distance from the centre |
| Regional analysis of an N–S strip | Transverse Mercator / UTM zone | Low distortion in a narrow zone |
| Web display only | Web Mercator | Tile compatibility, not for measurement |
2.4 Distances on the sphere and ellipsoid
The great-circle distance between points \((\varphi_1,\ell_1)\) and \((\varphi_2,\ell_2)\) on a sphere is computed stably with the haversine formula:
\[ a=\sin^{2}\!\frac{\Delta\varphi}{2}+\cos\varphi_1\cos\varphi_2\,\sin^{2}\!\frac{\Delta\ell}{2},\qquad d=2R\arcsin\!\sqrt{a}. \tag{2.10}\]
The spherical error can reach about 0.5 %. Geodesic distances on the ellipsoid (shortest paths) require iterative algorithms of Vincenty or Karney; modern libraries implement Karney’s method, which is accurate to nanometres.
2.5 Coordinate reference systems in practice
A coordinate reference system (CRS) combines a datum with a coordinate system (geographic or projected, with defined axes and units). Well-formed CRS definitions are identified by authority codes (for example, EPSG:4326 for geographic WGS 84, or the EPSG codes for individual UTM zones) and are exchanged in well-known-text or PROJ notations. Good practice:
- Store data in a clearly defined CRS and never rely on assumptions.
- Perform measurement and overlay in an appropriate projected CRS (equal-area for area, equidistant or local UTM for distances).
- Remember that reprojecting a raster involves resampling, which alters cell values (Chapter 6); reproject vectors freely but rasters sparingly.
- Beware of axis-order conventions (latitude–longitude versus longitude–latitude) and of silent datum mismatches that produce systematic offsets.
Exercises
Exercise 2.1 Determine the UTM zone number and central meridian for (a) longitude \(36.82^{\circ}\)E and (b) longitude \(122.4^{\circ}\)W. Explain why a study area that straddles a zone boundary poses a problem for analysis in UTM coordinates, and suggest two ways to handle it.
Exercise 2.2 On a sphere of radius \(R=6371\) km, compute (a) the length of one degree of latitude, and (b) the length of one degree of longitude at latitude \(60^{\circ}\). Comment on the implication for a grid defined in degrees.
Exercise 2.3 Compute the linear and areal scale factors of the Mercator projection at latitude \(60^{\circ}\). A lake has a true area of 250 km² at that latitude; what area would it appear to have if measured on a Mercator map?
Exercise 2.4 A GNSS receiver at point \(P\) records an ellipsoidal height \(h=1683.2\) m. A geoid model gives \(N=-18.5\) m at \(P\). Compute the orthometric height \(H\). What would be the error if \(h\) were used in place of \(H\) in a flood-inundation model?
Exercise 2.5 Using \(R=6371.0\) km and Equation 2.10, compute the great-circle distance between \(A=(1.2921^{\circ}\text{S},\,36.8219^{\circ}\text{E})\) and \(B=(0.5143^{\circ}\text{N},\,35.2698^{\circ}\text{E})\). Would you expect the ellipsoidal geodesic distance to be larger or smaller? By roughly how much?
Exercise 2.6 Design problem. A team compiles national land-cover statistics from a mosaic of imagery delivered in several UTM zones, and reports the area of each class. Describe the sequence of CRS operations you would use, the projection you would choose, and how you would demonstrate that the areas are correct.
3 Spatial Data Models and Data Structures
3.1 From reality to representation
A spatial database is a model of reality built through a chain of abstraction: from the real world, to a conceptual model (objects and fields), to a logical model (vector, raster, network, TIN), and finally to a physical implementation (file formats, database tables). At each step, choices of resolution, generalisation, and classification introduce uncertainty.
3.1.1 Levels of measurement
Attribute values are measured on one of four scales (Stevens’ typology), which determines the valid operations:
- Nominal: categories without order (land-cover class). Only counts and modes are meaningful.
- Ordinal: ordered categories (soil-erosion risk: low, medium, high). Medians and ranks are meaningful, differences are not.
- Interval: equal intervals without a natural zero (temperature in °C). Differences and means are meaningful.
- Ratio: equal intervals with an absolute zero (rainfall, elevation above a datum, population). Ratios are meaningful.
The two principal digital representations are contrasted in Figure 3.1. The same three features (a forest polygon, a river line, and two wells) are stored either as coordinate geometries with attributes (vector) or as values in a regular grid of cells (raster). The raster version is simpler to compute with but is limited by its cell size: the river is a chain of cells, and features smaller than a cell can vanish or be displaced.
3.2 The vector data model
Vector data represent geographic objects with points (0-dimensional), polylines (1-dimensional) and polygons (2-dimensional), each linked to an attribute record. The Open Geospatial Consortium’s Simple Features model standardises these geometries, including multi-part variants and collections.
3.2.1 Topology
Topology describes properties that are invariant under continuous deformation: connectivity (which lines join at which nodes), adjacency (which polygons share a border), and containment (which features lie inside which). Topological data structures reduce redundancy (shared boundaries are stored once), allow validation (no gaps, no overlaps, no dangling arcs), and support network analysis. Many modern systems store non-topological geometries and enforce topological rules on demand.
3.2.2 Geometric measures
The area of a simple polygon with vertices \((x_i,y_i)\), \(i=1,\dots,n\), listed in order and with \((x_{n+1},y_{n+1})=(x_1,y_1)\), follows the shoelace formula:
\[ A=\frac{1}{2}\left|\sum_{i=1}^{n}\bigl(x_i\,y_{i+1}-x_{i+1}\,y_i\bigr)\right|, \tag{3.1}\]
and its centroid \((C_x,C_y)\) is
\[ C_x=\frac{1}{6A_s}\sum_{i=1}^{n}(x_i+x_{i+1})(x_i y_{i+1}-x_{i+1}y_i),\quad C_y=\frac{1}{6A_s}\sum_{i=1}^{n}(y_i+y_{i+1})(x_i y_{i+1}-x_{i+1}y_i), \tag{3.2}\]
where \(A_s\) is the signed area (without the absolute value in Equation 3.1). These formulas are valid in a planar, projected coordinate system, which is why area calculations should not be performed on unprojected degrees.
3.2.3 Topological relations
The relationship between two geometries \(A\) and \(B\) is rigorously described by the dimensionally extended nine-intersection model (DE-9IM), a \(3\times3\) matrix recording the dimension of the intersection of the interior, boundary, and exterior of each geometry:
\[ \mathrm{DE\text{-}9IM}(A,B)= \begin{pmatrix} \dim(I_A\cap I_B) & \dim(I_A\cap B_B) & \dim(I_A\cap E_B)\\ \dim(B_A\cap I_B) & \dim(B_A\cap B_B) & \dim(B_A\cap E_B)\\ \dim(E_A\cap I_B) & \dim(E_A\cap B_B) & \dim(E_A\cap E_B) \end{pmatrix}. \tag{3.3}\]
Named predicates such as equals, disjoint, touches, within, contains, overlaps, and crosses are patterns of this matrix and underlie spatial queries and joins.
3.3 The raster data model
A raster is a regular tessellation of the plane into cells (pixels), each holding a value. It is the natural model for fields and for images. The link between cell indices (column \(c\), row \(r\)) and map coordinates \((x,y)\) is an affine transformation:
\[ \begin{pmatrix}x\\y\end{pmatrix}= \begin{pmatrix}x_0\\y_0\end{pmatrix}+ \begin{pmatrix}\Delta x & \beta_x\\ \beta_y & -\Delta y\end{pmatrix} \begin{pmatrix}c\\r\end{pmatrix}, \tag{3.4}\]
where \((x_0,y_0)\) is the coordinate of the upper-left corner, \(\Delta x\) and \(\Delta y\) are the cell sizes, and the off-diagonal terms \(\beta_x,\beta_y\) vanish for north-up rasters. A raster is therefore defined by its extent, cell size, CRS, data type (bit depth), number of bands, and a NoData value that flags missing observations.
The uncompressed size of a raster is
\[ S=\frac{n_{\text{rows}}\;n_{\text{cols}}\;n_{\text{bands}}\;b}{8}\ \ \text{bytes}, \tag{3.5}\]
with \(b\) the bits per sample. Practical systems reduce this cost through lossless compression (LZW, DEFLATE), tiling (for random access), overviews (image pyramids for fast display), and cloud-optimised layouts that permit reading only the required window over a network.
Discrete versus continuous rasters. Categorical rasters (land cover) must be resampled with nearest-neighbour rules to avoid inventing nonexistent classes, whereas continuous rasters (reflectance, elevation) may be interpolated. This distinction returns in Chapter 6.
3.4 Surfaces: grids, TINs, and point clouds
Elevation may be stored as a digital elevation model (DEM) grid, a triangulated irregular network (TIN), or a point cloud. A TIN connects irregularly spaced points into non-overlapping triangles, most commonly according to the Delaunay criterion: the circumcircle of every triangle contains no other vertex, which maximises the minimum interior angle and avoids thin slivers. TINs adapt density to terrain complexity, whereas grids offer simplicity and compatibility with map algebra. A digital surface model (DSM) includes vegetation and buildings; a digital terrain model (DTM) represents bare ground.
3.5 Attributes, databases, and spatial indexes
Attribute data are typically stored in relational tables linked to geometries by keys. The relational operations of selection, projection, and join (one-to-one, many-to-one, one-to-many) combine tables, while spatial joins combine tables through topological predicates (for example, assigning each health-facility point the attributes of the district polygon that contains it).
Large datasets are made searchable by spatial indexes: R-trees group nearby objects into nested bounding rectangles, quadtrees recursively subdivide space into quadrants, and grid or geohash indexes map locations to cell identifiers. An index does not change the answer to a query; it changes how fast it can be found.
3.6 Data quality, metadata, and interoperability
The international standard for geographic-information quality (ISO 19157) describes five elements: positional accuracy, thematic (attribute) accuracy, logical consistency, completeness (omission and commission), and temporal quality, plus lineage (the history of processing). Metadata document these properties and enable discovery. For satellite archives, SpatioTemporal Asset Catalogs (STAC) provide a common language for describing imagery, while the FAIR principles (findable, accessible, interoperable, reusable) guide data stewardship.
3.7 Generalisation and scale dependence
Cartographic generalisation adapts detail to scale through operations such as simplification, aggregation, displacement, and exaggeration. Analytic results are scale dependent: the measured length of a coastline increases as the measurement unit decreases, a phenomenon captured by the fractal-dimension relation \(L(\varepsilon)\propto\varepsilon^{\,1-D}\), where \(D>1\) is the fractal dimension and \(\varepsilon\) the measurement step. Therefore lengths and areas should be reported with the scale and resolution at which they were measured.
Exercises
Exercise 3.1 A four-band, 16-bit multispectral raster has 10 000 rows and 8 000 columns with 30 m cells. (a) Compute its uncompressed size using Equation 3.5, in megabytes (decimal) and mebibytes. (b) Compute the ground area covered in km². (c) If the cell size were refined to 10 m over the same extent, by what factor does storage increase?
Exercise 3.2 A parcel has vertices \((1,1)\), \((6,2)\), \((5,6)\), and \((2,5)\) in metres. Use Equation 3.1 to compute its area. Then explain why the same calculation performed on longitude–latitude values would be meaningless.
Exercise 3.3 Two polygons share a common edge but do not overlap. Write the pattern of the DE-9IM matrix that distinguishes “touches” from “overlaps,” and describe a practical query in a health-facility catchment analysis in which the difference matters.
Exercise 3.4 Discuss the advantages and disadvantages of representing a river network as (a) a raster of flow-accumulation values, (b) a set of polylines with attributes, and (c) a topological network graph. For each, name one analysis that is easy and one that is difficult.
Exercise 3.5 The measured length of a coastline is 1200 km using a 10 km step and 1890 km using a 2 km step. Estimate the fractal dimension \(D\) from \(L(\varepsilon)\propto\varepsilon^{\,1-D}\). What does the result imply for comparing coastline lengths from two sources?
4 Physics of Remote Sensing
4.1 Electromagnetic radiation
Electromagnetic (EM) radiation transports energy as waves that travel at the speed of light \(c\approx 2.998\times10^{8}\ \mathrm{m\,s^{-1}}\) in vacuum. Wavelength \(\lambda\), frequency \(f\), and photon energy \(E_{\gamma}\) are related by
\[ c=\lambda f,\qquad E_{\gamma}=h f=\frac{h c}{\lambda}, \tag{4.1}\]
where \(h=6.626\times10^{-34}\ \mathrm{J\,s}\) is Planck’s constant. Shorter wavelengths therefore carry more energy per photon, which is why gamma rays are damaging and why the detection of long-wave microwave emission requires large collecting areas.
Remote sensing exploits several regions of the spectrum (Figure 4.1): the visible (0.4–0.7 µm), near-infrared (NIR, 0.7–1.3 µm), shortwave infrared (SWIR, 1.3–3 µm), mid- and thermal infrared (3–14 µm) and the microwave region (about 1 mm to 1 m). The choice of region determines what the sensor can perceive: pigments and surface colour in the visible, leaf structure in the NIR, water content and minerals in the SWIR, temperature in the thermal infrared, and roughness, structure, and moisture in the microwave.
4.2 Radiometric quantities
Remote sensing measures radiance, but many applications need reflectance. The following quantities form the vocabulary of radiometry.
| Quantity | Symbol | Definition | Unit |
|---|---|---|---|
| Radiant flux | \(\Phi\) | energy per unit time | W |
| Irradiance | \(E\) | flux incident on a surface, per unit area | W m\(^{-2}\) |
| Radiant exitance | \(M\) | flux leaving a surface, per unit area | W m\(^{-2}\) |
| Radiance | \(L\) | flux per unit projected area per unit solid angle | W m\(^{-2}\) sr\(^{-1}\) |
| Spectral radiance | \(L_\lambda\) | radiance per unit wavelength | W m\(^{-2}\) sr\(^{-1}\) µm\(^{-1}\) |
| Reflectance | \(\rho\) | ratio of reflected to incident flux | unitless |
A surface is Lambertian if its radiance is independent of viewing direction. For such a surface, exitance and radiance are related by \(M=\pi L\), and the reflectance is
\[ \rho=\frac{\pi\,L}{E}. \tag{4.2}\]
Real surfaces are not Lambertian. Their directional behaviour is described by the bidirectional reflectance distribution function (BRDF),
\[ f_r(\theta_i,\phi_i;\theta_v,\phi_v)=\frac{\mathrm{d}L_v(\theta_v,\phi_v)}{\mathrm{d}E_i(\theta_i,\phi_i)}\quad[\mathrm{sr^{-1}}], \tag{4.3}\]
which is the reason that the same forest looks brighter when viewed from the direction of the Sun (the hotspot) than from the opposite direction, and why sensors with wide fields of view exhibit view-angle effects.
4.3 Blackbody radiation and thermal emission
Every object above absolute zero emits radiation. A blackbody, a perfect absorber and emitter, at absolute temperature \(T\) emits spectral exitance given by Planck’s law:
\[ M_\lambda(T)=\frac{2\pi h c^{2}}{\lambda^{5}}\;\frac{1}{\exp\!\left(\dfrac{h c}{\lambda k T}\right)-1}, \tag{4.4}\]
with Boltzmann’s constant \(k=1.381\times10^{-23}\ \mathrm{J\,K^{-1}}\). Two consequences follow. Wien’s displacement law gives the wavelength of maximum emission,
\[ \lambda_{\max}=\frac{2898\ \mathrm{\mu m\,K}}{T}, \tag{4.5}\]
and the Stefan–Boltzmann law gives the total exitance,
\[ M=\sigma T^{4},\qquad \sigma=5.670\times10^{-8}\ \mathrm{W\,m^{-2}\,K^{-4}}. \tag{4.6}\]
The Sun (\(T\approx5800\) K) peaks near 0.5 µm, in the visible; the Earth’s surface (\(T\approx300\) K) peaks near 10 µm, in the thermal infrared. The two emissions barely overlap, which is why remote sensing of reflected sunlight (0.4–3 µm) and remote sensing of emitted heat (3–14 µm) can be treated separately.
Real materials emit less than a blackbody by a factor called emissivity \(\varepsilon_\lambda\in[0,1]\), so that \(L_\lambda=\varepsilon_\lambda B_\lambda(T)\), where \(B_\lambda=M_\lambda/\pi\). By Kirchhoff’s law, for an opaque surface in thermal equilibrium, \(\varepsilon_\lambda=1-\rho_\lambda\): good reflectors are poor emitters. The temperature that a blackbody would need to produce the observed radiance is the brightness temperature. Inverting Equation 4.4 for a narrow band gives the form used operationally,
\[ T_b=\frac{K_2}{\ln\!\left(\dfrac{K_1}{L_\lambda}+1\right)}, \tag{4.7}\]
where \(K_1\) and \(K_2\) are sensor-specific calibration constants. Land-surface temperature requires additionally correcting for emissivity and atmospheric effects, most commonly through split-window algorithms that combine two adjacent thermal bands (approximately 10.8 and 12 µm) in a linear function of their brightness temperatures.
4.4 Interactions with the atmosphere
Radiation crossing the atmosphere is attenuated by absorption (chiefly water vapour, carbon dioxide, ozone, and oxygen) and redirected by scattering:
- Rayleigh scattering by molecules much smaller than the wavelength varies as \(\lambda^{-4}\), which explains why the sky is blue and why short-wave (blue) bands suffer the greatest haze.
- Mie scattering by particles comparable in size to the wavelength (aerosols, smoke, fine dust) varies approximately as \(\lambda^{-1}\) to \(\lambda^{-2}\).
- Non-selective scattering by particles much larger than the wavelength (cloud droplets) scatters all visible wavelengths about equally, making clouds white.
Direct transmission through a homogeneous layer follows the Beer–Lambert law:
\[ \tau=\exp\!\left(-\frac{\delta}{\cos\theta}\right)=\exp(-\delta\,m), \tag{4.8}\]
where \(\delta\) is the vertical optical depth, \(\theta\) the zenith angle of the path, and \(m=1/\cos\theta\) the air mass in the plane-parallel approximation. Oblique paths are longer and more strongly attenuated.
4.4.1 The at-sensor radiance
For a flat, uniform Lambertian target of surface reflectance \(\rho\) and ignoring adjacency and coupling terms, the radiance reaching a sensor in a reflective band may be written
\[ L_{\text{sensor}}=L_{\text{path}}+\tau_v\,\frac{\rho}{\pi}\,E_g,\qquad E_g=E_0\cos\theta_s\,\tau_s+E_{\text{diff}}, \tag{4.9}\]
where \(L_{\text{path}}\) is the radiance scattered by the atmosphere directly into the sensor (path radiance), \(\tau_v\) and \(\tau_s\) are the transmittances along the view and Sun paths, \(E_0\) is the exo-atmospheric solar irradiance, \(\theta_s\) the solar zenith angle, and \(E_{\text{diff}}\) the diffuse sky irradiance. Atmospheric correction (Chapter 6) is the inversion of Equation 4.9 for \(\rho\):
\[ \rho=\frac{\pi\,(L_{\text{sensor}}-L_{\text{path}})}{\tau_v\,E_g}. \tag{4.10}\]
4.5 Interactions with the target
Incident energy is partitioned among reflection, absorption, and transmission: \[ E_I(\lambda)=E_R(\lambda)+E_A(\lambda)+E_T(\lambda). \tag{4.11}\]
The spectral reflectance curve of a material, its spectral signature, is the basis of discrimination. Figure 4.2 shows the classic signatures of three cover types.
Vegetation. In the visible, chlorophyll and carotenoids absorb blue (about 0.45 µm) and red (about 0.67 µm) light, leaving a small green reflectance peak near 0.55 µm. At approximately 0.70–0.75 µm the reflectance rises abruptly, the red edge, to a high NIR plateau caused by multiple scattering at the interfaces between cell walls and air spaces in the leaf mesophyll. In the SWIR, reflectance is governed by leaf water, with absorption features near 1.45 and 1.95 µm. Stress, senescence, or disease reduce chlorophyll and water content and shift the red edge to shorter wavelengths.
Soils. Soil spectra are smoother and generally increase with wavelength. Moisture lowers reflectance throughout, most strongly in the water-absorption bands; organic matter darkens the visible; iron oxides impart the characteristic red–yellow colour and absorb in the blue and near-UV; clay minerals exhibit an absorption near 2.2 µm.
Water. Clear water absorbs strongly beyond 0.7 µm and is nearly black in the NIR and SWIR, which makes those bands ideal for delineating water bodies. In the visible, the reflectance depends on suspended sediments (raising green and red), chlorophyll (absorbing blue and red), and dissolved organic matter (absorbing blue).
Microwave interaction. The backscatter of radar depends on the target’s geometry and dielectric properties, not on pigment or temperature. A surface appears smooth to radar if its root-mean-square height \(h_{\mathrm{rms}}\) satisfies the Rayleigh-type criterion \[ h_{\mathrm{rms}}<\frac{\lambda}{8\cos\theta_i}, \tag{4.12}\] where \(\theta_i\) is the incidence angle. Smooth surfaces reflect energy away from the sensor (dark), rough surfaces scatter diffusely (bright). Because water has a high dielectric constant, wet soils and flooded vegetation are strong reflectors (Chapter 9).
Exercises
Exercise 4.1 Compute the peak emission wavelength for (a) the Sun (\(T=5800\) K) and (b) the Earth’s surface (\(T=300\) K). Using the Stefan–Boltzmann law, compute the exitance of a blackbody at 300 K, and the ratio of the Sun’s exitance to the Earth’s.
Exercise 4.2 Compute the photon energy, in joules and electron-volts (\(1\ \mathrm{eV}=1.602\times10^{-19}\) J), at \(\lambda=0.55\) µm. What is the frequency, in GHz, of a C-band radar with wavelength 5.6 cm?
Exercise 4.3 The vertical optical depth of an atmospheric column at 0.55 µm is \(\delta=0.2\). Compute the transmittance for (a) a vertical path, and (b) a path at a \(60^{\circ}\) zenith angle. By what percentage is the transmitted signal at \(60^{\circ}\) lower than the vertical value?
Exercise 4.4 A thermal band has \(K_1=774.89\) and \(K_2=1321.08\) (in the units used for radiance in \(\mathrm{W\,m^{-2}\,sr^{-1}\,\mu m^{-1}}\) and kelvin, respectively). Compute the brightness temperature for a measured radiance of 10.0. Explain why the land-surface temperature is generally different from \(T_b\).
Exercise 4.5 A false-colour composite assigns NIR to red, red to green, and green to blue. Predict, and explain using Figure 4.2, the colours of (a) dense healthy forest, (b) bare dry soil, and (c) clear deep water.
Exercise 4.6 An L-band SAR (\(\lambda=24\) cm) images a ploughed field with \(h_{\mathrm{rms}}=2\) cm at \(\theta_i=35^{\circ}\); a C-band SAR (\(\lambda=5.6\) cm) images the same field. Using Equation 4.12, determine whether the field appears smooth or rough at each wavelength, and state the consequence for backscatter.
5 Sensors, Platforms, and Image Acquisition
5.1 Platforms and orbits
Sensors are carried on ground-based masts, unmanned aerial vehicles (UAVs), aircraft, balloons, and satellites. The choice of platform trades spatial detail against coverage, flexibility against repeatability, and cost against scale.
5.1.1 Orbital mechanics
For a circular orbit of semi-major axis \(a=R_E+H_o\) (Earth radius plus altitude), Kepler’s third law gives the period
\[ T=2\pi\sqrt{\frac{a^{3}}{\mu}},\qquad \mu=GM_E=3.986\times10^{14}\ \mathrm{m^{3}\,s^{-2}} . \tag{5.1}\]
The principal orbit classes are shown in Figure 5.1.
- Geostationary orbit (GEO), at about 35 786 km altitude and zero inclination, has a period equal to Earth’s rotation, so the satellite remains above one point on the equator. It provides continuous monitoring (minute-scale repeat) of a hemisphere but at coarse resolution and poor viewing geometry at high latitudes. Weather and rapid-change monitoring use GEO.
- Sun-synchronous polar orbits at 500–850 km altitude are the workhorse of land observation. Their inclination is chosen so that the orbital plane precesses eastward at \(360^{\circ}\) per tropical year (\(\approx 0.9856^{\circ}\,\mathrm{day^{-1}}\)), matching the Earth’s motion around the Sun. This is achieved by exploiting the Earth’s equatorial bulge (the \(J_2\) term): \[ \dot{\Omega}=-\frac{3}{2}\,J_2\left(\frac{R_E}{a(1-e^{2})}\right)^{2}n\cos i,\qquad n=\sqrt{\mu/a^{3}}, \tag{5.2}\] so that \(\dot\Omega=+0.9856^{\circ}\,\mathrm{day^{-1}}\) requires a retrograde inclination \(i\approx 97\)–\(99^{\circ}\). A satellite in such an orbit crosses the equator at the same local solar time on every pass (for example, about 10:00–10:30 for many optical missions), giving consistent illumination that is essential for time-series analysis.
- Non-sun-synchronous low Earth orbits (for example, that of the International Space Station, inclination about \(52^{\circ}\)) sample all local times and are used by some LiDAR and thermal missions.
The revisit time is the interval between successive observations of the same location; it depends on the orbit’s repeat cycle, the sensor’s swath width, and any off-nadir pointing capability.
5.2 Sensor architectures
- Framing systems (film or digital frame cameras) record an entire scene at an instant; geometry is that of a central perspective and is described by the collinearity equations of photogrammetry.
- Across-track (whiskbroom) scanners use a rotating mirror to sweep a small number of detectors across the swath; each pixel is observed in sequence.
- Along-track (pushbroom) scanners use a linear array of thousands of detectors perpendicular to the flight direction, “pushed” along the track. Because each detector views its ground cell for longer, pushbrooms achieve better signal-to-noise ratio and geometric stability. Most modern optical satellites use pushbroom designs.
- Microwave radiometers passively measure thermal microwave emission (soil moisture, sea-surface temperature, snow water equivalent).
- Active systems: imaging radar (Chapter 9), LiDAR (Chapter 9), and radar altimeters.
5.3 The four resolutions
The information content of imagery is characterised by four resolutions (Figure 5.2).
5.3.1 Spatial resolution
The instantaneous field of view (IFOV) is the angular cone of visibility of a single detector element. For a sensor at altitude \(H_o\) with IFOV \(\beta\) (radians) looking at nadir, the ground sample distance (GSD) is
\[ \mathrm{GSD}=2H_o\tan\!\left(\frac{\beta}{2}\right)\approx H_o\,\beta . \tag{5.3}\]
For a camera with detector pitch \(p\) and focal length \(f\), \(\beta=p/f\) and thus \(\mathrm{GSD}=H_o\,p/f\). The swath width is set by the total field of view \(\Omega\): \(W=2H_o\tan(\Omega/2)\). The effective resolution is also limited by the optics (diffraction: the smallest resolvable angle is \(\theta\approx1.22\,\lambda/D\) for an aperture \(D\)), the platform motion, and atmospheric blurring, and is described by the modulation transfer function (MTF). GSD and effective resolution are therefore not the same.
The finer the resolution, the smaller the swath and the larger the data volume; moreover, finer resolution is not always better, since within-class variance increases and the pixel-based classification of heterogeneous classes such as urban land may become harder (Chapter 8).
5.3.2 Spectral resolution
Spectral resolution is the number and width of the bands, quantified by the full width at half maximum (FWHM) of each band’s response function. Multispectral sensors record 3–15 broad bands; hyperspectral sensors record tens to hundreds of narrow (about 5–10 nm), contiguous bands, sufficient to resolve diagnostic absorption features.
5.3.3 Radiometric resolution
Radiometric resolution is the sensitivity to small differences in radiance, determined by the bit depth \(b\) and the noise. A sensor quantising to \(b\) bits records \(2^{b}\) digital numbers (DN); for example, 8 bits give 256 levels and 12 bits give 4096. The signal-to-noise ratio (SNR), \(\mathrm{SNR}=\bar{S}/\sigma_n\), is the more informative quantity: quantisation beyond the noise floor adds no information.
5.3.4 Temporal resolution
The temporal resolution is the interval at which the same target is observed. Higher temporal resolution supports monitoring of vegetation phenology, floods, and disturbance, but at the coarse spatial scales that wide swaths allow. Constellations of many small satellites, and the combination of multiple missions (“virtual constellations”), reduce the trade-off.
5.4 Representative missions
Missions and their operational status change over time; Table 5.1 lists representative characteristics and should be verified against current mission documentation.
| System | Type | Spatial resolution | Swath | Revisit | Typical use |
|---|---|---|---|---|---|
| Landsat 8/9 (OLI/TIRS) | Multispectral, thermal | 15 m pan, 30 m MS, 100 m thermal | 185 km | 16 d each; 8 d combined | Land cover and change, long records |
| Sentinel-2 (MSI) | Multispectral, 13 bands | 10, 20, 60 m | 290 km | about 5 d (constellation) | Vegetation, agriculture, land cover |
| MODIS / VIIRS | Moderate-resolution | 250 m to 1 km (MODIS); 375–750 m (VIIRS) | about 2300–3000 km | daily | Global monitoring, fire, phenology |
| Sentinel-1 | C-band SAR | about 5 × 20 m (IW mode) | 250 km | 6–12 d | Floods, deformation, all-weather mapping |
| PlanetScope | Multispectral cubesat constellation | about 3–4 m | strips | near-daily | Rapid change monitoring |
| Very-high-resolution commercial (e.g., WorldView class) | Multispectral, pan | sub-metre pan | 10–20 km | 1–4 d (tasked) | Urban, infrastructure, detailed mapping |
| Geostationary imagers (e.g., ABI class) | Multispectral | 0.5–2 km | full disk | 10–15 min | Weather, fires, rapid change |
5.5 Data processing levels and analysis-ready data
Satellite products are commonly distributed at increasing levels of processing: Level 0 (raw telemetry), Level 1 (radiometrically calibrated, geolocated radiance or top-of-atmosphere reflectance), Level 2 (geophysical variables such as surface reflectance or land-surface temperature), Level 3 (spatially and temporally aggregated, gridded composites), and Level 4 (model-assimilated or gap-filled fields). Analysis-ready data (ARD) are Level-2 products that have been consistently calibrated, atmospherically and geometrically corrected, gridded, and accompanied by quality masks, so that analysts can proceed directly to analysis.
Exercises
Exercise 5.1 Compute the orbital period of a satellite at 705 km altitude (\(R_E=6378\) km) and the number of orbits per day. The Earth rotates about \(24.7^{\circ}\) during one orbit, so successive ground tracks at the equator are separated by that longitude difference (1° ≈ 111.3 km). Compute the separation, the fraction of it covered by a 185 km swath, and hence roughly how many days are needed before the swaths of successive days fill the gaps.
Exercise 5.2 A pushbroom sensor at 705 km altitude has an IFOV of 42.5 µrad. (a) Compute the GSD. (b) The sensor swath is 185 km; compute the total field of view in degrees. (c) A UAV camera has a pixel pitch of 2.4 µm, a focal length of 8.8 mm, and flies at 100 m above ground: compute its GSD.
Exercise 5.3 Explain why it is physically difficult to build a satellite that has simultaneously a very fine spatial resolution, a very wide swath, many narrow spectral bands, and daily global revisit. Which of these are typically relaxed in (a) a geostationary weather imager; (b) a global land-monitoring mission; (c) a commercial very-high-resolution satellite?
Exercise 5.4 A 30 m pixel contains 40 % dense forest (\(\rho_{\text{NIR}}=0.50\)) and 60 % bare soil (\(\rho_{\text{NIR}}=0.30\)). Assuming linear mixing, compute the pixel’s NIR reflectance. How would the pixel change if the true composition were 50 %/50 %, and what does this imply for the accuracy of forest-area estimates from thresholding coarse pixels?
Exercise 5.5 Design problem. Select and justify an imaging system (platform, bands, resolution, revisit) to monitor smallholder maize fields (typical size 0.5 ha) through a growing season under frequent cloud cover. State what you would give up.
6 Image Preprocessing
Raw sensor data are not measurements of the surface. Between the surface and the analyst lie the atmosphere, the sensor’s response, the platform’s motion, the terrain, and the illumination geometry. Preprocessing removes or reduces these effects so that pixel values represent the same physical quantity, in the same place, at all times.
6.1 Radiometric calibration
6.1.1 From digital numbers to radiance
Sensor detectors record digital numbers linearly related to at-sensor spectral radiance,
\[ L_\lambda=G\cdot \mathrm{DN}+B, \tag{6.1}\]
where the gain \(G\) and offset (bias) \(B\) are supplied in the metadata, derived from pre-launch and on-orbit calibration. Equivalent forms use the minimum and maximum radiance and the range of quantised values: \(L_\lambda=\frac{L_{\max}-L_{\min}}{Q_{\max}-Q_{\min}}(\mathrm{DN}-Q_{\min})+L_{\min}\). Systematic detector artefacts (striping, dropped lines, bad pixels) are removed by destriping (matching the detector histograms) or by interpolation.
6.1.2 From radiance to top-of-atmosphere reflectance
Comparison across dates and sensors requires normalising for the solar irradiance and geometry. The top-of-atmosphere (TOA) reflectance is
\[ \rho_{\text{TOA}}=\frac{\pi\,L_\lambda\,d^{2}}{E_{\mathrm{sun},\lambda}\cos\theta_s}, \tag{6.2}\]
where \(d\) is the Earth–Sun distance in astronomical units, \(E_{\mathrm{sun},\lambda}\) the mean exo-atmospheric solar irradiance in the band, and \(\theta_s\) the solar zenith angle (90° minus the sun elevation). A widely used approximation for \(d\) on day of year \(D\) is \(d\approx 1-0.01672\cos\bigl(0.9856\,(D-4)^{\circ}\bigr)\).
Example 6.1 (Converting a digital number to TOA reflectance) A pixel in a red band has \(\mathrm{DN}=112\). The calibration gives \(G=0.80\) and \(B=-4.0\) (radiance units \(\mathrm{W\,m^{-2}\,sr^{-1}\,\mu m^{-1}}\)). The image was acquired in mid-January (\(d=0.9833\) AU) with solar zenith angle \(\theta_s=40^{\circ}\), and \(E_{\mathrm{sun}}=1550\ \mathrm{W\,m^{-2}\,\mu m^{-1}}\).
Radiance: \(L=0.80(112)-4.0=85.6\).
Reflectance: \[\rho_{\text{TOA}}=\frac{\pi(85.6)(0.9833)^{2}}{1550\cos 40^{\circ}}=\frac{260.0}{1187.4}\approx 0.219 .\]
6.2 Atmospheric correction
TOA reflectance still contains the effects of atmospheric scattering and absorption (Equation 4.9). Atmospheric correction retrieves surface reflectance, the property of the target itself. Methods form a spectrum of increasing rigour:
- Image-based (relative) methods. Dark-object subtraction (DOS) assumes that some pixels (deep clear water, shadow, dense dark vegetation) should have near-zero reflectance, so the minimum value in each band estimates the additive path radiance, which is subtracted: \(\rho_{\text{corr}}=\rho_{\text{TOA}}-\rho_{\text{haze}}\). It is simple but neglects transmittance, and it is unreliable where dark objects are absent.
- Empirical line calibration. Field-measured reflectance \(\rho\) of at least two calibration targets (one bright, one dark) is regressed against image values, \(\rho=a\,\mathrm{DN}+b\) per band, and the fitted line is applied to the whole image.
- Radiative-transfer-based methods. Codes such as 6S, MODTRAN, and operational processors (for example, LaSRC for Landsat and Sen2Cor or MAJA for Sentinel-2) solve Equation 4.9 using estimates of the aerosol optical depth, column water vapour, ozone, and surface pressure, either from ancillary data or from the image itself, to invert Equation 4.10.
Corrections must also treat the adjacency effect (radiance from bright surroundings scattered into a dark pixel’s line of sight), which is important near shorelines and in heterogeneous landscapes. Absolute correction is essential when combining images from different dates or sensors, or when using physical models; it may be unnecessary for a single-date classification with training data drawn from the same image.
6.3 Topographic correction
In mountainous terrain, the illumination angle varies with slope \(S\) and aspect \(A\). The cosine of the local illumination angle \(i\) is
\[ \cos i=\cos\theta_s\cos S+\sin\theta_s\sin S\cos(\phi_s-A), \tag{6.3}\]
where \(\phi_s\) is the solar azimuth. Slopes facing the Sun (\(\cos i\) large) appear bright, those facing away appear dark, even for identical cover. Empirical corrections normalise reflectance to a flat, horizontal surface:
\[ \rho_{h}=\rho\,\frac{\cos\theta_s}{\cos i}\ \ \text{(cosine)},\qquad \rho_{h}=\rho\,\frac{\cos\theta_s+C}{\cos i+C}\ \ \text{(C-correction)},\qquad C=\frac{b}{m}, \tag{6.4}\]
where \(b\) and \(m\) are the intercept and slope of the regression \(\rho=b+m\cos i\) fitted per band and cover class. The simple cosine correction over-corrects poorly illuminated slopes; the \(C\) parameter damps this. Physical and semi-empirical alternatives (Minnaert, sun–canopy–sensor) are used where forest canopy geometry matters. A DEM of appropriate resolution and accuracy is prerequisite.
6.4 Geometric correction
6.4.1 Sources of geometric distortion
Distortions are systematic and predictable (Earth rotation during scanning, panoramic effects in wide-field scanners, sensor optics, orbit) or random (platform attitude, altitude, and velocity variations). Relief displacement is the radial displacement of an elevated point from its true planimetric position in a central-perspective image:
\[ d=\frac{r\,h}{H_p}, \tag{6.5}\]
where \(r\) is the radial distance from the nadir point to the displaced image point, \(h\) the height of the feature above the datum, and \(H_p\) the flying height above the datum. Displacement grows with distance from nadir and with height, and is removed by orthorectification using a DEM.
6.4.2 Polynomial and rigorous models
The classical procedure has three stages: (i) transformation model fitting from ground control points (GCPs), features identifiable in both the image and a reference (map, orthoimage, GNSS survey); (ii) accuracy evaluation; and (iii) resampling. A polynomial transformation of order \(n\) relates image coordinates \((x,y)\) to reference coordinates \((x',y')\):
\[ x'=\sum_{i=0}^{n}\sum_{j=0}^{n-i}a_{ij}\,x^{i}y^{j},\qquad y'=\sum_{i=0}^{n}\sum_{j=0}^{n-i}b_{ij}\,x^{i}y^{j}, \tag{6.6}\]
with \((n+1)(n+2)/2\) coefficients per axis (3, 6, 10 for orders 1, 2, 3). At least as many GCPs as coefficients are required, and more are needed for a robust least-squares solution \(\hat{\mathbf{a}}=(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{x}'\). High-order polynomials oscillate between control points and should be avoided. The fit is judged by the root-mean-square error (RMSE) of residuals \((\Delta x_k,\Delta y_k)\) at \(m\) control points,
\[ \mathrm{RMSE}=\sqrt{\frac{1}{m}\sum_{k=1}^{m}\left(\Delta x_k^{2}+\Delta y_k^{2}\right)}, \tag{6.7}\]
and, crucially, must be verified on independent check points not used in the fit. A registration RMSE of 0.5 pixel or better is the usual requirement for change detection. For high-resolution satellites, rational polynomial coefficients (RPCs), ratios of cubic polynomials in ground coordinates, provide a rigorous sensor model that can be refined with a few GCPs, and orthorectification combines the sensor model with a DEM.
6.4.3 Resampling
Geometric correction maps output pixel centres back to fractional positions in the input image, whose values must be estimated by resampling:
- Nearest neighbour copies the closest input value; it preserves original values (required for categorical rasters) but produces blocky, spatially jittered edges.
- Bilinear interpolation uses the four nearest pixels. For a point at fractional offset \((u,v)\in[0,1]^2\) within a cell whose corner values are \(f_{00},f_{10},f_{01},f_{11}\), \[ f(u,v)=(1-u)(1-v)f_{00}+u(1-v)f_{10}+(1-u)v\,f_{01}+uv\,f_{11}. \tag{6.8}\]
- Cubic convolution uses a \(4\times4\) neighbourhood; it gives sharper results but may overshoot at edges and take values outside the original range.
6.5 Mosaicking, compositing, and cloud masking
Adjacent scenes are mosaicked after radiometric harmonisation (matching overlap statistics) and by seam-line placement to minimise visible discontinuities. Clouds and cloud shadows are identified by thresholding brightness, temperature, and spectral tests (for example, Fmask), by machine-learning classifiers, or from the quality bands accompanying ARD products, and must be excluded before analysis. Temporal compositing builds a nearly cloud-free image from many acquisitions by selecting, per pixel, the observation that satisfies a criterion (maximum NDVI, minimum view angle) or by taking a robust statistic such as the median.
Exercises
Exercise 6.1 A green band has \(G=0.75\), \(B=-3.5\), \(E_{\mathrm{sun}}=1830\ \mathrm{W\,m^{-2}\,\mu m^{-1}}\). A pixel has \(\mathrm{DN}=96\). The image was taken on 21 June (\(d\approx1.0163\) AU) at solar elevation \(62^{\circ}\). Compute the TOA reflectance. (Use \(d^2\) in Equation 6.2.)
Exercise 6.2
- How many GCPs are the minimum for a second-order and a third-order polynomial transformation? (b) Four GCPs have residuals \((\Delta x,\Delta y)\) of \((0.3,0.4)\), \((-0.5,0.2)\), \((0.1,-0.6)\) and \((-0.2,-0.1)\) pixels. Compute the RMSE. (c) Why is RMSE computed at the GCPs an overly optimistic estimate of accuracy?
Exercise 6.3 Four neighbouring pixel values are \(f_{00}=10\), \(f_{10}=20\), \(f_{01}=30\), \(f_{11}=40\). Compute the bilinear value at \((u,v)=(0.25,0.75)\). Which resampling method would you use for a land-cover map, and why?
Exercise 6.4 A vertical aerial photograph is taken from 3000 m above the datum. A radio tower 100 m tall appears at a radial distance of 80 mm from the nadir point. Compute the displacement of the tower’s top from its base in the photograph.
Exercise 6.5 On a slope of \(25^{\circ}\) facing away from the Sun, with \(\theta_s=45^{\circ}\), the solar azimuth differs from the aspect by \(180^{\circ}\). Compute \(\cos i\) and the multiplicative factor of the cosine correction. Why might this correction produce unrealistically bright pixels?
Exercise 6.6 Conceptual. For which of the following tasks is absolute atmospheric correction necessary, and for which is it unnecessary? Justify each. (a) Single-date supervised classification of land cover. (b) Comparing NDVI between two dates. (c) Deriving chlorophyll concentration from water-leaving reflectance. (d) Mosaicking two adjacent same-day scenes.
7 Image Enhancement and Transformations
Enhancement improves the visual interpretability of imagery or derives new, more informative layers from the original bands. Unlike the corrections of Chapter 6, enhancement can be irreversible and is chosen to serve a specific analytical or visual purpose.
7.1 Histograms and contrast enhancement
The histogram of a band records the frequency of each DN. Because sensors are designed to record the full range of radiances on Earth, a scene often occupies only a narrow part of the available range and appears dull. Linear contrast stretching maps the range \([\mathrm{DN}_{\min},\mathrm{DN}_{\max}]\) (or, more robustly, the 2nd to 98th percentiles) to the display range \([0,255]\):
\[ \mathrm{DN}'=255\;\frac{\mathrm{DN}-\mathrm{DN}_{\min}}{\mathrm{DN}_{\max}-\mathrm{DN}_{\min}} . \tag{7.1}\]
Histogram equalisation redistributes DNs to approximate a uniform histogram, using the cumulative distribution of the input. For \(L\) grey levels and normalised frequency \(p_r(r_j)=n_j/n\), the output level for input \(r_k\) is
\[ s_k=(L-1)\sum_{j=0}^{k}p_r(r_j). \tag{7.2}\]
Equalisation maximises contrast in the most populated ranges (often at the expense of rare but important values). All such displays alter appearance only; analytical operations should use the unstretched data.
7.2 Spatial filtering
Spatial filters exploit the values of neighbouring pixels. A linear filter convolves the image \(f\) with a kernel \(w\) of size \((2a+1)\times(2b+1)\):
\[ g(i,j)=\sum_{m=-a}^{a}\sum_{n=-b}^{b}w(m,n)\,f(i+m,\,j+n). \tag{7.3}\]
Low-pass filters smooth noise and emphasise regional trends. The \(3\times3\) mean and a Gaussian approximation are, respectively,
\[ \mathbf{W}_{\text{mean}}=\frac{1}{9}\begin{pmatrix}1&1&1\\1&1&1\\1&1&1\end{pmatrix},\qquad \mathbf{W}_{\text{Gauss}}=\frac{1}{16}\begin{pmatrix}1&2&1\\2&4&2\\1&2&1\end{pmatrix}. \]
The median filter (non-linear) replaces each pixel by the median of its neighbourhood; it removes salt-and-pepper noise while preserving edges. High-pass filters emphasise local detail; a second-derivative example is the Laplacian,
\[ \mathbf{W}_{\text{Lap}}=\begin{pmatrix}0&-1&0\\-1&4&-1\\0&-1&0\end{pmatrix}. \]
Edge detection uses first derivatives. The Sobel kernels estimate the horizontal and vertical gradients,
\[ \mathbf{W}_x=\begin{pmatrix}-1&0&1\\-2&0&2\\-1&0&1\end{pmatrix},\quad \mathbf{W}_y=\begin{pmatrix}-1&-2&-1\\0&0&0\\1&2&1\end{pmatrix},\quad |\nabla f|=\sqrt{G_x^{2}+G_y^{2}}, \tag{7.4}\]
with edge direction \(\arctan(G_y/G_x)\).
7.3 Spectral indices and band ratios
Band ratios and indices combine bands to emphasise a target property while suppressing illumination and topographic effects that scale all bands by a common factor. The best known is the Normalised Difference Vegetation Index,
\[ \mathrm{NDVI}=\frac{\rho_{\mathrm{NIR}}-\rho_{\mathrm{Red}}}{\rho_{\mathrm{NIR}}+\rho_{\mathrm{Red}}}, \tag{7.5}\]
which is bounded in \([-1,1]\) and increases with green-leaf area and chlorophyll. Because it depends on the ratio of reflectances, NDVI is largely insensitive to a multiplicative illumination change but sensitive to additive atmospheric effects and to the soil background.
Example 7.1 (Interpreting a pixel) A pixel has \(\rho_{\text{Red}}=0.06\) and \(\rho_{\text{NIR}}=0.42\). Then \(\mathrm{NDVI}=(0.42-0.06)/(0.42+0.06)=0.75\), indicating dense, healthy vegetation. Scaling both reflectances by 0.7 (for example, in shadow) leaves NDVI unchanged at 0.75.
The behaviour of NDVI is illuminated by the red–NIR feature space (Figure 7.1). Bare soils fall along a soil line \(\rho_{\mathrm{NIR}}=a\,\rho_{\mathrm{Red}}+b\); vegetation moves upward and to the left of the line as cover increases; lines of constant NDVI radiate from the origin.
Indices designed to reduce the soil influence or atmospheric effects include:
| Index | Formula | Target |
|---|---|---|
| SAVI | \(\dfrac{(1+L)(\rho_{\mathrm{NIR}}-\rho_{\mathrm{Red}})}{\rho_{\mathrm{NIR}}+\rho_{\mathrm{Red}}+L}\), with \(L\approx0.5\) | vegetation, reduced soil effect |
| EVI | \(G\dfrac{\rho_{\mathrm{NIR}}-\rho_{\mathrm{Red}}}{\rho_{\mathrm{NIR}}+C_1\rho_{\mathrm{Red}}-C_2\rho_{\mathrm{Blue}}+L_e}\); \(G=2.5,C_1=6,C_2=7.5,L_e=1\) | dense vegetation, less saturation, reduced aerosol effect |
| PVI | \(\dfrac{\rho_{\mathrm{NIR}}-a\rho_{\mathrm{Red}}-b}{\sqrt{1+a^{2}}}\) | vegetation, relative to the soil line |
| NDWI (McFeeters) | \(\dfrac{\rho_{\mathrm{Green}}-\rho_{\mathrm{NIR}}}{\rho_{\mathrm{Green}}+\rho_{\mathrm{NIR}}}\) | open water |
| MNDWI | \(\dfrac{\rho_{\mathrm{Green}}-\rho_{\mathrm{SWIR1}}}{\rho_{\mathrm{Green}}+\rho_{\mathrm{SWIR1}}}\) | water, suppresses built-up noise |
| NDMI | \(\dfrac{\rho_{\mathrm{NIR}}-\rho_{\mathrm{SWIR1}}}{\rho_{\mathrm{NIR}}+\rho_{\mathrm{SWIR1}}}\) | vegetation water content |
| NBR | \(\dfrac{\rho_{\mathrm{NIR}}-\rho_{\mathrm{SWIR2}}}{\rho_{\mathrm{NIR}}+\rho_{\mathrm{SWIR2}}}\) | burn scars |
| NDBI | \(\dfrac{\rho_{\mathrm{SWIR1}}-\rho_{\mathrm{NIR}}}{\rho_{\mathrm{SWIR1}}+\rho_{\mathrm{NIR}}}\) | built-up areas |
| NDSI | \(\dfrac{\rho_{\mathrm{Green}}-\rho_{\mathrm{SWIR1}}}{\rho_{\mathrm{Green}}+\rho_{\mathrm{SWIR1}}}\) | snow |
NDVI saturates at high biomass, is affected by soil colour and atmosphere, and its relationship with biophysical variables (leaf area index, fraction of absorbed radiation, biomass) is empirical and vegetation-type dependent. Indices from different sensors differ because band positions and widths differ, and must be cross-calibrated before being combined in time series.
7.4 Principal component analysis
Adjacent bands are highly correlated, so much of the data is redundant. Principal component analysis (PCA) rotates the \(B\)-dimensional feature space to a new orthogonal basis whose axes, the principal components, are ordered by the variance they explain. Let \(\boldsymbol{\Sigma}\) be the covariance matrix of the bands, with eigen-decomposition
\[ \boldsymbol{\Sigma}\,\mathbf{v}_k=\lambda_k\,\mathbf{v}_k,\qquad \lambda_1\ge\lambda_2\ge\dots\ge\lambda_B\ge0. \tag{7.6}\]
The \(k\)-th component of a pixel \(\mathbf{x}\) is the projection \[ y_k=\mathbf{v}_k^{\top}(\mathbf{x}-\boldsymbol{\mu}), \tag{7.7}\] and the proportion of the total variance explained is \(\lambda_k/\sum_j\lambda_j\). For multispectral imagery, PC1 typically expresses overall brightness or albedo, PC2 and PC3 capture contrasts such as vegetation versus soil, and the last components mostly noise. PCA is thus used for dimensionality reduction, noise suppression, and change detection. Because PCA is scene dependent, components are not directly comparable between images. The Minimum Noise Fraction (MNF) transform, which orders components by signal-to-noise ratio rather than variance, is preferred for hyperspectral data.
The Tasselled Cap transformation is a fixed, sensor-specific linear transformation \(\mathbf{y}=\mathbf{W}\mathbf{x}\) that rotates the bands to axes with agronomic meaning: brightness (soil), greenness (green vegetation), and wetness (soil and canopy moisture).
7.5 Pan-sharpening (image fusion)
Optical satellites often deliver a high-resolution panchromatic band and lower-resolution multispectral bands. Pan-sharpening injects spatial detail from the former into the latter. Common component-substitution approaches include:
- Brovey transform: each band is scaled by the ratio of the panchromatic value to the sum of the multispectral bands, \[ \mathrm{MS}_i'=\frac{\mathrm{MS}_i}{\sum_{j}\mathrm{MS}_j}\;\mathrm{PAN}. \tag{7.8}\]
- IHS, PCA, and Gram–Schmidt methods replace the intensity, first principal component, or a simulated panchromatic component with the high-resolution band.
- Multiresolution (wavelet, high-pass filter) methods add the high-frequency content of the panchromatic band to resampled multispectral bands, generally preserving spectra better.
Fused images are excellent for visual interpretation but their radiometry is altered; they should be used cautiously for quantitative analysis. Quality is judged by indices that compare fused with reference images (spectral angle, ERGAS, Q-index) both spectrally and spatially.
7.6 Texture
Texture, the spatial variation of tone within a neighbourhood, discriminates classes that share similar mean spectra (for example, orchards versus natural woodland). The grey-level co-occurrence matrix (GLCM) \(p(i,j\mid d,\theta)\) records how often a pixel with grey level \(i\) occurs at displacement \(d\) and direction \(\theta\) from a pixel with level \(j\). Haralick statistics summarise it:
\[ \text{Contrast}=\sum_{i,j}(i-j)^{2}p(i,j),\quad \text{Homogeneity}=\sum_{i,j}\frac{p(i,j)}{1+|i-j|},\quad \text{Entropy}=-\sum_{i,j}p(i,j)\ln p(i,j). \tag{7.9}\]
7.7 Fourier analysis
The two-dimensional discrete Fourier transform,
\[ F(u,v)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)\,e^{-\mathrm{i}\,2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right)}, \tag{7.10}\]
represents an image in terms of spatial frequencies. Periodic noise (such as striping) appears as isolated bright peaks in the frequency domain and can be removed by masking those peaks and applying the inverse transform (a notch filter). By the convolution theorem, convolution in the spatial domain equals multiplication in the frequency domain, which is how large kernels are applied efficiently.
Exercises
Exercise 7.1 Compute NDVI, SAVI (\(L=0.5\)), and EVI for a pixel with \(\rho_{\text{Blue}}=0.03\), \(\rho_{\text{Red}}=0.06\) and \(\rho_{\text{NIR}}=0.42\). Also compute NDWI for a pixel with \(\rho_{\text{Green}}=0.10\) and \(\rho_{\text{NIR}}=0.04\), and state what its sign implies.
Exercise 7.2
- Apply the \(3\times3\) mean filter to the window \(\begin{pmatrix}10&10&10\\10&50&10\\10&10&10\end{pmatrix}\) and give the new central value. (b) Apply the Sobel operators to \(\begin{pmatrix}1&2&3\\4&5&6\\7&8&9\end{pmatrix}\) and compute the gradient magnitude at the centre. (c) Which filter would you select to remove isolated erroneous pixels while preserving field boundaries?
Exercise 7.3 The covariance matrix of two bands is \(\boldsymbol{\Sigma}=\begin{pmatrix}4&2\\2&3\end{pmatrix}\). Find the eigenvalues, the percentage of variance explained by the first component, and the first eigenvector (normalised to unit length).
Exercise 7.4 Show algebraically that NDVI is unchanged when both \(\rho_{\mathrm{NIR}}\) and \(\rho_{\mathrm{Red}}\) are multiplied by a constant \(k>0\), but changes when a constant \(c\) is added to both. Relate this to the effect of atmospheric haze on NDVI.
Exercise 7.5 A pixel has multispectral values (Blue, Green, Red) = (40, 60, 80) and a panchromatic value of 120. Use the Brovey transform to compute the sharpened values. What is the physical limitation of ratio-based fusion for quantitative analysis?
Exercise 7.6 Design problem. You are asked to map open surface water in a landscape with dark building roofs and steep terrain shadows. Choose a set of indices, thresholds strategy, and ancillary layers, and explain how each addresses a specific source of false positives.
8 Image Classification and Accuracy Assessment
Classification converts continuous spectral measurements into thematic categories (land cover, crop type, mineral, disturbance class). It is the central inference problem of quantitative remote sensing, and its rigour is measured by how honestly its accuracy is assessed.
8.1 Framework and concepts
A classification scheme must be mutually exclusive, totally exhaustive, and, where possible, hierarchical. It distinguishes information classes (categories defined by the user, such as “maize”) from spectral classes (groups of pixels with similar spectra); the two are related many-to-many, and a fundamental task is to connect them.
A pixel is described by a feature vector \(\mathbf{x}=(x_1,\dots,x_B)^{\top}\) that may contain spectral bands, indices, texture, temporal metrics, and terrain covariates. The classifier is a function \(\hat{y}=g(\mathbf{x})\) assigning a label to each vector. Methods are categorised as:
- Unsupervised (clustering; no training labels) or supervised (learned from labelled samples);
- Parametric (assume a statistical distribution) or non-parametric;
- Hard (one label per pixel) or soft/fuzzy (membership or proportion of each class);
- Pixel-based or object-based.
The overall workflow is shown in Figure 8.1.
8.2 Unsupervised classification
\(k\)-means partitions \(n\) pixels into \(K\) clusters by minimising the within-cluster sum of squares,
\[ J=\sum_{k=1}^{K}\sum_{\mathbf{x}_i\in C_k}\lVert\mathbf{x}_i-\boldsymbol{\mu}_k\rVert^{2}, \tag{8.1}\]
by iterating assignment of pixels to the nearest centroid and recomputation of centroids until convergence. ISODATA extends this by splitting clusters with large variance and merging close clusters. The analyst then labels the clusters using reference data or interpretation. Unsupervised methods are valuable for exploring data and discovering the spectral classes present, but the resulting clusters need not correspond to meaningful information classes.
8.3 Supervised parametric classifiers
Training samples of known class define each class in feature space. The minimum-distance-to-means classifier assigns \(\mathbf{x}\) to the class with the nearest mean \(\boldsymbol{\mu}_k\) in Euclidean distance. It ignores variance and correlation, so it performs poorly for classes with different spreads.
The Mahalanobis distance accounts for the covariance \(\boldsymbol{\Sigma}_k\) of class \(k\): \[ d_M^{2}(\mathbf{x},k)=(\mathbf{x}-\boldsymbol{\mu}_k)^{\top}\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}-\boldsymbol{\mu}_k). \tag{8.2}\]
The Gaussian maximum-likelihood classifier (MLC) is derived from Bayes’ rule. Assigning \(\mathbf{x}\) to the class with the highest posterior probability \(P(\omega_k\mid\mathbf{x})\propto p(\mathbf{x}\mid\omega_k)P(\omega_k)\), and assuming multivariate normal class-conditional densities, the decision rule is to choose the class maximising the discriminant function
\[ g_k(\mathbf{x})=\ln P(\omega_k)-\tfrac12\ln\lvert\boldsymbol{\Sigma}_k\rvert-\tfrac12(\mathbf{x}-\boldsymbol{\mu}_k)^{\top}\boldsymbol{\Sigma}_k^{-1}(\mathbf{x}-\boldsymbol{\mu}_k). \tag{8.3}\]
MLC is optimal when its assumptions hold: normal class distributions, accurate estimates of \(\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k\), and correct priors. Estimating \(\boldsymbol{\Sigma}_k\) requires at least \(B+1\) (preferably \(10B\) to \(30B\)) samples per class; with too few, the covariance is singular or unstable. This is the origin of the Hughes phenomenon: for a fixed training-set size, accuracy first improves and then declines as the number of features rises.
8.4 Non-parametric and machine-learning classifiers
Decision trees recursively split feature space with axis-aligned thresholds chosen to maximise the reduction in impurity, for example the Gini impurity \(G=1-\sum_k p_k^{2}\). A random forest grows many trees, each on a bootstrap sample of the training data and using a random subset of features at each split, and assigns the majority vote. Random forests are robust to noise, handle many correlated features, need little tuning, and provide out-of-bag error and variable-importance estimates. They do not extrapolate beyond the training feature space and can inherit biases in the training data.
Support vector machines (SVMs) find the hyperplane \(\mathbf{w}^{\top}\mathbf{x}+b=0\) with maximum margin between two classes. For non-separable data the soft-margin problem is
\[ \min_{\mathbf{w},b,\xi}\ \tfrac12\lVert\mathbf{w}\rVert^{2}+C\sum_{i}\xi_i\quad\text{s.t.}\quad y_i(\mathbf{w}^{\top}\mathbf{x}_i+b)\ge1-\xi_i,\ \xi_i\ge0, \tag{8.4}\]
where \(C\) controls the penalty of violations. Non-linear boundaries are obtained with a kernel, most commonly the radial basis function \(K(\mathbf{x},\mathbf{x}')=\exp(-\gamma\lVert\mathbf{x}-\mathbf{x}'\rVert^{2})\). SVMs perform well with small training sets in high-dimensional spaces, but require careful tuning of \(C\) and \(\gamma\).
Deep learning learns features directly from image patches. Convolutional neural networks (CNNs) apply learned convolution kernels (the operation of Equation 7.3 with trainable weights \(w\)) followed by non-linearities and pooling; semantic-segmentation architectures such as encoder–decoder networks assign a class to every pixel and exploit spatial context that pixel-based classifiers ignore. Deep models can achieve excellent accuracy given sufficient labelled data, but require many training samples, are prone to overfit to the region and season of the training data, and are difficult to interpret. Transfer learning and large pre-trained foundation models for Earth observation reduce the labelled-data requirement.
8.5 Object-based image analysis (OBIA)
OBIA first partitions the image into homogeneous segments (objects) and then classifies the objects using spectral, geometric, textural, and contextual features. Multiresolution segmentation merges neighbouring pixels while minimising a heterogeneity criterion of the form
\[ f=w_{\text{color}}\,h_{\text{color}}+(1-w_{\text{color}})\,h_{\text{shape}}, \tag{8.5}\]
subject to a scale parameter that limits the growth of the objects. OBIA suppresses the salt-and-pepper effect of pixel-based classification, is well suited to very-high-resolution imagery where a single object (a tree crown, a building) comprises many pixels, and allows rules based on geometry and context (“a water object adjacent to a river network”). It is sensitive to the segmentation scale, which should be tuned to the target object size.
8.6 Sub-pixel classification and spectral unmixing
In coarse imagery, many pixels are mixtures. The linear spectral mixture model represents the observed spectrum as a weighted sum of \(K\) pure endmember spectra \(\mathbf{e}_k\):
\[ \mathbf{x}=\sum_{k=1}^{K}f_k\,\mathbf{e}_k+\boldsymbol{\varepsilon}=\mathbf{E}\mathbf{f}+\boldsymbol{\varepsilon},\qquad \sum_{k}f_k=1,\quad f_k\ge0, \tag{8.6}\]
where \(f_k\) are fractional abundances. The unconstrained least-squares solution is \(\hat{\mathbf{f}}=(\mathbf{E}^{\top}\mathbf{E})^{-1}\mathbf{E}^{\top}\mathbf{x}\); the constraints require a constrained (fully constrained least squares) solver. The residual \(\lVert\boldsymbol{\varepsilon}\rVert\) indicates whether the endmember set explains the pixel. The model is exact only for checkerboard-type (areal) mixing; non-linear models are needed where multiple scattering occurs (for example, in vegetation canopies over soil).
Example 8.1 (Two-endmember unmixing) Take two bands, endmembers \(\mathbf{e}_{\text{veg}}=(0.05,\,0.45)^{\top}\) and \(\mathbf{e}_{\text{soil}}=(0.25,\,0.25)^{\top}\), and a pixel \(\mathbf{x}=(0.10,\,0.40)^{\top}\). With \(f_{\text{soil}}=1-f_{\text{veg}}\), band 1 gives \(0.05f+0.25(1-f)=0.10\Rightarrow f_{\text{veg}}=0.75\). Band 2 confirms: \(0.45(0.75)+0.25(0.25)=0.40\). The pixel is 75 % vegetation and 25 % soil.
The spectral angle mapper (SAM) compares the direction (not the magnitude) of two spectra, \[ \alpha=\arccos\!\left(\frac{\mathbf{x}\cdot\mathbf{y}}{\lVert\mathbf{x}\rVert\,\lVert\mathbf{y}\rVert}\right), \tag{8.7}\] and is therefore insensitive to a multiplicative change in illumination.
8.7 Accuracy assessment
A classification is a model, and its map must be accompanied by a quantified statement of its reliability, computed from reference data of higher quality than the map and independent of the training data.
8.7.1 Sampling design
Common probability designs are simple random, systematic, and stratified random sampling, the last usually with strata defined by mapped classes so that rare classes receive enough samples. Cluster sampling reduces field costs but complicates variance estimation. A rule of thumb for a binomial proportion is \[ n=\frac{z^{2}\,p(1-p)}{e^{2}}, \tag{8.8}\] where \(p\) is the expected accuracy, \(e\) the allowable half-width of the confidence interval, and \(z\) the standard normal quantile (1.96 for 95 %). Because probability sampling underpins valid inference, samples chosen “where it is convenient” or by the same analyst who chose the training sites must not be used.
8.7.2 The error matrix
The error (confusion) matrix cross-tabulates map labels (rows) against reference labels (columns). Table 8.1 gives an example with \(n=300\) samples.
| Map Reference | Forest | Cropland | Urban | Water | Row total | User’s accuracy |
|---|---|---|---|---|---|---|
| Forest | 85 | 8 | 2 | 0 | 95 | 89.5 % |
| Cropland | 10 | 78 | 5 | 2 | 95 | 82.1 % |
| Urban | 3 | 7 | 60 | 0 | 70 | 85.7 % |
| Water | 2 | 2 | 3 | 33 | 40 | 82.5 % |
| Column total | 100 | 95 | 70 | 35 | 300 | |
| Producer’s accuracy | 85.0 % | 82.1 % | 85.7 % | 94.3 % |
Let \(n_{ij}\) be the entry in row \(i\) and column \(j\), \(n_{i+}\) the row total, \(n_{+j}\) the column total, and \(n\) the grand total. The standard measures are:
- Overall accuracy: \(\mathrm{OA}=\dfrac{1}{n}\sum_i n_{ii}\).
- User’s accuracy (precision, \(1-\) commission error): \(\mathrm{UA}_i=\dfrac{n_{ii}}{n_{i+}}\): the probability that a pixel labelled \(i\) on the map really is \(i\).
- Producer’s accuracy (recall, \(1-\) omission error): \(\mathrm{PA}_j=\dfrac{n_{jj}}{n_{+j}}\): the probability that a reference \(j\) location is mapped as \(j\).
- F1 score: \(F_1=\dfrac{2\,\mathrm{UA}\cdot\mathrm{PA}}{\mathrm{UA}+\mathrm{PA}}\).
- Cohen’s kappa, the agreement corrected for chance:
\[ \kappa=\frac{p_o-p_e}{1-p_e},\qquad p_o=\frac{1}{n}\sum_i n_{ii},\quad p_e=\frac{1}{n^{2}}\sum_i n_{i+}\,n_{+i}. \tag{8.9}\]
Example 8.2 (Accuracy measures from Table 8.1) The diagonal sums to \(85+78+60+33=256\), so \(\mathrm{OA}=256/300=0.853\). Chance agreement is \(p_e=(95\cdot100+95\cdot95+70\cdot70+40\cdot35)/300^{2}=24\,825/90\,000=0.276\), hence \(\kappa=(0.853-0.276)/(1-0.276)\approx0.797\). Water has high producer’s accuracy (94.3 %) but only 82.5 % user’s accuracy: nearly all real water is found, but some mapped water is not water.
Kappa has been widely criticised: it can be difficult to interpret, depends on the marginals, and provides little information beyond overall accuracy. Best practice is to report the full error matrix, class-wise user’s and producer’s accuracies with confidence intervals, and, optionally, decompositions of disagreement into quantity and allocation components.
8.7.3 Area estimation and bias correction
Maps are not unbiased estimators of class area: pixel counting confounds omission and commission errors. With a stratified sample in which map class \(i\) (proportion of mapped area \(W_i\)) contains \(n_{i+}\) samples, the estimated proportion of the area that is truly class \(j\) is
\[ \hat{p}_{j}=\sum_{i}W_i\,\frac{n_{ij}}{n_{i+}},\qquad \hat{A}_j=A_{\mathrm{tot}}\,\hat{p}_j, \tag{8.10}\]
with an accompanying standard error derived from the stratified-estimator variance. Reporting the adjusted area with its confidence interval is essential, for example, in national forest-loss statistics, where the rare class (change) makes the map-based area unreliable.
8.7.4 Common pitfalls
- Spatial autocorrelation and leakage. Neighbouring training and test pixels are correlated; random pixel-level cross-validation then greatly overstates accuracy. Use spatially blocked or independent-site validation.
- Class imbalance. Overall accuracy is dominated by common classes; rare classes must be evaluated separately.
- Reference-data error. Interpretation of high-resolution imagery, field data, and crowd-sourced labels have their own errors; when substantial, they bias accuracy and area estimates.
- Temporal mismatch between imagery and reference data.
- Extrapolation. A model validated in one region and season is not necessarily valid in another.
Exercises
Exercise 8.1 Class means are \(\boldsymbol{\mu}_1=(0.08,\,0.35)^{\top}\) and \(\boldsymbol{\mu}_2=(0.20,\,0.25)^{\top}\) in (red, NIR). Assign the pixel \(\mathbf{x}=(0.10,\,0.30)^{\top}\) using the minimum-distance classifier. Give an example of when this classifier fails but the MLC succeeds.
Exercise 8.2 A two-class error matrix (rows: map; columns: reference) is \(\begin{pmatrix}45&5\\10&40\end{pmatrix}\). Compute the overall accuracy, both user’s and producer’s accuracies, the \(F_1\) score for class 1, and kappa.
Exercise 8.3 How many samples are needed to estimate overall accuracy with a 95 % confidence half-width of 5 %, if the expected accuracy is 85 %? What if the half-width is reduced to 2.5 %? Discuss why an allocation that is proportional to mapped area is a poor design for a rare class.
Exercise 8.4 Two spectra are \(\mathbf{x}=(0.1,0.2,0.4)\) and \(\mathbf{y}=(0.2,0.4,0.8)\); a third is \(\mathbf{z}=(0.3,0.2,0.1)\). Compute the spectral angle between \(\mathbf{x}\) and \(\mathbf{y}\), and between \(\mathbf{x}\) and \(\mathbf{z}\), and interpret the results.
Exercise 8.5 A pixel has spectrum \((0.12,\,0.36)\) in two bands. The endmembers are vegetation \((0.05,\,0.45)\) and soil \((0.25,\,0.25)\). Solve for the fractions using the sum-to-one constraint and the first band only, then check the second band. What does the discrepancy imply?
Exercise 8.6 Conceptual. A study reports 97 % accuracy for a random-forest crop classifier, validated by a random 70/30 split of pixels sampled from 20 fields. Explain the source of bias and propose a validation design that would give a defensible estimate.
9 Active Sensing, LiDAR, and Imaging Spectroscopy
9.1 Imaging radar
Radar (radio detection and ranging) transmits microwave pulses and records the echo. Its all-weather and day-and-night capability, sensitivity to structure and moisture, and phase information make it a distinct source of information from optical sensors.
9.1.1 Frequency bands
| Band | Wavelength | Frequency | Typical penetration and sensitivity |
|---|---|---|---|
| X | 2.4–3.8 cm | 8–12 GHz | canopy top, fine surface roughness |
| C | 3.8–7.5 cm | 4–8 GHz | crops, leaves, soil roughness, ocean |
| S | 7.5–15 cm | 2–4 GHz | intermediate |
| L | 15–30 cm | 1–2 GHz | branches, trunks, soil moisture, deeper into canopy |
| P | 30–100 cm | 0.3–1 GHz | trunks, biomass, sub-surface in very dry soils |
9.1.2 The radar equation and backscatter
The power received from a point target at range \(R\) is given by the radar equation:
\[ P_r=\frac{P_t\,G^{2}\,\lambda^{2}\,\sigma}{(4\pi)^{3}R^{4}}, \tag{9.1}\]
where \(P_t\) is the transmitted power, \(G\) the antenna gain, and \(\sigma\) the target’s radar cross-section. For distributed targets such as terrain, the geometry-independent quantity is the backscattering coefficient \(\sigma^{0}=\sigma/A\), the cross-section per unit ground area, usually expressed in decibels:
\[ \sigma^{0}_{\mathrm{dB}}=10\log_{10}\sigma^{0}. \tag{9.2}\]
Backscatter depends on the radar parameters (wavelength, polarisation, incidence angle) and target parameters (roughness relative to \(\lambda\), dielectric constant, and geometry). Three main scattering mechanisms are recognised: surface scattering (bare soil, water), volume scattering (vegetation canopies, dry snow), and double-bounce (ground–trunk or ground–wall corner reflections, strong in flooded forests and urban areas).
9.1.3 Geometry and resolution
Radar is a side-looking system (Figure 9.1). It measures the time delay of echoes, so targets are ordered by slant range \(R\), not by ground position.
The slant-range resolution depends on the signal bandwidth \(B\) (not the pulse length after pulse compression), and projects onto the ground with the incidence angle \(\theta\):
\[ \rho_{sr}=\frac{c}{2B},\qquad \rho_{gr}=\frac{c}{2B\sin\theta}. \tag{9.3}\]
Resolution in azimuth for a real-aperture radar is \(\lambda R/L\), which degrades with range and is impractical from orbit. Synthetic aperture radar (SAR) synthesises a long antenna by coherently combining echoes collected along the flight path; the result, remarkably, is an azimuth resolution of about \(L/2\) (half the physical antenna length \(L\)), independent of range and wavelength.
9.1.4 Geometric distortions and speckle
Because of the side-looking geometry, steep terrain produces foreshortening (slopes facing the sensor are compressed), layover (tops of tall features displaced toward the sensor, overlapping the base), and shadow (areas hidden from the illumination, with no signal). These effects mandate specific processing (radiometric terrain correction with a DEM, and geometry-aware masks) and careful selection of the orbit (ascending vs. descending) in mountainous areas.
Speckle is a granular pattern arising from the coherent interference of many scatterers within a resolution cell. It is multiplicative, \(I=\sigma^{0}\cdot n\), with \(n\) having unit mean and a variance that decreases with the number of independent looks \(L_{\text{looks}}\) as \(1/L_{\text{looks}}\). Multilooking averages looks at the expense of resolution; adaptive filters (Lee, refined Lee, Gamma-MAP) and multitemporal filters reduce speckle while preserving edges. Speckle filtering is not “noise removal” in the optical sense and should be applied deliberately.
9.1.5 Polarimetry
Radar systems transmit and receive in horizontal (H) or vertical (V) polarisation, producing HH, VV, and cross-polarised (HV, VH) channels. Cross-polarisation arises mainly from volume scattering and is sensitive to vegetation structure and biomass, while co-polarised ratios are sensitive to surface roughness and moisture. Fully polarimetric data allow target decompositions (Freeman–Durden, Cloude–Pottier entropy/anisotropy/alpha) that separate the scattering mechanisms.
9.2 Interferometric SAR (InSAR)
A SAR image is complex-valued: it records the amplitude and the phase of the echo. The phase is proportional to the two-way range in wavelengths:
\[ \phi=-\frac{4\pi}{\lambda}R+\phi_{\text{scatt}}. \tag{9.4}\]
An interferogram is formed by multiplying one image by the complex conjugate of another acquired from a slightly different position or time, so that the scattering phase cancels and the phase difference remains:
\[ \Delta\phi=\underbrace{\phi_{\text{flat}}}_{\text{reference surface}}+\underbrace{\phi_{\text{topo}}}_{\text{height}}+\underbrace{\phi_{\text{defo}}}_{\text{motion}}+\underbrace{\phi_{\text{atm}}}_{\text{delay}}+\phi_{\text{noise}}. \tag{9.5}\]
After removing the flat-earth and topographic terms (using a DEM), the residual phase measures line-of-sight (LOS) displacement:
\[ d_{\text{LOS}}=-\frac{\lambda}{4\pi}\,\Delta\phi_{\text{defo}}. \tag{9.6}\]
One full fringe (\(2\pi\)) thus corresponds to \(\lambda/2\) of LOS motion (about 2.8 cm for C-band). Quality is measured by the coherence \(\gamma\in[0,1]\), the normalised complex correlation between the two images; it is lost through temporal decorrelation (changes in the surface, especially vegetation), geometric (baseline) decorrelation, and noise. Time-series techniques (persistent scatterer InSAR, small-baseline subset) use many interferograms to separate deformation from atmospheric artefacts, reaching millimetre-per-year precision for subsidence, landslides, volcanic inflation, and fault motion. Phase is measured modulo \(2\pi\); phase unwrapping recovers the absolute phase and is a major source of error.
9.3 Light detection and ranging (LiDAR)
LiDAR measures the two-way travel time \(t\) of a laser pulse to derive the range:
\[ R=\frac{c\,t}{2}. \tag{9.7}\]
A discrete-return system records one to several returns per pulse (first, last, intermediate); a full-waveform system digitises the entire backscattered intensity profile, which characterises the vertical structure of vegetation. Airborne laser scanning (ALS) combines the laser with a GNSS receiver and an inertial measurement unit (IMU); the direct georeferencing equation gives the ground position of each return,
\[ \mathbf{x}_g=\mathbf{x}_{\mathrm{GNSS}}+\mathbf{R}_{\mathrm{INS}}\left(\mathbf{R}_{\mathrm{bore}}\,\mathbf{r}_{s}+\mathbf{l}\right), \tag{9.8}\]
where \(\mathbf{r}_s\) is the range vector in the scanner frame (from \(R\) and the scan angle), \(\mathbf{R}_{\mathrm{bore}}\) the boresight rotation, \(\mathbf{l}\) the lever-arm offset between the GNSS antenna and the scanner, and \(\mathbf{R}_{\mathrm{INS}}\) the attitude rotation. Errors in position, attitude, boresight, and range all propagate to the point cloud.
The result is a point cloud. Its processing includes ground filtering (progressive TIN densification, cloth simulation), interpolation of ground points to a DTM, and of first returns to a DSM. The canopy height model (CHM) is
\[ \mathrm{CHM}=\mathrm{DSM}-\mathrm{DTM}. \tag{9.9}\]
Point density (points per m²) controls the detail retrievable; typical airborne surveys of 5–20 pts/m² support DTMs with vertical RMSE of 5–15 cm on open ground. Forest metrics (height percentiles, canopy cover, gap fraction, leaf-area density) are computed from the vertical distribution of returns and calibrated to field plots to estimate biomass and carbon. Terrestrial laser scanning gives centimetre-scale trees and structures, and spaceborne LiDAR missions sample forest height and ice-sheet and vegetation profiles as sparse footprints or tracks rather than as images.
9.4 Imaging spectroscopy (hyperspectral remote sensing)
Imaging spectrometers acquire hundreds of narrow, contiguous bands, providing a nearly continuous spectrum for every pixel. Diagnostic absorption features, which are invisible in broad multispectral bands, allow the identification of minerals, leaf biochemicals (chlorophyll, water, nitrogen, lignin), and water constituents.
Continuum removal normalises a spectrum by a convex hull \(C(\lambda)\) fitted across the reflectance peaks, isolating the absorption feature:
\[ \mathrm{CR}(\lambda)=\frac{\rho(\lambda)}{C(\lambda)},\qquad D=1-\mathrm{CR}(\lambda_{\min}), \tag{9.10}\]
with \(D\) the band depth, which is related to the abundance of the absorbing material, and the position \(\lambda_{\min}\) diagnostic of its identity.
Processing includes atmospheric correction with water-vapour retrieval, noise-aware dimensionality reduction (MNF), endmember extraction (pixel purity index, N-FINDR, vertex component analysis), and identification by spectral matching (SAM, spectral-feature fitting) against spectral libraries. The main challenges are the curse of dimensionality (Chapter 8), high data volume, low SNR in narrow bands, and the difficulty of moving from lab spectra to pixels affected by mixing, scale, and illumination.
9.5 UAV photogrammetry and structure-from-motion
Small unmanned aircraft carrying consumer-grade or multispectral cameras have opened centimetre-scale mapping. Structure-from-motion (SfM) photogrammetry reconstructs three-dimensional structure and camera poses from overlapping photographs by automatically matching features across images and solving a bundle adjustment. Standard products are a dense point cloud, a DSM, and an orthomosaic. The GSD follows Equation 5.3, \(\mathrm{GSD}=H_o\,p/f\). Typical designs use 75–80 % forward and 65–70 % side overlap, well-distributed ground control points (or RTK/PPK positioning), and a standardised radiometric calibration (reflectance panels, irradiance sensors) for multispectral work. SfM works poorly on homogeneous or moving surfaces (water, dense uniform crops, tree canopy in wind).
Exercises
Exercise 9.1 A SAR has a chirp bandwidth \(B=150\) MHz and a 10 m antenna. (a) Compute the slant-range resolution and the ground-range resolution at \(\theta=30^{\circ}\). (b) What is the best azimuth resolution attainable? (c) Explain why increasing the antenna length worsens azimuth resolution in SAR, contrary to real-aperture radar.
Exercise 9.2 An interferogram at C-band (\(\lambda=5.6\) cm) shows 3.5 fringes of deformation phase across a subsidence bowl. Compute the maximum LOS displacement. Explain why the measurement is only of the LOS component and how using ascending and descending orbits helps.
Exercise 9.3
- A laser pulse returns 6.67 µs after emission. Compute the range. (b) An airborne DSM has an elevation of 1427.6 m at a tree and the DTM at the same location has 1409.9 m. Compute the canopy height. (c) Name three sources of error in the resulting CHM.
Exercise 9.4 Convert \(\sigma^{0}=0.05\) to decibels, and \(-20\) dB to a linear value. A flooded field changes from \(-8\) dB to \(-19\) dB. By what factor does the linear backscatter change? Why does calm open water appear dark in SAR imagery?
Exercise 9.5 Explain, using the concept of spectral resolution, why two minerals with similar broadband spectra may be discriminated by a hyperspectral sensor but not by a multispectral one. What is the role of continuum removal?
Exercise 9.6 Design problem. Choose between optical and SAR data (or a fusion) for mapping seasonal flooding in a persistently cloudy river delta, and separately for estimating forest biomass in a tropical forest. Discuss band/frequency, polarisation, temporal sampling, and limitations.
10 Spatial Analysis and Geostatistics
Spatial analysis is the collection of techniques for describing, explaining, and predicting phenomena by exploiting their locations. Its defining feature is that location matters: the results change if the locations are permuted.
10.1 Map algebra and overlay
Tomlin’s map algebra provides a small, closed set of raster operations that can be composed into arbitrarily complex models. They are classified by the extent of the neighbourhood that contributes to each output cell (Figure 10.1).
- Local operations compute \(O(i,j)=f\bigl(A(i,j),B(i,j),\dots\bigr)\): arithmetic, reclassification, Boolean logic, and index calculation (for example, NDVI is a local operation).
- Focal operations compute an output from a neighbourhood: smoothing, slope, density, majority filters. Convolution (Equation 7.3) is a focal operation.
- Zonal operations summarise a value layer within each zone of a categorical layer: mean rainfall by district, land-cover area by watershed.
- Global operations depend on every cell: Euclidean and cost distance, flow accumulation, viewsheds.
Overlay combines layers. In Boolean overlay, criteria are combined with AND/OR/NOT to select suitable areas. In weighted linear combination (weighted overlay), standardised criterion scores \(s_i\in[0,1]\) or \([0,10]\) are combined with weights \(w_i\),
\[ S=\sum_{i=1}^{n}w_i\,s_i,\qquad \sum_i w_i=1, \tag{10.1}\]
which permits trade-offs among criteria and is the core of multi-criteria decision analysis (MCDA). Weights are often derived from pairwise comparisons in the Analytic Hierarchy Process (AHP): the principal eigenvector of the comparison matrix gives the weights, and the consistency ratio \(\mathrm{CR}=\mathrm{CI}/\mathrm{RI}\), with \(\mathrm{CI}=(\lambda_{\max}-n)/(n-1)\) and \(\mathrm{RI}\) a tabulated random index (0.58, 0.90, 1.12 for \(n=3,4,5\)), should be below about 0.10. Vector overlay operations (intersect, union, clip, erase, dissolve) and buffering (the set of points within distance \(d\) of a feature) are the vector analogues.
10.2 Proximity, cost, and networks
The Euclidean distance to the nearest source is a global operation; its generalisation, cost distance, accumulates a friction (travel time, slope, land-cover resistance) as the path crosses cells, and the least-cost path between two locations is found by Dijkstra-type algorithms. In network analysis the street network is a graph \(G=(V,E)\) with edge impedances; the shortest path from an origin minimises total impedance; a service area is the set of network locations reachable within a threshold; and location–allocation problems such as the \(p\)-median problem choose \(p\) facility sites to minimise total weighted travel,
\[ \min\ \sum_{i}\sum_{j}w_i\,d_{ij}\,x_{ij}\quad\text{s.t.}\quad\sum_{j}x_{ij}=1\ \forall i,\ \ x_{ij}\le y_j,\ \ \sum_j y_j=p,\ \ x_{ij},y_j\in\{0,1\}, \tag{10.2}\]
where \(w_i\) is the demand at \(i\), \(d_{ij}\) the network distance, \(x_{ij}=1\) if demand \(i\) is assigned to facility \(j\), and \(y_j=1\) if a facility is located at \(j\). Such models underpin accessibility studies for schools, water points, and health services, and their conclusions are sensitive to the network’s completeness and to the travel-speed assumptions.
10.3 Spatial autocorrelation
Given values \(z_i\) at \(n\) locations and a spatial weights matrix \(\mathbf{W}=(w_{ij})\) encoding neighbourhood (contiguity, \(k\) nearest neighbours, distance band, inverse distance), Moran’s \(I\) measures global spatial autocorrelation:
\[ I=\frac{n}{S_0}\;\frac{\sum_{i}\sum_{j}w_{ij}(z_i-\bar z)(z_j-\bar z)}{\sum_{i}(z_i-\bar z)^{2}},\qquad S_0=\sum_{i}\sum_{j}w_{ij}. \tag{10.3}\]
Under spatial randomness \(E[I]=-1/(n-1)\approx0\). Values significantly above the expectation indicate clustering of similar values, values below it a checkerboard-like dispersion. Significance is tested analytically (normal or randomisation assumption) or by permutation. Moran’s \(I\) is a global measure and describes the overall tendency only.
Local indicators of spatial association (LISA) decompose the global statistic into contributions of each location; the local Moran’s \(I_i=\frac{z_i-\bar z}{m_2}\sum_j w_{ij}(z_j-\bar z)\) identifies clusters (high–high, low–low) and spatial outliers (high–low, low–high). The Getis–Ord \(G_i^{*}\) statistic identifies hot spots and cold spots by comparing the sum of values in a neighbourhood with its expectation. With many simultaneous local tests, multiple-testing correction (false-discovery-rate control) is required to avoid spurious hot spots.
10.4 Point-pattern analysis
A point pattern is a set of event locations (disease cases, tree stems, fires) within a study region. The first-order property is the intensity \(\lambda(\mathbf{s})\), the expected number of events per unit area; the second-order properties describe dependence between events. The reference model is complete spatial randomness (CSR), a homogeneous Poisson process.
For CSR with intensity \(\lambda\), the expected mean nearest-neighbour distance is \(\bar d_{\mathrm{exp}}=1/(2\sqrt{\lambda})\); the nearest-neighbour ratio \(R=\bar d_{\mathrm{obs}}/\bar d_{\mathrm{exp}}\) is below 1 for clustering and above 1 for regularity. A scale-dependent description is provided by Ripley’s \(K\)-function,
\[ K(r)=\frac{1}{\lambda}\,E\bigl[\text{number of events within distance }r\text{ of a typical event}\bigr],\qquad L(r)=\sqrt{K(r)/\pi}-r, \tag{10.4}\]
for which CSR gives \(K(r)=\pi r^{2}\) and \(L(r)=0\); \(L(r)>0\) indicates clustering at scale \(r\). Kernel density estimation smooths events into a continuous intensity surface,
\[ \hat\lambda(\mathbf{s})=\frac{1}{h^{2}}\sum_{i=1}^{n}K\!\left(\frac{\lVert\mathbf{s}-\mathbf{s}_i\rVert}{h}\right), \tag{10.5}\]
where the bandwidth \(h\) controls the trade-off between bias (oversmoothing) and variance and must be reported. Because event maps merely reproduce the population distribution, epidemiological analyses compare case and control (or population-at-risk) densities, for example through the ratio of kernel densities, or use the spatial scan statistic, which searches for windows whose likelihood ratio of inside versus outside risk is greatest.
10.5 Geostatistics and spatial interpolation
Geostatistics treats a spatially continuous variable as a realisation of a random field and estimates values at unsampled locations from those observed. Under the intrinsic hypothesis, the expected differences are zero and their variance depends only on the separation vector, so that the semivariogram may be estimated from data:
\[ \hat\gamma(h)=\frac{1}{2N(h)}\sum_{i=1}^{N(h)}\bigl[z(\mathbf{x}_i)-z(\mathbf{x}_i+\mathbf{h})\bigr]^{2}, \tag{10.6}\]
where \(N(h)\) is the number of pairs separated by lag \(h\). A theoretical model, such as the spherical model,
\[ \gamma(h)=\begin{cases}c_0+c\left[\dfrac{3h}{2a}-\dfrac{1}{2}\left(\dfrac{h}{a}\right)^{3}\right],&0<h<a,\\[2ex] c_0+c,&h\ge a,\end{cases} \tag{10.7}\]
is fitted to the empirical values (Figure 10.2). Its parameters have physical meaning: the nugget \(c_0\) (microscale variation and measurement error), the sill \(c_0+c\) (the total variance), and the range \(a\) (the distance beyond which observations are uncorrelated). Exponential and Gaussian models are common alternatives; anisotropy (directional dependence) is examined with directional variograms.
Ordinary kriging predicts \(\hat z(\mathbf{x}_0)=\sum_i\lambda_i z(\mathbf{x}_i)\) with weights that minimise the estimation variance subject to unbiasedness (\(\sum_i\lambda_i=1\)). The weights solve the kriging system
\[ \sum_{j=1}^{n}\lambda_j\,\gamma(\mathbf{x}_i,\mathbf{x}_j)+\mu=\gamma(\mathbf{x}_i,\mathbf{x}_0),\quad i=1,\dots,n,\qquad \sum_{j}\lambda_j=1, \tag{10.8}\]
where \(\mu\) is a Lagrange multiplier, and the kriging variance is \(\sigma_K^{2}(\mathbf{x}_0)=\sum_i\lambda_i\gamma(\mathbf{x}_i,\mathbf{x}_0)+\mu\). Kriging is the best linear unbiased predictor under its assumptions and, unlike other interpolators, gives a measure of uncertainty for each prediction. Variants include universal kriging and regression kriging (a regression on covariates plus kriged residuals), co-kriging, block kriging, and indicator kriging for probabilities.
Simpler interpolators are inverse-distance weighting (IDW),
\[ \hat z(\mathbf{x}_0)=\frac{\sum_i w_i\,z_i}{\sum_i w_i},\qquad w_i=d_i^{-p}, \tag{10.9}\]
which is exact at the data, has no uncertainty estimate, and yields “bull’s-eyes” around samples; splines minimise curvature. All interpolators should be evaluated by cross-validation (leave-one-out or spatially blocked), with reporting of the mean error, RMSE, and calibration of the predicted uncertainty.
10.6 Terrain analysis
From a DEM, local geometry is derived with a \(3\times3\) window. Label the elevations
\[ \begin{pmatrix}a&b&c\\d&e&f\\g&h&i\end{pmatrix} \]
with the top row toward north. Horn’s method estimates the gradient components with weighted finite differences over a cell size \(\Delta\),
\[ \frac{\partial z}{\partial x}=\frac{(c+2f+i)-(a+2d+g)}{8\Delta},\qquad \frac{\partial z}{\partial y}=\frac{(a+2b+c)-(g+2h+i)}{8\Delta}, \tag{10.10}\]
from which the slope and aspect (downslope direction, clockwise from north) are
\[ S=\arctan\sqrt{\left(\frac{\partial z}{\partial x}\right)^{2}+\left(\frac{\partial z}{\partial y}\right)^{2}},\qquad A=\left[\operatorname{atan2}\!\left(\frac{\partial z}{\partial x},\frac{\partial z}{\partial y}\right)+180^{\circ}\right]\bmod 360^{\circ}. \tag{10.11}\]
Curvature (profile and plan) describes the acceleration and convergence of flow. In hydrological analysis, the D8 algorithm routes water from each cell to the steepest-descending of its eight neighbours; the count of upstream cells (flow accumulation) delineates streams and watersheds after depressions (“sinks”) are filled or breached. The topographic wetness index \(\mathrm{TWI}=\ln\bigl(a_s/\tan\beta\bigr)\), where \(a_s\) is the specific catchment area and \(\beta\) the slope, indicates where soil is likely to be saturated. Viewsheds, terrain-ruggedness indices, and solar radiation models are further global or focal derivatives. The DEM resolution and vertical accuracy govern the reliability of every one of these derivatives.
10.7 Spatial regression and modelling
Ordinary least squares assumes independent errors. When residuals are spatially autocorrelated, OLS coefficient estimates remain unbiased but their standard errors are underestimated, inflating significance. Diagnose residual autocorrelation (Moran’s \(I\) of the residuals) before trusting a model. Two classical specifications are:
\[ \text{Spatial lag: }\ \mathbf{y}=\rho\,\mathbf{W}\mathbf{y}+\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon},\qquad \text{Spatial error: }\ \mathbf{y}=\mathbf{X}\boldsymbol{\beta}+\mathbf{u},\ \ \mathbf{u}=\lambda_{e}\mathbf{W}\mathbf{u}+\boldsymbol{\varepsilon}. \tag{10.12}\]
The lag model represents substantive spillover (the outcome in one place depends on neighbours’ outcomes), the error model represents autocorrelated omitted influences. Geographically weighted regression (GWR) addresses non-stationarity by fitting a separate weighted regression at each location \((u_i,v_i)\):
\[ \hat{\boldsymbol\beta}(u_i,v_i)=\bigl(\mathbf{X}^{\top}\mathbf{W}_i\mathbf{X}\bigr)^{-1}\mathbf{X}^{\top}\mathbf{W}_i\,\mathbf{y},\qquad w_{ij}=\exp\!\left[-\tfrac12\left(\frac{d_{ij}}{b}\right)^{2}\right], \tag{10.13}\]
with bandwidth \(b\) chosen by cross-validation or information criteria; coefficient surfaces should be interpreted with caution because of local multicollinearity and multiple testing. For small-area disease mapping, hierarchical Bayesian models such as the Besag–York–Mollié model, \(\log\theta_i=\alpha+u_i+v_i\) with a spatially structured effect \(u_i\) (conditional autoregressive) and an unstructured effect \(v_i\), borrow strength from neighbours to stabilise noisy relative risks in areas with small populations.
10.8 Pitfalls of spatial analysis
- The modifiable areal unit problem and the ecological fallacy: inferences about individuals cannot be made from area-level associations.
- Edge (boundary) effects: features outside the study region influence statistics for those inside; neighbourhood-based statistics are biased near edges.
- Non-stationarity: a single global model may describe no place well.
- Scale mismatch: predictors, responses, and processes operate at different scales.
- Positional and attribute uncertainty propagates through every operation (Chapter 12).
Exercises
Exercise 10.1 Three samples have values 10, 20, and 30 at distances of 2, 4, and 5 km from a target location. Compute the IDW estimate for (a) \(p=2\) and (b) \(p=1\). What is the effect of increasing \(p\)?
Exercise 10.2 A \(3\times3\) DEM window with 30 m cells is \(\begin{pmatrix}100&102&104\\98&100&102\\96&98&100\end{pmatrix}\) (top row north). Compute the slope (in degrees) and the aspect by Horn’s method.
Exercise 10.3 Four cells lie in a row with values 1, 2, 3, 4 and rook (binary, adjacent) contiguity. Compute Moran’s \(I\) and its expectation under randomness. Is the pattern clustered?
Exercise 10.4 In a survey, five pairs of observations separated by approximately 100 m have squared-difference components with differences of 1, 2, 0, 3, and 2 units. Compute the empirical semivariance at that lag. What would happen to the estimate if only one pair were available?
Exercise 10.5 Twenty-five clinics are located in a 1 km² district. The mean observed nearest-neighbour distance is 80 m. Compute the expected value under CSR and the ratio \(R\), and interpret. Name two reasons a health-facility pattern would not be expected to be CSR.
Exercise 10.6 A site has criterion scores of 8, 4, and 6 (on a 0–10 scale) for slope, distance to road, and soil quality, with weights 0.5, 0.3, and 0.2. Compute the weighted suitability. Explain why a site with a score of 0 for a constraint (for example, a protected area) should be removed by a Boolean mask rather than by the weighted sum.
Exercise 10.7 Design problem. A ministry wants to allocate five mobile clinics to reduce travel time to care in a rural district. Specify the data you need, the analytical model (including the objective function), the sensitivity analyses, and the caveats of your conclusions.
11 Time Series and Change Detection
Change is often the phenomenon of interest: deforestation, urban growth, crop phenology, flood, drought, glacier retreat. Detecting change from imagery requires that real surface change be separated from apparent change due to illumination, atmosphere, phenology, sensor differences, and misregistration.
11.1 Prerequisites for reliable change detection
- Geometric co-registration to a fraction of a pixel (a registration error of one pixel can produce spurious change of the same order as the true signal in heterogeneous landscapes).
- Radiometric consistency: absolute atmospheric correction to surface reflectance, or relative normalisation using pseudo-invariant features (dark and bright, temporally stable targets such as deep water, bare rock, and pavement).
- Anniversary or phenology-matched dates, so that seasonal vegetation dynamics are not confused with land-cover change.
- Comparable sensors and viewing geometry, or explicit cross-calibration.
- A clear definition of change (thematic, magnitude, minimum mapping unit, and time span).
11.2 Bi-temporal methods
Image differencing subtracts the images pixelwise, \(\Delta=x_{t_2}-x_{t_1}\), and classifies as change those pixels whose difference departs from the mean by more than \(k\) standard deviations,
\[ \Delta<\mu_\Delta-k\,\sigma_\Delta\quad\text{or}\quad \Delta>\mu_\Delta+k\,\sigma_\Delta . \tag{11.1}\]
Image ratioing is used to reduce multiplicative effects (or, for SAR, because speckle is multiplicative). Change vector analysis (CVA) treats the pair of spectra as a vector in feature space and describes each change by a magnitude and direction,
\[ \lVert\Delta\rVert=\sqrt{\sum_{b=1}^{B}\bigl(x_{b,t_2}-x_{b,t_1}\bigr)^{2}}, \qquad \alpha=\arctan\!\left(\frac{\Delta_2}{\Delta_1}\right)\ \ \text{(for two features)}, \tag{11.2}\]
where magnitude indicates how much change occurred and direction indicates what type (for example, vegetation loss versus gain). Principal-component and multivariate alteration detection (MAD) methods find linear combinations of bands in which change stands out and are invariant to linear radiometric differences.
In post-classification comparison, each date is classified independently and the maps are cross-tabulated into a from–to transition matrix. This gives explicit thematic change information and does not require radiometric normalisation, but errors accumulate: if the two maps have independent accuracies of \(A_1\) and \(A_2\), the probability that both classifications are correct at a pixel is about \(A_1A_2\), so change maps are often much less accurate than either map. Errors of 5 % in each date can give change errors comparable to the true change rate. Hence direct, bi-temporal change classification, and stratified sampling that focuses on the change strata, are recommended when change is rare.
Example 11.1 (Burn severity) The differenced Normalised Burn Ratio is \(\mathrm{dNBR}=\mathrm{NBR}_{\text{pre}}-\mathrm{NBR}_{\text{post}}\). For a pixel with \(\mathrm{NBR}_{\text{pre}}=0.65\) and \(\mathrm{NBR}_{\text{post}}=-0.05\), \(\mathrm{dNBR}=0.70\), which falls in the high-severity class under the commonly used thresholds (about \(>0.66\)). Thresholds are empirically derived and should be calibrated to local field data.
11.3 Time-series analysis
With frequent acquisitions, the temporal profile of each pixel can be analysed as a signal. After cloud masking, gap filling (interpolation, or fitting smooth curves), and smoothing (Savitzky–Golay or Whittaker filters), the series can be decomposed into trend, seasonal, and remainder components. A common seasonal model is harmonic regression:
\[ y(t)=\beta_0+\beta_1t+\sum_{k=1}^{K}\bigl[a_k\cos(2\pi k t)+b_k\sin(2\pi k t)\bigr]+\varepsilon(t), \tag{11.3}\]
where \(t\) is time in years. For the first harmonic, the amplitude is \(\sqrt{a_1^{2}+b_1^{2}}\) and the phase \(\operatorname{atan2}(b_1,a_1)\) locates the seasonal peak within the year; these parameters serve as compact features for classification, since crops and natural vegetation differ in their seasonal signatures.
Phenology metrics extracted from a growing-season curve include start of season, end of season, length of season, peak value and timing, and integrated greenness; start and end of season are typically defined as the dates at which the smoothed curve crosses a fraction (for example 50 %) of its seasonal amplitude.
Breakpoint and disturbance detection algorithms test whether a model fitted to the series remains valid. BFAST decomposes a series and detects abrupt changes in trend and season; LandTrendr segments annual composites into linear pieces; CCDC (continuous change detection and classification) fits harmonic models and flags observations that depart from prediction by more than a multiple of the model’s RMSE for several consecutive observations. They can label disturbance and recovery dates at pixel level.
Monotonic trends are commonly assessed with the non-parametric Mann–Kendall test, whose statistic counts concordant and discordant pairs,
\[ S=\sum_{k=1}^{n-1}\sum_{j=k+1}^{n}\operatorname{sgn}(x_j-x_k), \tag{11.4}\]
with \(\mathrm{Var}(S)=n(n-1)(2n+5)/18\) in the absence of ties, and the trend magnitude by the Theil–Sen slope, the median of all pairwise slopes \((x_j-x_k)/(j-k)\). Both are robust to outliers and do not assume normality, but temporal autocorrelation (serial correlation) inflates significance and must be handled (pre-whitening).
11.4 Thresholding by Otsu’s method
Many mapping tasks (flooded versus non-flooded from SAR backscatter, burned versus unburned) reduce to selecting a threshold in a bimodal histogram. Otsu’s method chooses the threshold that maximises the between-class variance,
\[ \sigma_B^{2}=\omega_0\,\omega_1\,(\mu_0-\mu_1)^{2}, \tag{11.5}\]
where \(\omega_0,\omega_1\) are the class proportions and \(\mu_0,\mu_1\) the class means. It works well when the two classes have comparable sizes and fails when one class is very small (a small flood in a large scene), for which local or tiled thresholding is preferable.
11.5 Cloud computing and data cubes
The analysis of decades of imagery at planetary scale has become feasible through data cubes (multidimensional arrays indexed by x, y, time, and band) and cloud platforms that bring computation to the archive rather than the data to the analyst. Their enabling ingredients are analysis-ready data, cloud-optimised formats, spatio-temporal catalogues (STAC), and parallel processing. Analysts should be aware of trade-offs: the ease of computation makes it tempting to skip validation, the platforms differ in the processing of the underlying data, and reproducibility demands recording exact collection versions, parameters, and dates.
Exercises
Exercise 11.1 The difference image of two NDVI dates has mean \(0.02\) and standard deviation \(0.05\). Using \(k=2\), what are the thresholds for change? If the classified change area is 3 % of the scene, discuss whether that value is plausible given a Gaussian \(\pm2\sigma\) rule, and what does that imply for the choice of \(k\)?
Exercise 11.2 Between two dates, a pixel’s red and NIR reflectances change by \(\Delta_{\text{red}}=+0.03\) and \(\Delta_{\text{NIR}}=-0.04\). Compute the CVA magnitude and direction (angle measured from the red axis). Interpret the direction (vegetation loss or gain?).
Exercise 11.3 A harmonic model for NDVI has \(a_1=0.20\) and \(b_1=0.15\) (with \(t\) in years). Compute the amplitude and the phase in degrees, and convert the peak to an approximate day of year.
Exercise 11.4 For the annual series \(2,3,5,4,6\): (a) compute the Mann–Kendall statistic \(S\) and the associated \(Z=(S-1)/\sqrt{\mathrm{Var}(S)}\); (b) compute the Theil–Sen slope.
Exercise 11.5 A SAR scene’s histogram is split by a candidate threshold into non-water pixels (60 %, mean \(-10\) dB) and water pixels (40 %, mean \(-22\) dB). Compute the between-class variance. Convert both means to linear power units and comment on why the threshold should be applied in dB.
Exercise 11.6 Design problem. Design a workflow to estimate, with confidence intervals, the area of forest loss in a country between 2018 and 2024 from optical and SAR imagery. Address the choice of method, the treatment of clouds, the sampling design for validation, and the reporting.
12 Uncertainty, Ethics, and Reproducibility
12.1 Error, accuracy, precision, and uncertainty
Accuracy is closeness to the truth; precision is the closeness of repeated estimates to each other; bias is systematic error, and uncertainty is the quantified doubt about a result that remains after all known corrections. A spatial dataset without a statement of uncertainty is incomplete. Uncertainty enters at every step: sensor noise, calibration, atmospheric correction, geometric registration, resampling, classification, generalisation, and analysis.
12.1.1 Positional accuracy
Positional accuracy is measured by comparing coordinates to an independent source of higher accuracy at well-defined check points. With errors \(\Delta x_k,\Delta y_k\) at \(m\) points,
\[ \mathrm{RMSE}_x=\sqrt{\tfrac1m\sum\Delta x_k^{2}},\quad \mathrm{RMSE}_y=\sqrt{\tfrac1m\sum\Delta y_k^{2}},\quad \mathrm{RMSE}_r=\sqrt{\mathrm{RMSE}_x^{2}+\mathrm{RMSE}_y^{2}}. \tag{12.1}\]
The US National Standard for Spatial Data Accuracy (NSSDA) reports horizontal accuracy at the 95 % confidence level as \(1.7308\,\mathrm{RMSE}_r\) when \(\mathrm{RMSE}_x\approx\mathrm{RMSE}_y\), and vertical accuracy as \(1.9600\,\mathrm{RMSE}_z\), provided that the errors are approximately normal, unbiased, and at least 20 check points are used.
12.1.2 Propagation of uncertainty
For a derived quantity \(y=f(x_1,\dots,x_n)\) computed from independent inputs with standard uncertainties \(\sigma_i\), the first-order (Taylor) propagation law is
\[ \sigma_y^{2}\approx\sum_{i=1}^{n}\left(\frac{\partial f}{\partial x_i}\right)^{2}\sigma_i^{2}, \tag{12.2}\]
with covariance terms added if the inputs are correlated. For NDVI, \(\partial\mathrm{NDVI}/\partial\rho_{\mathrm{NIR}}=2\rho_{\mathrm{Red}}/(\rho_{\mathrm{NIR}}+\rho_{\mathrm{Red}})^{2}\) and \(\partial\mathrm{NDVI}/\partial\rho_{\mathrm{Red}}=-2\rho_{\mathrm{NIR}}/(\rho_{\mathrm{NIR}}+\rho_{\mathrm{Red}})^{2}\); for high-vegetation pixels, where red is small and NIR is large, the uncertainty is dominated by the red band’s error. Where the model is non-linear, spatial, or has correlated errors, Monte Carlo simulation (repeatedly perturbing inputs according to their error distributions, recomputing the output, and summarising the ensemble) gives a more realistic uncertainty, at the cost of computing time. For spatial data, the errors must be simulated with their spatial correlation; independent per-pixel noise underestimates the uncertainty of aggregated results.
12.2 Validation principles
The validation of any product should be independent (reference data not used for calibration or training), representative (probability-based sample covering the range of conditions and classes), adequately sized, spatially and temporally consistent with the product, and transparent, reported with the sampling design, the reference-data quality, and confidence intervals. The Committee on Earth Observation Satellites (CEOS) working group on calibration and validation defines hierarchical validation stages for satellite products, from preliminary comparison at a few sites to systematic evaluation across the global range of conditions. For classified maps, the estimators of Equation 8.10 yield area estimates with their uncertainty; for continuous variables, use bias, RMSE, and \(R^{2}\) (or better, the coefficient of agreement and prediction intervals) against independent reference values, and never evaluate on the training data.
12.3 Ethics in geospatial practice
Geospatial data and analysis can benefit and harm.
- Location privacy. Precise locations of individuals, households, or health events are personal data. Geomasking (random displacement, “donut” masking that enforces a minimum displacement, aggregation to areal units) reduces re-identification risk but also degrades analyses; the choice must be documented and its effect on spatial statistics evaluated. Small-area health counts require suppression rules and \(k\)-anonymity thresholds.
- Surveillance and dual use. High-resolution imagery, mobility traces, and machine learning can enable surveillance of individuals and communities. Consider proportionality, consent, and governance.
- Bias and representation. Training data are unevenly distributed across regions, land-cover types, and social groups; models trained in one context may fail, silently and unfairly, in another. Report where a model was validated and where it was not.
- Data sovereignty and participation. Indigenous and local communities have rights concerning data about their lands and lives. Participatory mapping and community engagement should be planned from the outset.
- Communication. Maps persuade. Choices of projection, classification breaks, colour scheme, and area units can mislead; uncertainty should be visible.
- Environmental footprint of large-scale computation is real and should be weighed against the benefit.
12.4 Reproducibility and open science
A result is reproducible if another analyst, given the same data and description, can obtain it. Practical elements include: a complete written workflow and parameters; versioned inputs with persistent identifiers (dataset version, acquisition IDs, download dates); provenance metadata recording each processing step; scripted rather than manual processing where feasible; literate documents (such as the Quarto documents in which this book is written) that couple the narrative to the analysis; and open licences with proper citation of data and software. The FAIR principles apply equally to data and to models.
12.5 Frontiers
Active areas include foundation models pre-trained on petabytes of imagery and adapted to many tasks; hyperspectral, thermal, and SAR constellations that increase spectral and temporal coverage; the fusion of physical models with machine learning (“physics-informed” approaches); explainable and uncertainty-aware GeoAI; causal inference with observational spatial data; digital twins of cities and ecosystems; and low-cost sensors and UAVs producing citizen-generated data. Each raises the same enduring questions: What is measured? At what scale? With what uncertainty? For whom?
Exercises
Exercise 12.1 For \(\rho_{\mathrm{NIR}}=0.42\) and \(\rho_{\mathrm{Red}}=0.06\), each with standard uncertainty 0.02, compute the standard uncertainty of NDVI by first-order propagation. Which band contributes more, and why?
Exercise 12.2 A check-point survey gives \(\mathrm{RMSE}_x=3\) m and \(\mathrm{RMSE}_y=4\) m. Compute \(\mathrm{RMSE}_r\) and the NSSDA horizontal accuracy at the 95 % level.
Exercise 12.3 A forest-loss map covers 1 000 000 ha. Mapped loss occupies 5 % (stratum 2) and stable forest 95 % (stratum 1). A stratified sample yields: of 200 samples in the mapped-stable stratum, 195 are stable and 5 are loss in the reference; of 100 samples in the mapped-loss stratum, 80 are loss and 20 are stable in the reference. Using Equation 8.10, estimate the true proportion and area of loss, and the user’s and producer’s accuracy for loss. What does the result say about relying on pixel counting?
Exercise 12.4 A district health office wants to publish a map of individual-level case locations of a stigmatised infectious disease at 10 m precision. Identify the risks, propose two mitigation strategies, and describe the effect each would have on the utility of the map.
Exercise 12.5 Capstone project. Choose a spatial problem relevant to your field. Write a two-page protocol that specifies: the question; the conceptual model (object or field); data sources and their CRS, resolution, and uncertainty; preprocessing; analytical method with equations; validation design and metrics; uncertainty propagation; ethical considerations; and the reproducibility plan.
Quick-Reference Formulas
| Topic | Relation |
|---|---|
| Wave relation | \(c=\lambda f\); \(E_\gamma=hc/\lambda\) |
| Wien / Stefan–Boltzmann | \(\lambda_{\max}=2898/T\ (\mu\mathrm{m})\); \(M=\sigma T^{4}\) |
| Lambertian reflectance | \(\rho=\pi L/E\) |
| TOA reflectance | \(\rho=\pi L d^{2}/(E_{\mathrm{sun}}\cos\theta_s)\) |
| Beer–Lambert | \(\tau=e^{-\delta/\cos\theta}\) |
| GSD | \(\mathrm{GSD}\approx H_o\,\beta=H_o\,p/f\) |
| Orbital period | \(T=2\pi\sqrt{a^{3}/\mu}\) |
| UTM zone | \(Z=\lfloor(\ell+180)/6\rfloor+1\); \(\ell_0=6Z-183\) |
| Height relation | \(h=H+N\) |
| Mercator scale | \(k=\sec\varphi\) |
| NDVI | \((\rho_{\mathrm{NIR}}-\rho_{\mathrm{Red}})/(\rho_{\mathrm{NIR}}+\rho_{\mathrm{Red}})\) |
| SAR range resolution | \(\rho_{gr}=c/(2B\sin\theta)\); azimuth \(\approx L/2\) |
| InSAR displacement | \(d_{\mathrm{LOS}}=-\lambda\Delta\phi/(4\pi)\) |
| LiDAR range | \(R=ct/2\) |
| Overall accuracy / kappa | \(\mathrm{OA}=\sum n_{ii}/n\); \(\kappa=(p_o-p_e)/(1-p_e)\) |
| Moran’s \(I\) | see Equation 10.3; \(E[I]=-1/(n-1)\) |
| Semivariance | \(\gamma(h)=\frac{1}{2N(h)}\sum[z(x_i)-z(x_i+h)]^{2}\) |
| Error propagation | \(\sigma_y^{2}=\sum(\partial f/\partial x_i)^{2}\sigma_i^{2}\) |
Answers to Selected Exercises
Answers are given for the computational exercises; conceptual and design problems have no single correct answer and are best discussed in class.
Chapter 1. Exercise 1.2: (a) 1.6 km; (b) 2 km × 1.25 km = 2.5 km² = 250 ha; (c) 0.64 cm.
Chapter 2. Exercise 2.1: (a) zone 37, central meridian \(39^{\circ}\)E; (b) zone 10, central meridian \(123^{\circ}\)W. Exercise 2.2: (a) 111.19 km; (b) 55.60 km (\(=111.19\cos60^{\circ}\)). Exercise 2.3: linear scale 2, areal scale 4; apparent area 1000 km². Exercise 2.4: \(H=h-N=1701.7\) m; using \(h\) overstates elevation by 18.5 m. Exercise 2.5: about 264.8 km on the sphere; the ellipsoidal geodesic is slightly shorter, about 264.1 km (near the equator a degree of latitude is shorter than the spherical value).
Chapter 3. Exercise 3.1: (a) \(6.4\times10^{8}\) bytes = 640 MB = 610.4 MiB; (b) 300 km × 240 km = 72 000 km²; (c) factor \((30/10)^{2}=9\). Exercise 3.2: 16 m². Exercise 3.5: \(D=1-\ln(1890/1200)/\ln(2/10)\approx1.28\); lengths are meaningless without stating the step or scale.
Chapter 4. Exercise 4.1: (a) 0.50 µm; (b) 9.66 µm; \(M=459\) W m\(^{-2}\); ratio \((5800/300)^{4}\approx1.4\times10^{5}\). Exercise 4.2: \(3.61\times10^{-19}\) J = 2.25 eV; 5.36 GHz. Exercise 4.3: (a) 0.819; (b) 0.670, which is 18.1 % lower. Exercise 4.4: \(T_b\approx302.8\) K. Exercise 4.6: threshold \(\lambda/(8\cos\theta_i)\) is 3.66 cm for L-band, so the 2 cm field is smooth (dark), and 0.85 cm for C-band, so it is rough (bright).
Chapter 5. Exercise 5.1: \(T\approx98.9\) min; about 14.6 orbits per day; successive tracks are separated by \(24.7^{\circ}\times111.3\approx2750\) km at the equator, so a 185 km swath covers about 6.7 % of it per day, and about 15 days are needed to fill the gaps (consistent with a 16-day repeat cycle). Exercise 5.2: (a) 29.96 m (\(\approx30\) m); (b) \(2\arctan(92.5/705)\approx14.95^{\circ}\); (c) 2.73 cm. Exercise 5.4: \(\rho_{\mathrm{NIR}}=0.38\) for the 40/60 mix and \(0.40\) for the 50/50 mix; a 0.02 difference is comparable to the sensor noise and the within-class variability, so threshold-based forest-area estimates from coarse pixels are unstable.
Chapter 6. Exercise 6.1: \(L=68.5\), \(\theta_s=28^{\circ}\), \(\rho_{\mathrm{TOA}}\approx0.138\). Exercise 6.2: (a) 6 and 10; (b) RMSE \(\approx0.49\) pixel. Exercise 6.3: 27.5. Exercise 6.4: \(d=80\times100/3000\approx2.7\) mm. Exercise 6.5: \(\cos i=\cos(45^{\circ}+25^{\circ})=0.342\); the correction factor is \(0.707/0.342\approx2.07\), so the values are doubled; a small \(\cos i\) in the denominator makes the correction unstable, and it amplifies noise and error in the DEM.
Chapter 7. Exercise 7.1: NDVI \(=0.75\); SAVI \(\approx0.551\); EVI \(\approx0.579\); NDWI \(\approx0.429\) (positive, hence open water is likely). Exercise 7.2: (a) 14.4; (b) \(G_x=8\), \(G_y=24\), magnitude \(\sqrt{640}\approx25.3\); (c) median filter. Exercise 7.3: \(\lambda=(7\pm\sqrt{17})/2=5.56,\,1.44\); the first component explains 79.4 %; \(\mathbf{v}_1\approx(0.788,\,0.615)^{\top}\). Exercise 7.5: sum is 180; sharpened values \((26.7,\,40.0,\,53.3)\); the ratio operation conserves relative band proportions but not absolute radiometry, since the sum of the outputs is forced to equal the panchromatic value (120).
Chapter 8. Exercise 8.1: \(d_1=0.054\), \(d_2=0.112\), so class 1. Exercise 8.2: OA = 0.85; UA = 0.90 and 0.80; PA = 0.818 and 0.889; \(F_1=0.857\) for class 1; \(\kappa=0.70\). Exercise 8.3: 196 samples (5 %), 784 samples (2.5 %). Exercise 8.4: \(0^{\circ}\) (parallel spectra differing only in brightness), and about \(50.1^{\circ}\). Exercise 8.5: \(0.05f+0.25(1-f)=0.12\) gives \(f_{\text{veg}}=0.65\), so band 2 predicts \(0.38\) versus \(0.36\) observed: a residual of 0.02, showing that the two-endmember model explains the pixel closely but not perfectly.
Chapter 9. Exercise 9.1: (a) 1.0 m slant, 2.0 m ground; (b) 5 m (\(L/2\)). Exercise 9.2: 3.5 fringes \(\times\lambda/2=9.8\) cm. Exercise 9.3: (a) \(\approx1000.5\) m; (b) 17.7 m. Exercise 9.4: \(-13.0\) dB; 0.01; \(-8\) dB is 0.158 and \(-19\) dB is 0.0126, a factor of about 12.6 (11 dB) decrease.
Chapter 10. Exercise 10.1: (a) 14.04; (b) 16.84; a larger \(p\) gives more weight to the nearest sample. Exercise 10.2: \(\partial z/\partial x=\partial z/\partial y=0.0667\); slope \(=5.39^{\circ}\); the surface rises toward the north-east, so the aspect is \(225^{\circ}\) (south-west). Exercise 10.3: \(I=(4/6)(2.5/5)=0.333\); \(E[I]=-0.333\); positive autocorrelation (a gradient) though 4 cells are too few for a significance test. Exercise 10.4: \(\hat\gamma=(1+4+0+9+4)/(2\cdot5)=1.8\); a single pair gives an extremely unreliable estimate. Exercise 10.5: expected 100 m; \(R=0.8\) indicates a tendency toward clustering; facilities follow population, roads, and administrative decisions. Exercise 10.6: \(S=0.5(8)+0.3(4)+0.2(6)=6.4\).
Chapter 11. Exercise 11.1: thresholds \(-0.08\) and \(0.12\); a Gaussian rule flags about 4.6 % of pixels as “change” even without real change, so a 3 % change area is consistent with noise alone and \(k\) should be increased or the method reconsidered. Exercise 11.2: magnitude 0.05; angle \(-53^{\circ}\) (reduced NIR, increased red) indicates vegetation loss. Exercise 11.3: amplitude 0.25; phase \(36.9^{\circ}\), peak at \(0.102\) yr \(\approx\) day 37. Exercise 11.4: \(S=8\), \(\mathrm{Var}(S)=16.7\), \(Z\approx1.71\); Sen slope \(=1.0\) per step. Exercise 11.5: \(\sigma_B^{2}=0.6\cdot0.4\cdot12^{2}=34.6\ \mathrm{dB^{2}}\); \(-10\) dB \(=0.10\) and \(-22\) dB \(=0.0063\) in linear power: the linear distribution is highly skewed, whereas dB is approximately symmetric and better suited to threshold selection.
Chapter 12. Exercise 12.1: \(\sigma_{\mathrm{NDVI}}\approx0.074\); the red band contributes \(0.073\) versus \(0.010\) for NIR because \(\lvert\partial\mathrm{NDVI}/\partial\rho_{\mathrm{Red}}\rvert=3.65\) is much larger than \(0.52\). Exercise 12.2: \(\mathrm{RMSE}_r=5\) m; accuracy at 95 % \(=1.7308\times5\approx8.65\) m. Exercise 12.3: \(\hat p_{\text{loss}}=0.95(5/200)+0.05(80/100)=0.06375\), so \(\hat A_{\text{loss}}=63\,750\) ha (versus 50 000 ha mapped); UA (loss) \(=80/100=0.80\); PA (loss) \(=0.04/0.06375\approx0.63\). Pixel counting underestimates loss by about 22 % here, because the omission of loss in the very large stable stratum outweighs the commission errors in the small loss stratum.
Glossary
Albedo. Fraction of incident solar radiation reflected by a surface, integrated over the solar spectrum and all directions.
Analysis-ready data (ARD). Satellite data processed to a minimum set of requirements (calibration, atmospheric and geometric correction, gridding, quality flags) that permit immediate analysis.
Backscatter. Radiation scattered back toward the radar antenna.
BRDF. Bidirectional reflectance distribution function; describes reflectance as a function of illumination and viewing geometry.
Coherence. In interferometry, the magnitude of the complex correlation between two SAR images.
Datum. A reference frame specifying the ellipsoid and its position and orientation relative to the Earth.
Digital number (DN). The integer value recorded by a sensor for a pixel and band.
Endmember. A spectrally pure material used as a component in a mixture model.
Geoid. The equipotential surface of Earth’s gravity field that best fits mean sea level.
Ground sample distance (GSD). Distance on the ground between adjacent pixel centres.
Hyperspectral. Imagery with tens to hundreds of narrow, contiguous spectral bands.
Kriging. Best linear unbiased interpolation based on the modelled spatial covariance.
Lambertian. Ideal surface with radiance independent of view angle.
LiDAR. Light detection and ranging; active laser measurement of range.
Mixed pixel. A pixel whose area contains more than one land-cover class.
Moran’s I. Statistic of global spatial autocorrelation.
Nugget. Semivariance at vanishing lag; microscale variation and measurement error.
Orthorectification. Geometric correction that removes the effects of relief and sensor tilt.
Panchromatic. Single, broad band covering the visible (and often near-infrared) region.
Radiance. Radiant flux per unit projected area per unit solid angle.
Reflectance. The fraction of incident radiation reflected by a surface.
Resampling. Estimation of pixel values at new positions after geometric transformation.
Semivariogram. Function giving half the expected squared difference between values at a given separation.
Spectral signature. Characteristic variation of reflectance or emission with wavelength.
Speckle. Multiplicative granular noise in coherent imaging systems.
Synthetic aperture radar (SAR). Radar that synthesises a long antenna via platform motion to achieve fine azimuth resolution.
Topology. Spatial relationships (connectivity, adjacency, containment) preserved under continuous deformation.
Further Reading and References
Core textbooks
- Lillesand, T. M., Kiefer, R. W., and Chipman, J. W. Remote Sensing and Image Interpretation. Wiley. A comprehensive treatment of the physics of remote sensing, sensor systems, and interpretation of satellite and aerial imagery.
- Green, K., Congalton, R. G., and Tukman, M. Imagery and GIS: Best Practices for Extracting Information from Imagery. Esri Press. The integration of imagery into GIS workflows for real-world spatial analysis.
- Zandbergen, P. A. Python Scripting for ArcGIS Pro. Esri Press. Automation of geoprocessing workflows; recommended for laboratory work that accompanies this text.
- Tomlin, C. D. Geographic Information Systems and Cartographic Modeling. Prentice-Hall. The foundational text of map algebra and raster-based spatial analysis.
- Burrough, P. A., and McDonnell, R. A. Principles of Geographical Information Systems. Oxford University Press. Spatial data structures, interpolation, fuzzy logic, and geostatistics.
Complementary references
- Campbell, J. B., and Wynne, R. H. Introduction to Remote Sensing. Guilford Press.
- Jensen, J. R. Introductory Digital Image Processing: A Remote Sensing Perspective. Pearson.
- Richards, J. A. Remote Sensing Digital Image Analysis. Springer.
- Schowengerdt, R. A. Remote Sensing: Models and Methods for Image Processing. Academic Press.
- Longley, P. A., Goodchild, M. F., Maguire, D. J., and Rhind, D. W. Geographic Information Science and Systems. Wiley.
- Snyder, J. P. Map Projections: A Working Manual. U.S. Geological Survey Professional Paper 1395.
- Isaaks, E. H., and Srivastava, R. M. An Introduction to Applied Geostatistics. Oxford University Press.
- Fotheringham, A. S., Brunsdon, C., and Charlton, M. Geographically Weighted Regression. Wiley.
- Anselin, L. Spatial Econometrics: Methods and Models. Kluwer.
- Oliver, C., and Quegan, S. Understanding Synthetic Aperture Radar Images. SciTech.
- Congalton, R. G., and Green, K. Assessing the Accuracy of Remotely Sensed Data: Principles and Practices. CRC Press.
- Olofsson, P., Foody, G. M., Herold, M., Stehman, S. V., Woodcock, C. E., and Wulder, M. A. (2014). Good practices for estimating area and assessing accuracy of land change. Remote Sensing of Environment, 148, 42–57.
- Pontius, R. G., and Millones, M. (2011). Death to Kappa: birth of quantity disagreement and allocation disagreement for accuracy assessment. International Journal of Remote Sensing, 32(15), 4407–4429.
- Tobler, W. R. (1970). A computer movie simulating urban growth in the Detroit region. Economic Geography, 46, 234–240.