Remote Sensing and Earth Observation

Remote sensing is the measurement of the Earth’s surface from a distance, mostly from satellites and aircraft. Earth observation, often used as a synonym, usually means observation from space. For computer vision it is an unusual imaging domain. The sensor looks down, usually close to vertically, from hundreds of kilometres away. Each pixel stands for a known patch of ground with known coordinates. Images come in many spectral bands, several of them outside visible light, or from radar that sees through cloud and works at night. The same place is imaged again every few days for decades. The tasks are familiar: image classification, semantic segmentation, object detection and change detection. The data under them are not.

This article is the overview of the field for vision practitioners. It covers the problem and its history, from Landsat 1 in 1972 to self-supervised models trained on satellite archives. It then describes the data: platforms and orbits, passive optical and multispectral sensing, synthetic aperture radar (SAR) in outline, the four kinds of resolution, and the processing levels that turn raw counts into surface reflectance. The two public missions most used in research, Landsat and Sentinel, are described from their operators’ pages. The vision tasks and benchmarks follow, then a worked example on a synthetic scene: a vegetation index, what happens to it when pixels grow coarser, and a simple cloud mask. The article ends with evaluation, failure modes and deployment.

Problem

Remote sensing tries to turn images of the Earth into maps and measurements: what covers the land, where it changed, where the buildings, roads, ships and aircraft are, how much vegetation is growing, how far a flood has spread. Zhu et al. (2017), in a review of deep learning for remote sensing, list what makes satellite images different from photographs, and their list is a good summary of the problem.

• Many modalities. Optical sensors (multispectral and hyperspectral) and SAR differ completely in imaging geometry and content, and combining them is a research problem of its own.
• Geolocation. Each pixel corresponds to a coordinate on the ground. This makes it possible to fuse images with maps, other sensors and geotagged data.
• Measurements, not pictures. The data are calibrated geodetic and radiometric measurements of controlled quality, so knowledge of the sensor matters more than in purely data-driven vision.
• Time. Continuous acquisition over decades shifts the analysis from single images to time series.
• Volume. Data arrive in very large and growing volumes, often for the whole planet, so algorithms must be fast and must transfer from one region to the next. The authors cite about 25 petabytes acquired by the Sentinel satellites by 2017.
• Physical quantities. Many applications estimate a physical or biochemical quantity rather than detect or classify objects, and process models and expert priors exist for them.

Two further constraints run through everything below. Clouds hide the ground in optical images. And the pixel is often larger than the objects of interest, so one pixel mixes several surfaces. The example below works through both.

History

Landsat 1. The archive that most of the field is built on starts with a launch on 23 July 1972. NASA describes the satellite as the first mission dedicated to monitoring Earth’s land surfaces. It went up as the Earth Resources Technology Satellite (ERTS-1) and was later renamed Landsat 1 (NASA, Landsat 1). It carried two instruments. The Return Beam Vidicon (RBV) was a set of three television-tube cameras and was meant to be the primary sensor. A power surge on 6 August 1972 led controllers to shut it down, and it was permanently switched off 15 days after launch. The secondary instrument, the Multispectral Scanner System (MSS), became the primary one. It swept the ground with an oscillating mirror and recorded four bands, green, red and two near-infrared, as digital numbers with 6-bit radiometric resolution. Its resolution is given as 80 m, commonly resampled to 60 m, with a ground sampling interval of 57 × 79 m and a scene of 170 × 185 km. The satellite covered the whole planet every 18 days, operated until 6 January 1978 and collected more than 300,000 scenes (NASA, Landsat 1).

Vegetation indices. The MSS data soon gave the field one of its lasting tools. Rouse et al. (1974), working on rangelands of the U.S. Great Plains with ground measurements at ten test sites, computed what they called a Band Ratio Parameter: the difference between the radiances in ERTS-1 bands 5 (red) and 7 (near infrared), divided by their sum. The normalization was meant to remove seasonal differences in sun angle and to reduce the effect of the atmosphere. They also defined a Transformed Vegetation Index as the square root of that parameter plus 0.5, to avoid negative values. They showed that it was correlated with above-ground green biomass, and reported a preliminary indication that it was insensitive below a threshold, so that sparsely vegetated areas could not be characterized with it. Tucker’s later account of this work calls the 1974 quantity the vegetation index (VI) (Tucker, 1979). Tucker compared it with other combinations of red and near-infrared radiance measured with field spectrometers over grass plots. He found the normalized difference, its transformed version and the square root of the infrared-to-red ratio the most significant, and he explained why they work. Chlorophyll absorbs strongly in the red, so red radiance falls as green biomass grows and saturates quickly. In the near infrared, leaves absorb little and scatter light within and between them, so radiance rises with biomass. The normalized difference of the two is today usually called the normalized difference vegetation index (NDVI).

More bands, then free data. NASA records that Landsat added shortwave and thermal infrared bands in 1982. The Thematic Mapper on Landsats 4 and 5 had seven bands, Landsat 7’s Enhanced Thematic Mapper Plus eight, and Landsats 8 and 9 have eleven. A change of policy mattered as much as new sensors: in 2008 the Landsat archive was opened free of charge. NASA describes it as a foundational shift: before, only well-funded institutions and agencies could afford the data (NASA, Landsat Data Overview). Europe’s Copernicus programme followed with a free, full and open data policy for its Sentinel satellites, set in EU law (European Commission, legal notice). Sentinel-1, a radar mission, was first launched in 2014, and Sentinel-2, an optical one, in 2015.

Radar. SAR has its own history. Moreira et al. (2013) trace the principle to an invention by Carl Wiley in 1951. Airborne side-looking radar until then had an azimuth resolution that worsened with range, and synthesizing a long antenna from the platform’s motion removed that limit. Early development in the 1950s and 1960s was mostly military. The launch of Seasat in 1978, the first civilian SAR satellite, together with advances in digital technology, opened the way to civilian Earth observation by radar.

From features to deep learning. Early automated analysis classified each pixel from its spectral values. The STATLOG experiment described in the k-nearest neighbours article is one example: pixels of a Landsat image classified from four spectral bands and their neighbours. With very high-resolution aerial images, classifying whole image patches by land use became a vision problem in its own right. Yang and Newsam (2010) built the UC Merced dataset from one-foot-resolution USGS orthoimagery: 21 land-use classes with 100 images of 256 × 256 pixels each. They classified it with bag-of-visual-words models built on local features. Deep learning reached remote sensing a few years after it transformed mainstream vision, as Zhu et al. (2017) describe: convolutional networks for classification and segmentation, then detectors for overhead objects. With them came benchmarks built from the free archives and from commercial very high-resolution imagery: EuroSAT (first released in 2017), Functional Map of the World, DOTA, xView and SpaceNet (2018), and BigEarthNet (2019). Since about 2022, self-supervised models pretrained on unlabelled satellite archives, such as SatMAE (Cong et al., 2022), have extended the trend that foundation models set in general vision.

Visual Data

Platforms and orbits

Most Earth observation satellites fly in low, near-polar, sun-synchronous orbits. Landsat 8 orbits at 705 km with an inclination of 98.2° and crosses the equator at about 10:12 local time on its descending pass, completing about 14.5 orbits a day (NASA, Landsat 8). Sentinel-2 flies at a mean altitude of 786 km and crosses at 10:30. ESA explains the reason for a sun-synchronous orbit. Each pass sees the ground under nearly the same sun angle, apart from small seasonal variations, so shadows and illumination stay consistent across a time series. The morning time is a compromise between enough sunlight and less cloud (ESA, Sentinel-2 Mission). Aircraft and drones fly lower. They give much finer pixels over smaller areas, at times their operators choose, and overlapping images from them can be turned into 3D models with structure from motion and stereo vision.

Optical satellite sensors build an image line by line as the satellite moves. The MSS used a scanning mirror. Landsat 8’s Operational Land Imager replaced the earlier whiskbroom design with a pushbroom sensor (NASA, Landsat 8). Sentinel-2’s Multi-Spectral Instrument is also pushbroom. It collects rows across its 290 km swath and relies on the satellite’s forward motion to add new rows (ESA, Sentinel-2 Mission). Each row is read from long detector arrays, whose physics the image sensor article describes. One consequence affects vision models. On Sentinel-2 the detectors for different bands sit at different places on the focal plane, so each band sees a given patch of ground at a slightly different time. ESA documents the offsets between bands. A moving object such as a ship or an aircraft should therefore appear displaced from band to band (our reasoning from those offsets).

Passive optical and multispectral sensing

Optical sensors are passive. They measure sunlight reflected from the surface, together with the light scattered by the atmosphere along the way. A multispectral sensor records a handful of bands, each a few tens of nanometres wide, chosen for what they reveal. Sentinel-2’s ESA description is a good example of the reasoning:

• four bands at 10 m: blue, green and red, and a broad near-infrared band;
• six bands at 20 m: four narrow bands in the red edge, where vegetation reflectance rises steeply, and two shortwave infrared bands for snow, ice and cloud detection and vegetation moisture;
• three bands at 60 m, mainly for screening clouds and correcting the atmosphere: aerosols at about 443 nm, water vapour at about 945 nm and cirrus at about 1374 nm.

Hyperspectral sensors record many more and narrower bands, so each pixel carries close to a full spectrum. Because a band is just another channel, a multispectral image is a digital image with more than three channels. The colour conventions of photographs do not apply to it, and an image displayed in “natural colour” is a choice of three bands, not a photograph (see colour spaces for the photographic case).

The physics behind the most used band pair is the one Tucker described: green vegetation is dark in the red and bright in the near infrared. Tucker’s spectra also show bare soil, whose reflectance rises much more gently from red to near infrared. Clouds, as ESA’s cloud-detection documentation puts it, are characterised by high reflectance in the blue (ESA, Sentinel-2 Processing). These differences are what simple band ratios and thresholds exploit.

Synthetic aperture radar in outline

SAR is active. It transmits microwave pulses sideways from the moving platform and records the amplitude and phase of the echoes (Moreira et al., 2013). Because it supplies its own illumination at wavelengths of centimetres, it images by day or night and through cloud. Image brightness depends on the geometry, the roughness and the electrical properties of the surface. Longer wavelengths penetrate further into vegetation, snow and dry soil. The flight direction is called azimuth and the line of sight slant range. Focusing the echoes into an image synthesizes an antenna far longer than the physical one, which gives fine azimuth resolution independent of range.

Two properties of SAR images matter for vision models. First, the side-looking geometry projects a three-dimensional scene onto slant range and azimuth. This causes radar shadow for areas hidden from the radar, and foreshortening and layover, which stretch and compress sloped terrain. Second, the many scatterers within each resolution cell add coherently, which produces speckle. Speckle is a strong, multiplicative, grain-like fluctuation from pixel to pixel whose variance grows with intensity. Moreira et al. stress that, although often called noise, speckle is a physical measurement of structure below the resolution cell. It cannot be reduced by transmitting more power. Multi-look processing averages the intensity non-coherently, improving interpretability at the cost of resolution. Noise models and denoisers designed for optical image noise therefore need rethinking before they are applied to SAR (our reasoning). Polarimetry and interferometry, which use the polarization and the phase of the echoes, are SAR’s further tools. ESA’s Sentinel-1 page gives polarimetric signatures as a way to tell volume scatterers such as forest canopy from surface scatterers such as the sea.

The four resolutions and ground sampling distance

Remote sensing describes a sensor by four resolutions. The definitions below follow ESA’s Sentinel-2 page.

• Spatial resolution is the size on the ground that one detector element represents. The related ground sample distance (GSD), also called ground sampling distance, is the spacing between pixel centres on the ground (a common usage, not ESA’s definition). The xView authors, for example, describe their imagery as having a ground sample distance of 0.3 m. The two need not be equal. Landsat 1’s MSS is described as having 80 m resolution, a ground sampling interval of 57 × 79 m and products resampled to 60 m (NASA, Landsat 1). Like the pixel pitch of a camera, GSD sets which objects can be resolved at all; sampling and aliasing gives the theory. A 3 m car spans 10 pixels at xView’s 0.3 m (Lam et al., 2018) and less than one pixel at Sentinel-2’s 10 m.
• Spectral resolution is the instrument’s ability to distinguish features in the electromagnetic spectrum, which the number, placement and width of the bands determine.
• Temporal resolution is the revisit interval: how often the same place is imaged. A single Sentinel-2 satellite revisits every 10 days, and two satellites phased 180° apart revisit every 5 days at the equator (ESA, Sentinel-2 Mission).
• Radiometric resolution is how finely intensity is quantized, usually given in bits. It ranges from the MSS’s 6 bits through Landsat 7’s 8 bits to 12 bits on Landsat 8 and Sentinel-2 and 14 bits on Landsat 9 (NASA; ESA).

The resolutions trade against each other. A sensor with finer pixels collects fewer photons per pixel and usually covers a narrower swath, so it sees a given place less often. This is one reason why public missions have kept moderate resolution (10–30 m) with wide swaths, while sub-metre imaging comes from commercial satellites that point at chosen targets (our reasoning, consistent with the swaths and resolutions given above).

Processing levels and corrections

Satellite images reach users after a chain of processing that operators divide into levels. Sentinel-2’s chain is documented in detail (ESA, Sentinel-2 Processing).

• Level-1B converts raw counts to calibrated top-of-atmosphere (TOA) radiance in the sensor’s own geometry.
• Level-1C applies radiometric and geometric corrections, including orthorectification with a digital elevation model and registration between bands. It resamples each band onto map-projected tiles and converts radiance to TOA reflectance. It also computes cloud, cirrus and snow masks. Geolocation is refined against a Global Reference Image, a set of mostly cloud-free reference images with a global geolocation uncertainty better than 6 m.
• Level-2A corrects for the atmosphere to give surface reflectance (formerly called bottom-of-atmosphere reflectance). Its Scene Classification layer has three cloud classes, including cirrus, and classes for shadows, cloud shadows, vegetation, non-vegetated land, water and snow.

Two details of the storage format catch out vision code. Sentinel-2 measures at 12 bits, but its products store reflectance as 16-bit integers. Reflectance is multiplied by a quantification value of 10,000 by default, and since processing baseline 04.00 (in operation from 25 January 2022) the stored values are also shifted by a band-dependent offset so that negative values are not truncated. The scale and offset are in each product’s metadata. ESA’s Collection-1 reprocessing has since brought the product format with the offset to the archive back to 2015 (ESA, Sentinel-2 Processing). A model trained on products without the offset and applied to products with it would see a shifted input distribution unless the offset is removed (our reasoning). The other detail is that bands come at different resolutions and are resampled to a common grid before most deep networks see them. The choice of resampling is therefore part of the model’s input (our reasoning).

These corrections are the remote sensing form of image registration and radiometric calibration. They matter more than in ordinary vision because analyses compare pixels across dates and sensors. A change in an uncorrected image may be a change in the atmosphere, the sun angle or the geolocation rather than on the ground.

Public Missions and Data

The two public families of satellites are summarized below from their operators’ pages, as of October 2026. Commercial satellites with sub-metre pixels supply the imagery behind benchmarks such as xView and SpaceNet, under their own licences.

Mission	Operator	Launches	Sensor	Bands and spatial resolution	Revisit	Swath
Landsat 8	NASA and USGS	11 Feb 2013	OLI (pushbroom, 12-bit) and TIRS	OLI: 8 bands at 30 m and a panchromatic band at 15 m; TIRS: 2 thermal bands at 100 m, resampled to 30 m	16 days; 8 days with a second Landsat	185 km
Landsat 9	NASA and USGS	27 Sep 2021	OLI-2 (14-bit) and TIRS-2	as Landsat 8	16 days; 8 days with Landsat 8	185 km
Sentinel-2	ESA, for the EU’s Copernicus	2A: 23 Jun 2015; 2B: 7 Mar 2017; 2C: 5 Sep 2024	MSI (pushbroom, 12-bit)	13 bands: 4 at 10 m, 6 at 20 m, 3 at 60 m	10 days per satellite; 5 days with two	290 km
Sentinel-1	ESA, for the EU’s Copernicus	1A: 3 Apr 2014; 1C: 5 Dec 2024; 1D: 4 Nov 2025	C-band SAR, 5.405 GHz (about 5.55 cm)	four imaging modes, resolution down to 5 m, coverage up to 400 km	6-day repeat with two satellites	up to 400 km

Landsat. The program is managed jointly by NASA and the U.S. Geological Survey (USGS). For Landsat 8, NASA was responsible for the space and launch segments and the USGS for the ground segment and the operation of the mission after launch (NASA, Landsat 8). Landsat 9 was launched into Landsat 7’s orbit and replaced it, keeping the 8-day combined coverage with Landsat 8. Its OLI-2 is a copy of Landsat 8’s OLI except that it quantizes to 14 bits instead of 12 (NASA, Landsat 9). NASA describes the whole archive as more than 4 million scenes and over 50 trillion 30 m pixels, a record that runs from 1972 (NASA, Landsat Data Overview). Landsat data and science products are free of charge. Since 2008, as noted above, that has included the full archive.

Sentinel-2. The mission was designed as two identical satellites in the same orbit, phased 180° apart. As of October 2026, ESA lists Sentinel-2B and Sentinel-2C as the operating pair. Sentinel-2A, nearing the end of its life after a decade, flies 36° from Sentinel-2B in an extension campaign that adds revisits over Europe and over cloudy tropical regions of Africa and South America (ESA, Sentinel-2 Mission). For S2A, ESA gives the equivalent wavelengths (the barycentre of each band’s spectral response) of the 10 m bands as 492.7 nm (B2, blue), 559.9 nm (B3, green), 664.6 nm (B4, red) and 832.8 nm (B8, near infrared), with bandwidths of 65, 35, 31 and 118 nm. The values for the other units differ slightly. The ESA pages give the radiometric accuracy requirement as better than 5%, with a goal of 3%.

Sentinel-1. ESA describes the constellation as two sun-synchronous satellites sharing an orbital plane 180° apart, imaging day and night in C-band regardless of weather. Sentinel-1B stopped delivering data in December 2021 after a power-supply anomaly, and ESA and the European Commission announced the end of its mission on 3 August 2022. Sentinel-1C, launched in 2024, then operated with Sentinel-1A, and ESA describes the constellation as moving to Sentinel-1C and Sentinel-1D, launched in 2025 (ESA, Sentinel-1 Mission).

Open data. The European Commission’s legal notice states that EU law provides users with free, full and open access to Copernicus Sentinel data, with reproduction, distribution, adaptation and combination with other data allowed, and requires a source notice (“Copernicus Sentinel data [Year]”) when the data are distributed (European Commission, legal notice). Open access shaped the benchmarks: EuroSAT’s authors chose Sentinel-2 because it is openly and freely accessible (Helber et al., 2019).

Vision Tasks and Methods

Scene classification and land cover

The simplest task assigns one label to an image patch: residential, forest, river, highway. This is image classification on overhead patches. Its benchmarks include UC Merced, EuroSAT and BigEarthNet (below). A patch often contains several land covers, so BigEarthNet labels each patch with between 1 and 12 classes, which makes the task multi-label classification (Sumbul et al., 2019). Land cover mapping proper produces a class for every pixel and is a form of semantic segmentation. Encoder–decoder networks such as U-Net are widely used for it, and Zhu et al. (2017) review the fully convolutional and patch-based approaches used on aerial benchmarks. Two things set remote sensing apart. The input has more channels than RGB, so a network pretrained on ImageNet needs its first layer adapted or the extra bands dropped; EuroSAT’s ImageNet-pretrained baselines, for example, took three channels at a time: the red, green and blue bands, other three-band combinations, or a single band copied to all three inputs (Helber et al., 2019). And labels often come from existing maps, such as the European CORINE Land Cover inventory behind BigEarthNet, rather than from annotators looking at the images.

Object detection

Overhead object detection finds vehicles, ships, aircraft, buildings and storage tanks in very high-resolution images. Xia et al. (2018) give the reasons why it lagged behind detection in natural images. Objects vary hugely in scale. They appear at arbitrary orientations, because the camera looks down rather than sideways, and gravity no longer keeps them upright in the image. And annotated data were scarce. DOTA therefore annotates objects with oriented quadrilaterals instead of axis-aligned boxes. The general detectors in the object detection article, such as Faster R-CNN and YOLO, are the starting point; DOTA’s baselines included a Faster R-CNN adapted to predict oriented boxes (Xia et al., 2018). xView’s authors stress small objects. Their objects range from 3 m (10 pixels) to more than 3 km (10,000 pixels) at 0.3 m GSD, and the dataset keeps the imbalance between frequent and rare classes found in real scenes (Lam et al., 2018). Buildings and roads are often extracted as polygons and graphs rather than boxes, the setting of SpaceNet (Van Etten et al., 2018).

Change detection and time series

Because the same place is imaged repeatedly, many applications compare dates: deforestation, urban growth, flood extent, damage after a disaster. The simplest comparison differences two co-registered images. Learned methods map both images into a joint representation in which changes stand out more clearly. Zhu et al. (2017) review approaches based on deep belief networks, autoencoders and recurrent networks; one of them detects the annual urban dynamics of four cities from 1984 to 2016 in Landsat data. All of them depend on precise co-registration and on radiometric consistency between dates, which is why the corrections above come first. With 5- and 8-day revisits, dense time series have also become inputs in their own right. Zhu et al. call for architectures that use time jointly with spatial and spectral information, and SatMAE’s temporal embedding, described below, is one response.

Super-resolution and pansharpening

Many satellites record a sharper panchromatic band alongside coarser multispectral bands: Landsat 8 and 9 at 15 m against 30 m (NASA). Pansharpening fuses the two to give multispectral images at the panchromatic resolution. Zhu et al. (2017) describe it as a special case of super-resolution and note that training data for learned pansharpening are abundant, since most such satellites carry a panchromatic band. The same reviewers note that pansharpening has been extended from panchromatic bands to thermal and hyperspectral images. A super-resolved image is an estimate, though, not a measurement. It should not be used where the analysis needs what the sensor actually recorded (our reasoning).

Cloud and shadow masking

Before any of these tasks, cloudy pixels must be found. Operators supply masks (ESA, Sentinel-2 Processing). Sentinel-2’s Level-1C masks are computed at 60 m. A threshold on blue reflectance (band B1 by default) detects opaque cloud, a threshold on the 1375 nm band B10 detects cirrus, and a snow index from bands B4 and B11 keeps bright snow from being flagged as cloud. The masks are then cleaned with a morphological opening. The Level-2A Scene Classification uses a longer series of threshold tests, including NDVI and a snow index, with morphological operations, elevation and land cover data, and the parallax between Sentinel-2’s detectors. ESA’s documentation says its Level-1C criteria are imperfect, with false detections from features bright in the blue band or from geometric errors. Thin cloud, haze and shadow over dark surfaces are the hard cases. Cloud masking can therefore also be posed as a segmentation task for learned models. The example below shows why a simple threshold leaves errors at cloud edges and in shadows.

Self-supervised and foundation models

Labels are expensive in remote sensing, while unlabelled imagery is abundant and free. SatMAE (Cong et al., 2022) adapts the masked autoencoder to satellite data. For time series it adds a temporal embedding and masks patches independently across time. For multispectral data it groups bands, with a separate patch embedding and spectral encoding for each group. The authors reported gains over earlier methods both in supervised performance on benchmarks and when the pretrained model was transferred to land cover classification and segmentation. This follows the general move to self-supervised learning and transfer learning. The specific twist is that the spectral and temporal structure of the data becomes part of the pretext task.

Example

The script below builds a synthetic multispectral scene with known ground truth and uses it to show three things. The first is the NDVI on a fine grid. The second is what happens to pixel classification when the same scene is imaged at a coarser GSD. The third is how a simple cloud mask fails at cloud edges and in shadows.

The scene is 1.2 km on a side on a 2 m grid. It contains fields of crop, sparse grass and bare soil, a river 20 m wide, a village, a road 6 m wide and 300 trees with 6 m crowns. The four bands are named after Sentinel-2’s 10 m bands B2, B3, B4 and B8. The reflectance values for each class, the cloud’s reflectance, the noise level and the NDVI threshold of 0.4 are illustrative values chosen for the example, not measured spectra or operational settings. The index is computed as

$$\mathrm{NDVI} = \frac{\rho_{\mathrm{NIR}} - \rho_{\mathrm{red}}}{\rho_{\mathrm{NIR}} + \rho_{\mathrm{red}}},$$

where $\rho$ is reflectance in each band. This is the normalized difference of Rouse et al. (1974) and Tucker (1979), applied here to reflectance rather than to radiance, which those authors used. It ranges from −1 to 1.

Coarser images are simulated by averaging blocks of fine pixels in each band. This idealized sensor integrates the ground evenly over each pixel, with no blur from optics or motion. Mixed pixels are then unmixed under the linear mixing model: a pixel’s spectrum is the area-weighted sum of the class spectra (“endmembers”). The fractions are solved by non-negative least squares with a sum-to-one row. The cloud is a soft disc of opacity $\alpha$ blended with the surface, and its shadow is the same disc displaced and darkening the surface by up to 70%. The script runs with NumPy 2.3.5 and SciPy 1.16.3, and two runs gave identical output.

import numpy as np
from scipy import ndimage
from scipy.optimize import nnls

rng = np.random.default_rng(0)

# Ground truth: a 1.2 km x 1.2 km scene on a 2 m grid (600 x 600 pixels).
WATER, SOIL, CROP, GRASS, BUILT = range(5)
names = ["water", "soil", "crop", "grass", "built"]
H = W = 600
cls = np.full((H, W), SOIL)
for _ in range(40):                                   # fields of random size and crop
    y, x = rng.integers(0, H - 20, 2)
    h, w = rng.integers(20, 160, 2)
    cls[y:y + h, x:x + w] = rng.choice([CROP, GRASS, SOIL], p=[0.5, 0.3, 0.2])
yy, xx = np.mgrid[0:H, 0:W]
river = np.abs(yy - (300 + 80 * np.sin(xx / 70.0))) < 10      # river 20 m wide
cls[river] = WATER
cls[400:520, 40:180] = BUILT                                  # a village block
cls[:, 452:455] = BUILT                                       # a road 6 m wide
trees = rng.integers(0, H, (300, 2))                          # trees of 6 m crown
for ty, tx in trees:
    disk = (yy[ty-3:ty+4, tx-3:tx+4] - ty)**2 + (xx[ty-3:ty+4, tx-3:tx+4] - tx)**2 <= 9
    patch = cls[ty-3:ty+4, tx-3:tx+4]
    patch[disk & (patch == SOIL)] = CROP

# Illustrative reflectances in four bands named after Sentinel-2's 10 m bands
# B2 (blue), B3 (green), B4 (red), B8 (near infrared); not measured spectra.
refl = np.array([[0.06, 0.05, 0.03, 0.02],    # water
                 [0.10, 0.14, 0.18, 0.25],    # bare soil
                 [0.03, 0.07, 0.04, 0.45],    # crop and trees
                 [0.06, 0.09, 0.12, 0.22],    # sparse grass
                 [0.14, 0.15, 0.16, 0.20]])   # roofs and roads
img = refl[cls] + rng.normal(0, 0.005, (H, W, 4))

def ndvi(b):
    return (b[..., 3] - b[..., 2]) / (b[..., 3] + b[..., 2])

def coarsen(a, k):              # ideal box aggregation to a k-times larger GSD
    s = a.shape
    return a.reshape(s[0] // k, k, s[1] // k, k, *s[2:]).mean(axis=(1, 3))

THR = 0.4                       # NDVI threshold for "vegetation", illustrative
true_veg = (cls == CROP).mean()
print("class NDVI:", ", ".join(f"{n} {ndvi(refl[i]):+.2f}" for i, n in enumerate(names)))
print(f"true crop/tree cover: {100 * true_veg:.1f}%")
print(" GSD  pure px  NDVI>thr cover  unmixed cover  |NDVI(mean)-mean(NDVI)| mean/max")
onehot = np.eye(5)[cls]
A = np.vstack([refl.T, 10 * np.ones(5)])       # 4 bands + sum-to-one row (weighted)
for gsd in (2, 10, 30, 60):
    k = gsd // 2
    b = coarsen(img, k)
    frac = coarsen(onehot, k)                  # true class fractions per pixel
    pure = (frac.max(axis=-1) > 0.999).mean()
    hard = (ndvi(b) > THR).mean()
    f = np.array([nnls(A, np.append(p, 10.0))[0] for p in b.reshape(-1, 4)])
    unmixed = f[:, CROP].mean()
    d = np.abs(ndvi(b) - coarsen(ndvi(img), k))
    print(f"{gsd:3d} m  {100 * pure:5.1f}%      {100 * hard:5.1f}%        "
          f"{100 * unmixed:5.1f}%          {d.mean():.3f} / {d.max():.3f}")

# A cloud and its shadow over the 10 m product. alpha is the cloud's opacity:
# 1 in the core, falling smoothly to 0 at thin edges.
b10 = coarsen(img, 5)
c10 = coarsen(onehot, 5)
y, x = np.mgrid[0:120, 0:120]
alpha = np.clip(1.6 - np.hypot(y - 45, x - 70) / 15, 0, 1)
alpha = ndimage.gaussian_filter(alpha, 3)
shadow = np.roll(alpha, (12, -8), axis=(0, 1))   # displaced by the sun's angle
cloud_refl = np.array([0.45, 0.45, 0.46, 0.48])  # bright and spectrally flat
obs = b10 * (1 - 0.7 * shadow[..., None])
obs = (1 - alpha[..., None]) * obs + alpha[..., None] * cloud_refl

affected = (alpha > 0.02) | (shadow > 0.02)      # at least 2% opacity or shading
cloud = obs[..., 0] > 0.20                       # bright in blue: cloud
cloud_d = ndimage.binary_dilation(cloud, iterations=3)
dark = obs[..., 3] < 0.15                        # dark in NIR: shadow (or water)
water = c10[..., WATER] > 0.5
print(f"\npixels affected by cloud or shadow: {affected.sum()} of {affected.size}")
print("mask                 kept  affected left  clear removed (water)  vegetation (cloud-free)")
for label, m in [("none", np.zeros_like(cloud)), ("blue > 0.20", cloud),
                 ("+ 3-px dilation", cloud_d), ("+ NIR < 0.15", cloud_d | dark)]:
    keep = ~m
    est = (ndvi(obs)[keep] > THR).mean()
    ref = (ndvi(b10)[keep] > THR).mean()
    wrong = m & ~affected                        # clear pixels removed
    print(f"{label:16s} {keep.sum():6d}  {(affected & keep).sum():11d}  "
          f"{wrong.sum():9d} ({(wrong & water).sum():4d})  {100 * est:9.1f}% ({100 * ref:4.1f}%)")

class NDVI: water -0.20, soil +0.16, crop +0.84, grass +0.29, built +0.11
true crop/tree cover: 30.4%
 GSD  pure px  NDVI>thr cover  unmixed cover  |NDVI(mean)-mean(NDVI)| mean/max
  2 m  100.0%       30.4%         31.0%          0.000 / 0.000
 10 m   86.5%       32.5%         30.5%          0.006 / 0.592
 30 m   57.1%       34.4%         30.5%          0.019 / 0.542
 60 m   30.5%       37.2%         30.4%          0.025 / 0.519

pixels affected by cloud or shadow: 2973 of 14400
mask                 kept  affected left  clear removed (water)  vegetation (cloud-free)
none              14400         2973          0 (   0)       30.9% (32.5%)
blue > 0.20       13293         1866          0 (   0)       33.5% (34.3%)
+ 3-px dilation   12950         1523          0 (   0)       34.2% (34.6%)
+ NIR < 0.15      12589         1427        265 ( 260)       35.0% (35.4%)

NDVI separates the classes. With these spectra, crop sits far from the other classes, and a threshold of 0.4 recovers the true cover exactly on the 2 m grid. Sparse grass falls below it, by the choice of spectra, a simple stand-in for the insensitivity to sparse vegetation that Rouse et al. reported.

Coarser pixels bias hard classification. The share of pure pixels falls from 86% at 10 m to 31% at 60 m, so most coarse pixels mix classes. Thresholding them overestimates vegetation, by about 2 percentage points at 10 m and 7 at 60 m. The reason lies in the index. A mixed pixel’s NDVI is not the area-weighted mean of its parts’ NDVIs, because each part’s contribution is weighted by its brightness, and vegetation is the brightest class in the near infrared. With these spectra, a crop–soil mixture crosses the 0.4 threshold when about 32% of it is crop (our calculation from the spectra above). So a pixel can be classed as vegetation while two-thirds of it is soil. The per-pixel gap between the NDVI of the mixture and the mean NDVI is small on average, but it reaches 0.5 to 0.6 at boundaries such as river banks, where dark water meets bright crop. Linear unmixing with the true endmembers recovers the cover to within 0.1 percentage points at 10 m and coarser. It does so because the simulated sensor obeys the linear mixing model exactly. At 2 m it reads 0.6 points high: the noise is not yet averaged down, and since fractions cannot be negative, small errors on pure pixels of other classes add a little crop rather than cancelling (our reasoning). Real pixels add the sensor’s blur, endmembers that vary within a class, and non-linear mixing in canopies, and the endmembers must themselves be estimated. In practice, the reported area of a class depends on the GSD and on how mixed pixels are handled, and a map’s resolution belongs in its description (our reading of the example).

A brightness test misses thin cloud and all shadow. The blue threshold mimics, in simplified form, the opaque-cloud test of Sentinel-2’s Level-1C masks. Without a mask, the cloud lowers the measured vegetation cover by 1.6 percentage points relative to the same pixels without cloud, because a bright, flat spectrum pulls NDVI towards zero. The blue test catches only the opaque core, 1,107 of the 2,973 affected pixels. The 3-pixel dilation reaches into the thin edges, and the residual bias falls to 0.4 points. The NIR test, meant for shadows, removes few extra affected pixels but 260 clear water pixels. Water and deep shadow are both dark in the near infrared. Operational schemes such as ESA’s Scene Classification bring in more tests, elevation data and geometry for this reason. The affected pixels still left after all three steps include faint cloud edges and shadow tails, and also shadowed pixels whose near-infrared reflectance stays above the threshold. They are what learned cloud masks and multi-date methods try to catch. Masking also changes which ground is measured. The cloud-free reference rises from 32.5% to 35.4% as masked pixels are removed, because the cloud and its shadow happened to lie over the river and over ground with little crop. A cloud mask therefore introduces a sampling bias whenever cloud correlates with the surface. In reality it does, with terrain, coasts and season (our reasoning).

Datasets and Benchmarks

Public benchmarks reflect the two sources of imagery: free moderate-resolution archives (Landsat, Sentinel) and commercial very high-resolution images. The table summarizes several from their papers.

Dataset	Year	Imagery	Task	Size and labels
UC Merced Land Use	2010	USGS aerial orthoimagery, one-foot pixels	Scene classification	21 classes × 100 images of 256 × 256 pixels
EuroSAT	2017 (journal 2019)	Sentinel-2, 13 bands	Scene classification	27,000 geo-referenced 64 × 64 patches in 10 classes, from cities in 34 European countries
Functional Map of the World	2018	Commercial satellite, multispectral, temporal sequences	Classification of building and land function	1,047,691 images from 207 countries; 63 categories including “false detection”; rich metadata
DOTA	2018	Google Earth, the GF-2 and JL-1 satellites, and CycloMedia aerial images	Oriented object detection	2,806 images of about 800 × 800 to 4000 × 4000 pixels; 188,282 instances in 15 categories, as quadrilaterals
xView	2018	WorldView-3, 0.3 m GSD	Object detection	Over 1 million objects in 60 classes over more than 1,400 km²; axis-aligned boxes
SpaceNet (challenges 1–3)	2018	WorldView-2 and -3, 30–50 cm	Building footprints, road networks	Rio de Janeiro, Las Vegas, Paris, Shanghai, Khartoum; polygons and road graphs
BigEarthNet	2019	Sentinel-2 Level-2A, 12 bands	Multi-label classification	590,326 patches (120 × 120 pixels at 10 m) from 125 tiles in 10 European countries; labels from CORINE Land Cover 2018

Some details from the papers show how these benchmarks were built. EuroSAT’s authors reported an overall accuracy of 98.57% with a ResNet-50 pretrained on ILSVRC-2012, on an 80/20 split of the RGB version of their dataset. After checking all 27,000 patches by hand, they removed patches that were mislabelled or hidden by cloud, snow or ice and dead pixels, and show examples (Helber et al., 2019, arXiv version). BigEarthNet’s tiles were chosen with less than 1% cloud cover. Even so, the authors identified 70,987 patches fully covered by seasonal snow, cloud or cloud shadow and suggest leaving them out. Winter is the least represented season because winter tiles are cloudier (Sumbul et al., 2019). Their baselines removed those patches and split the rest at random into 60% training, 20% validation and 20% test. A shallow network with three convolutional layers, trained from scratch on the RGB bands, scored higher on precision, recall, F1 and F2 than an ImageNet-pretrained Inception-v2 whose weights were frozen and whose added final layers were fine-tuned on the same bands, and a version using all bands scored higher still. The authors conclude that the remote sensing community can benefit from models pretrained on BigEarthNet rather than on computer vision archives. DOTA’s abstract gives the image size as about 4000 × 4000 pixels, while its text gives a range from about 800 × 800. The table follows the text (Xia et al., 2018, arXiv version). The image sources in the table are from DOTA’s dataset page.

The terms of use differ, and the datasets’ pages give them as follows, as of October 2026. EuroSAT is released under the MIT licence, with the Copernicus terms applying to the underlying Sentinel data. BigEarthNet is released under the Community Data License Agreement – Permissive, version 1.0; its page now describes a version 2.0 of 549,488 pairs of Sentinel-1 and Sentinel-2 patches from 115 Sentinel-2 tiles, with pixel-level reference maps in addition to the multi-labels. Functional Map of the World uses its own Functional Map of the World Challenge Public License. DOTA’s images and annotations may be used for academic purposes only, with commercial use prohibited, and Google Earth’s terms apply to its Google Earth images. xView is distributed to registered users under the terms and conditions of its challenge, and SpaceNet’s entry in the AWS open data registry gives its licence as various.

Challenges

• Mixed pixels and scale. At public-mission resolutions many objects are smaller than a pixel, and class boundaries cut through pixels. As the example shows, hard labels at coarse GSD bias area estimates. At sub-metre resolution the problem reverses: objects span hundreds of pixels at arbitrary orientations and very different scales (Xia et al., 2018; Lam et al., 2018).
• Clouds and illumination. Optical sensors see the ground only through clear sky. Cloud masks have errors at thin edges and in shadows, and cloudy seasons and regions are underrepresented in training data. BigEarthNet’s shortage of winter images is an example (Sumbul et al., 2019).
• Labels. Annotation needs specialists and is expensive (Cong et al., 2022). Labels taken from maps such as CORINE carry the map’s own level of detail, date and errors (our reasoning). DOTA’s objects were annotated by experts in aerial image interpretation (Xia et al., 2018).
• Distribution shift. Models trained in one region, season or sensor can degrade in another: the spectra of a class plausibly vary with climate, soil and crop calendar, and sensors differ in bands, resolution and processing (our reasoning). This is dataset bias in a geographic form. Xia et al. (2018), citing earlier work, note that sensor resolution and variety cause dataset biases, and collected DOTA from several sensors for that reason.
• Data volume and channels. Global coverage at 10 m is a large amount of data per date. Inputs have 4 to 13 or more bands at mixed resolutions, so architectures and pretraining built for three-channel photographs need adapting (Zhu et al., 2017).
• Physics and geometry. Atmospheric effects, viewing and sun angles, terrain relief and, for SAR, speckle and layover all change the measured values without any change on the ground.

Evaluation

Each task brings its usual metrics, with conventions specific to the field.

• Classification is scored by overall accuracy and the confusion matrix; EuroSAT reports overall accuracy (Helber et al., 2019). BigEarthNet’s baselines are scored with precision and recall, F1 and F2, which its authors describe as widely used for multi-label classification (Sumbul et al., 2019).
• Detection follows the general metrics. DOTA adopts the PASCAL VOC calculation of average precision and reports baselines against both horizontal boxes and oriented ground truth (Xia et al., 2018).
• Building footprints in SpaceNet count a proposal as a true positive when its intersection over union with a footprint is at least 0.5, then compute precision, recall and F1 per city and average the F1 over cities (Van Etten et al., 2018).
• Road networks show why pixel metrics can mislead. SpaceNet’s authors give an example in which a road mask with wrong widths but no breaks scores lower on pixel IoU and F1 than a mask with correct widths and a gap, although only the first could be used for routing. They therefore introduced the Average Path Length Similarity (APLS) metric, which compares shortest-path lengths between the true and proposed road graphs on a scale from 0 to 1 (Van Etten et al., 2018).

How the test set is chosen matters as much as the metric. Neighbouring pixels and patches are strongly correlated, and patches cut from the same scene share acquisition conditions. A random split of patches therefore tends to give more optimistic scores than a split by region or by date, which is what a deployed map faces (our reasoning). The baseline splits of both EuroSAT (80/20, applied class by class) and BigEarthNet (60/20/20) were random splits of patches (Helber et al., 2019; Sumbul et al., 2019). For area estimates, a classifier’s pixel accuracy is not the same as the accuracy of the area it reports, as the example shows.

Failure Modes and Risks

• Undetected cloud, haze and shadow pass into downstream products as false change, false loss of vegetation or missing objects. Overcautious masks instead remove valid ground, water in the example, and bias what is measured.
• Inconsistent preprocessing across dates or product versions creates spurious change. Examples are TOA against surface reflectance, different processing baselines, or Sentinel-2’s offset introduced in 2022.
• Geographic shift. A model can perform well on its benchmark and fail in regions, seasons or landscapes outside its training data, often silently, because no labels exist there to measure the drop.
• Resolution misread. Reporting areas or counts derived from coarse pixels as if they were exact. Super-resolved images are a related risk when treated as measurements.
• Sensitive uses. Very high-resolution imagery can reveal people’s homes, vehicles and activities. Detection datasets such as xView were developed with defence organizations, as their author affiliations show (Lam et al., 2018). Questions of privacy, dual use and the licences of commercial imagery accompany such work, and the rules differ by country.

Deployment Considerations

• Use the operator’s levels and masks. Start from surface reflectance products, such as Sentinel-2 Level-2A, where the analysis compares dates. Read the product metadata for scale factors and offsets, and treat the operators’ cloud masks as a baseline to evaluate, not as ground truth.
• Keep the geometry. Work in the product’s map projection and keep the georeferencing through the pipeline, so outputs remain maps that can be overlaid, validated and compared over time.
• Match revisit to the task. A 5-day nominal revisit gives fewer usable images in cloudy regions and seasons. Radar from Sentinel-1 can fill gaps where optical images are cloud-covered, at the cost of a different modality.
• Validate where it will be used. Hold out regions and dates, not random patches, and collect reference data in the target area before trusting a model there.
• Licences and attribution. Landsat and Sentinel data are free and open. Copernicus requires a source notice when its data are distributed. Commercial imagery and some datasets have restrictive terms.