Sentinel2 Annual GeoMAD¶
Date modified: 8 Apr 2021
Product overview¶
Background¶
GeoMAD is the the Digital Earth Africa (DE Africa) Sentinel2 annual geomedian and triple Median Absolute Deviation product. It is a cloudfree composite of Sentinel2 satellite data compiled for each calendar year.
GeoMAD has two main components:
Geomedian
Median Absolute Deviations (MADs)
The geomedian product combines measurements collected in each year to produce one representative, multispectral measurement for every 10 x 10 square metre unit of the African continent. The end result is a comprehensive dataset that can be used to generate truecolour images for visual inspection of anthropogenic or natural landmarks. The full spectral dataset can be used to develop more complex algorithms. The product leverages Sentinel2’s highfrequency flyovers, with the annual time interval allowing for 70 to 140 satellite passes of any given location. This helps scope out perpetually cloudy areas. For each pixel, invalid data is discarded, and remaining observations are mathematically summarised using the geomedian statistic.
Variations between the geomedian and the individual measurements are captured by the three Median Absolute Deviation (MAD) layers. These are higherorder statistical measurements calculating variation relative to the geomedian. The MAD layers can be used on their own or together with geomedian to gain insights about the land surface and understand change over time.
Specifications¶
The GeoMAD dataset is annual (one value per pixel per band per calendar year) and spans continental Africa. GeoMAD metadata can also be viewed on Digital Earth Africa’s Metadata Explorer.
Table 1: GeoMAD product specifications
Specification 


Product name 
Annual Sentinel2 GeoMAD 
Number of bands 
14 
Cell size  X (metres) 
10 
Cell size  Y (metres) 
10 
Coordinate reference system (CRS) 
EPSG: 6933 
Temporal resolution 
Annual 
Temporal range 
20170101 00:00:00 – 20201231 11:59:59 
Parent dataset 
Sentinel2 Level2A 
Update frequency 
Annual 
Figure 1: Sentinel2 Annual GeoMAD geographic extent
Digital Earth Africa Sentinel2 data is available for the regions shaded in blue. Specific temporal and geographic extents can be explored as an interactive map on the Digital Earth Africa Metadata Explorer.
Table 2: GeoMAD measurements
Band ID 
Description 
Value range 
Data type 
No data value 

B02 
Geomedian B02 (Blue) 



B03 
Geomedian B03 (Green) 



B04 
Geomedian B04 (Red) 



B05 
Geomedian B05 (Red edge 1) 



B06 
Geomedian B06 (Red edge 2) 



B07 
Geomedian B07 (Red edge 3) 



B08 
Geomedian B08 (Near infrared (NIR) 1) 



B8A 
Geomedian B8A (NIR 2) 



B11 
Geomedian B11 (Shortwave infrared (SWIR) 1) 



B12 
Geomedian B12 (SWIR 2) 



SMAD 
Spectral Median Absolute Deviation 



EMAD 
Euclidean Median Absolute Deviation 



BCMAD 
BrayCurtis Median Absolute Deviation 



COUNT 
Number of clear observations 



The GeoMAD dataset has fourteen bands of data which can be subdivided as follows:
Geomedian — 10 bands: The geomedian is calculated using ten spectral bands of Sentinel2 data collected during one year. This results in the 10band annual geomedian provided by the GeoMAD product. The ten geomedian surface reflectance bands are B02, B03, B04, B05, B06, B07, B08, B8A, B11, and B12. Surface reflectance values have been scaled between
1
and10000
to allow for more efficient data storage as unsigned 16bit integers (uint16
). Note the raw Sentinel2 product contains more bands, some of which are not used in GeoMAD.The geomedian band IDs correspond to bands in the parent Sentinel2 Level2A data. For example, the GeoMAD band
B02
contains the annual geomedian of the Sentinel2B02
band.Median Absolute Deviations (MADs) — 3 bands: Deviations from the annual geomedian are quantified through median absolute deviation calculations. The GeoMAD product utilises three MADs, each stored in a separate band: Euclidean MAD (EMAD), spectral MAD (SMAD), and BrayCurtis MAD (BCMAD). Each MAD is calculated using the same ten bands as in the geomedian. SMAD and BCMAD are normalised ratios, therefore they are unitless and their values always fall between
0
and1
. EMAD is a function of surface reflectance but is neither a ratio nor normalised, therefore its valid value range depends on the number of bands used in the geomedian calculation — ten in GeoMAD.Count — 1 band: The number of clear satellite measurements of a pixel for that calendar year. This is around 60 to 70 for most pixels, but doubles at areas of overlap between scenes. “Count” is not incorporated in either the geomedian or MADs calculations. It is intended for metadata analysis and data validation.
Processing¶
All clear observations for the given calendar year are collated from the parent dataset. Cloudy pixels are identified and excluded. The geomedian and MADs calculations are then performed by the hdstats package. 2017 – 2020 GeoMAD datasets use hdstats
version 0.2.
Media and example images¶
Image 1: Mangroves in GuineaBissau. 2019 geomedian, truecolour (RGB).
Credit: Contains modified Copernicus Sentinel data 2019, processed by Digital Earth Africa.
Image 2: Croplands in South Africa. 2019 triple MADs, plotted as RGB.
Credit: Contains modified Copernicus Sentinel data 2019, processed by Digital Earth Africa.
References¶
Roberts, D., Mueller, N., & Mcintyre, A. (2017). Highdimensional pixel composites from Earth observation time series. IEEE Transactions on Geoscience and Remote Sensing, 55(11), 6254–6264. https://doi.org/10.1109/TGRS.2017.2723896
Roberts, D., Dunn, B., Mueller, N. (2018). Open Data Cube products using highdimensional statistics of time series. 86478650. https://doi.org/10.1109/IGARSS.2018.8518312.
Acknowledgements¶
The highdimensional statistics algorithms incorporated in this product are the work of Dr Dale Roberts, Australian National University.
Data access¶
Amazon Web Services S3¶
GeoMAD data can be accessed from the associated S3 bucket.
Table 3: AWS data access details
AWS S3 details 


Bucket ARN 

Product name 

The bucket is in the AWS region afsouth1
(Cape Town). Additional region specifications can be applied as follows:
aws s3 ls region=afsouth1 s3://deafricaservices/gm_s2_annual/
The file paths follow the format:
<productname>/<version>/<x>/<y>/<timeperiod>/<x><y>_<timeperiod>_<band>.<extension>
Table 4: AWS file path convention
File path element 
Description 
Example 


Product name 


Product version 


Tile number in the 


Tile number in the 


Year of data collection followed by period of time and time unit in the format



File name. Combines 
Tile numbering uses 96 km grid squares based on the default CRS for this product, EPSG:6933
, with the origin set at the bottomleftmost corner of the valid EPSG:6933
region.
OGC Web Services (OWS)¶
This product is available through DE Africa’s OWS.
Table 5: OWS data access details
OWS details 


Name 

Web Map Services (WMS) URL 

Web Coverage Service (WCS) URL 

Layer name 

Digital Earth Africa OWS details can be found at https://ows.digitalearth.africa/.
For instructions on how to connect to OWS, see this tutorial.
Open Data Cube (ODC)¶
The GeoMAD dataset is not yet but will soon be accessible in the Digital Earth Africa ODC API, which is available through the Digital Earth Africa Sandbox.
ODC product name: gm_s2_annual
Specific bands of data can be called by using either the default names or any of a band’s alternative names, as listed in the table below. ODC datacube.Datacube.load
commands without specified bands will load all bands; see ODC documentation.
Table 6: GeoMAD ODC band names
Band name 
Alternative names 

B02 
band_02, blue 
B03 
band_03, green 
B04 
band_04, red 
B05 
band_05, red_edge_1 
B06 
band_06, red_edge_2 
B07 
band_07, red_edge_3 
B08 
band_08, nir, nir_1 
B8A 
band_8a, nir_narrow, nir_2 
B11 
band_11, swir_1, swir_16 
B12 
band_12, swir_2, swir_22 
SMAD 
smad, sdev, SDEV 
EMAD 
emad, edev, EDEV 
BCMAD 
bcmad, bcdev, BCDEV 
COUNT 
count 
Product names and band names are casesensitive.
For examples on how to use the ODC API, see the DE Africa example notebook repository.
Technical information¶
Geomedian¶
Pixel composites are created by compiling multiple satellite observations of the same area to form one representative image. They have become a staple to Earth observation research as combining measurements compensates for missing data, such as gaps caused by cloud cover. This results in comprehensive datasets that allow thorough analysis of large areas of interest.
There are many ways of forming this image. GeoMAD uses the geomedian summary statistic to combine a year of data into one scientificallyrigorous composite. It produces one multispectral observation for each pixel on continental Africa.
The geomedian statistic is sometimes referred to as the “geometric median”, as in the figure below. To remove ambiguity with other kinds of statistics with similar names, it is always referred to in DE Africa as the “geomedian”.
Figure 2: Illustrated explanation of the geomedian.
Figure 2 uses a threedimensional example to illustrate how the geomedian is calculated for a single pixel containing multiple measurements of red, blue and green data.
The dataset has \(N\) timesteps of satellite data. In the case of Figure 1, \(N=20\)
Each timestep contains \(p\) bands. In this case, \(p =3\) (one dimension each for red, green, blue)
Therefore each timestep can be represented by a 3dimensional vector \(\mathbf{x}\). A single vector at timestep \(t\), \(\mathbf{x}^{(t)}\), looks like this:
The entire dataset can be represented as the collection of the vectors at every timestep \(\mathbf{x}^{(1)}, \mathbf{x}^{(2)}, \dots \mathbf{x}^{(20)}\)
Projected onto a plane, it will look like the plot in Figure 2.2 — twenty points placed depending on their values of red, green and blue. The geomedian of this pixel is then the point where the Euclidean distance (straightline distance) between all data points is minimised. This is depicted in Figure 2.3; for further information on Euclidean distance, see the section below on the Euclidean MAD.
Let’s call the geomedian point in this example \(\mathbf{m}_\text{example}\). Euclidean distance between a point \(\mathbf{x}\) and the data points \(\mathbf{x}^{(t)}\) is given by \(\lVert \mathbf{x}  \mathbf{x}^{(t)} \rVert\). The point at which all twenty distances are minimised is found by taking the \(\mathrm{argmin}\); the argument of the minima. The geomedian can then be expressed by the following equation:
In this example, the resulting \(\mathbf{m}_\text{example}\), like any \(\mathbf{x}\), will be a 3dimensional vector, with a value each for red, green, and blue.
The geomedian point \(\mathbf{m}\) is not selected from one of the twenty data points. It is a separate vector with unique values that do not necessarily correspond to existing band values from the observation dataset.
The formula for the geomedian of a single pixel can be generalised for \(p\) bands and \(N\) timesteps, as given in Roberts et al, 2017:
In GeoMAD, the ten bands of the geomedian result in a 10dimensional space. It is very difficult to illustrate this example in ten dimensions — but we can instead provide the equation, and form of the result. Here, the number of timesteps \(N\) varies per pixel; pixels on the overlap between satellite swathes may have very large \(N\) (up to around 140 per year), while pixels in the centre of the path are only observed once per flyover, and have \(N\) closer to 70.
This calculation is repeated for every pixel in the specified extent.
The significance of the geomedian is that all bands are considered simultaneously. This maintains the spectral relationship between bands, providing the most representative value. Additionally, the geomedian statistic reduces spatial noise and improves colour balance compared to similar statistics such as the median or medoid.
The geomedian is therefore a vital foundational dataset for applications such as band index analysis on vegetation, water and urban area detection. It is also useful as a feature layer in machine learning algorithms.
DE Africa provides an interactive geomedian visualisation available as a downloadable Jupyter Notebook widget. It graphically compares the geomedian to a median for a 3D example. An example screenshot is shown below.
Figure 3: Geomedian visualisation widget.
Triple Median Absolute Deviations (MADs)¶
By definition, the geomedian smooths variations in the satellite data to pick the most central, representative value. Outlying extreme highs and lows are generally filtered out by the robustness of the geomedian calculation. However, it is still useful to know about the variation within the dataset. This results in a set of secondorder statistics: the Median Absolute Deviations (MADs).
MADs are change statistics based on the geomedian. They show how much variation each pixel underwent in the given timeframe.
Let us break down the acronym “MAD”; as in the title, it stands for median absolute deviation:
This implies we have a collection of measurements
We then find the deviation of each measurement from a baseline value, obtaining one deviation for every measurement
These deviations are all absolute values, so each deviation is equal to or greater than 0
We then find the median, or middle value, of these deviations
This gives us a median absolute deviation
In this case, our “collection of measurements” is satellite data from each flyover. Even for a single pixel, we have multiple measurements in the time axis: one from every pass.
The “baseline value” is the annual geomedian, which provides one multispectral result for each pixel.
The “deviations” here are three different distance or dissimilarity values. There are multiple ways of quantifying how change has occurred, so this product computes three different MADs for use in data analysis. We are calculating, in three separate ways, the deviation between the annual geomedian and a single flyover’s measurement. These three values have been chosen to reflect a range of changes that appear in Earth observation data, and hence this section of the dataset is often referred to as “triple MADs”.
The three MADs used in DE Africa are:
Euclidean MAD, EMAD (based on Euclidean distance)
Spectral MAD, SMAD (based on cosine distance)
BrayCurtis MAD, BCMAD (based on BrayCurtis dissimilarity)
Each will be explained in their own sections below. Example calculations with real numbers are in the appendix.
Note there are many other types of statistical distances and dissimilarities that can be used for median absolute deviation analysis (for example: Manhattan distance, Canberra distance, there are many and they can all be used to calculate a MAD). However, in DE Africa, the terms “triple MADs” or “MADs” are always specifically referring to the three MADs included in the GeoMAD dataset: EMAD, SMAD, and BCMAD.
Euclidean MAD (EMAD)¶
The most logical place to start thinking about any of the MADs is the Euclidean MAD (EMAD). This is because EMAD comes from Euclidean distance, and Euclidean distance can be explained with a physical analogy: it is how we measure straightline distances between points. In our threedimensional world, it may look like this:
Figure 4: Euclidean distance in three dimensions.
In the case of satellite data, we are measuring the Euclidean distance between a pixel’s geomedian value and a single multispectral measurement. The number of dimensions is equal to the number of bands in the data. In the illustration below, \(m\) is the geomedian value and \(\mathbf{x}\) the measured value. In real data, there will be multiple measurements over a time period, so \(t\) is the timestep number, otherwise noted in equations as superscript \((t)\).
For example, if we had three bands of data (red, green and blue), and three timesteps of data, then we can calculate the Euclidean distances as follows:
Figure 5: Euclidean distance in three dimensions over three timesteps.
Each timestep gives a separate Euclidean distance result. Then EMAD is the median of all those distances.
In most real life conditions, there will be more than three timesteps and more than three bands. A general expression of Euclidean distance for \(p\) bands is given as:
Then EMAD for \(N\) timesteps is given by Roberts et al, 2018, as the median of the Euclidean distances from all the timesteps.
In GeoMAD, the MADs are calculated from the same ten bands used in the geomedian, therefore \(p=10\). The result of \(\lVert \mathbf{x}^{(t)}  \mathbf{m} \rVert_{\mathbb{R}^p}\) is a positive scalar, so \(\text{EMAD}_\text{GeoMAD}\) is a positive scalar number. As in the geomedian, \(N\) is dependent on the number of satellite flyovers particular to that pixel.
The maximum possible value for EMAD depends on the value ranges for each of the bands in the dataset. In the case of GeoMAD, which uses annual timescales of ten bands of Sentinel2 data, valid EMAD values range from 0  31623
.
EMAD is useful for showing albedo shifts in satellite spectra.
Spectral MAD (SMAD)¶
The spectral MAD (SMAD) is based on the median absolute deviations in the cosine distance between the geomedian and individual measurements.
In two dimensions, cosine distance can be graphically compared to Euclidean distance by the following figure:
Figure 6: Relative relationships between Euclidean and cosine distances.
In a general sense, cosine distance is related to the angle between the two points \(\theta\), while Euclidean distance is related to the straightline distance between the two points \(d\). Like Euclidean distance, points are more similar when the cosine distance between them is small. The value of the cosine distance is smaller when \(\theta\) is small (i.e. close to 0) or when \(\theta\) is close to 180\(^{\circ}\).
Notice we could have a small cosine distance but a large Euclidean distance; for example, if the angle between the vectors is small, but one is much longer than the other. This is an important property of cosine distance (and thus SMAD) — unlike Euclidean distance, cosine distance is not skewed by the magnitude of the measurements.
Cosine distance is defined more formally as:
For more than two dimensions, we can generalise the cosine distance formula for a single pixel. For a multispectral measurement of \(p\) bands at timestep \(t\), \(\mathbf{x}^{(t)}\), and the geomedian at the same point \(\mathbf{m}\), the cosine distance is:
Then for \(N\) timesteps, SMAD is the median of the cosine distances.
As with the other distances and dissimilarities used in the MADs, this results in a positive scalar value, thus SMAD is a positive scalar. Valid values for SMAD fall between 0
– 1
.
In applications of Earth observation data, SMAD is useful for showing areas of land cover change. One reason is that SMAD is less affected by cloud; unlike EMAD, it is invariant to albedo changes, such as that caused by the diffusion of solar radiation. SMAD can also be used to track water bodies, as water has high variation in reflectance.
BrayCurtis MAD (BCMAD)¶
The BrayCurtis MAD (BCMAD) is calculated from the BrayCurtis dissimilarity. The BrayCurtis dissimilarity emphasises differences in each band between the measurement and the geomedian.
For a single band of satellite data, the BrayCurtis dissimilarity looks remarkably like a normalised band index. For example, if we only had red band data, it might look something like this:
It can be generalised to a multispectral dataset with \(p\) bands:
The BrayCurtis dissimilarity will be maximised at a value of 1
when the measurements in each band are completely different. Conversely, the value of the dissimilarity will be small where each band observation is similar to the geomedian of that band.
As with the other MADs, the BCMAD is found by taking the median of all the BrayCurtis dissimilarities from \(N\) timesteps. For GeoMAD, \(p=10\).
BCMAD takes on values from 0  1
.
Appendix¶
Example: calculating Euclidean distance¶
Let’s take a selection of bands from one pixel, from one timestep. For that pixel, we have both the measurements taken by the single satellite flyover, and the geomedian value.
Table A.1: Example multispectral data — one timestep, four bands
Band 
Surface reflectance measurement \(x^{(t)}\) 
Surface reflectance geomedian \(m\) 

Blue 
1028 
969 
Green 
1468 
1406 
Red 
2176 
2032 
Near Infrared (NIR) 1 
3090 
3078 
Then the Euclidean distance for this pixel at this timestep, \(t\), is:
To then find the EMAD of this dataset, the calculation for Euclidean distance would first need to be repeated for all the other timesteps (not provided in the example data).
Example: calculating cosine distance¶
Using the example multispectral data from Table A.1, we can manually calculate the value of cosine distance for timestep \(t\).
Example: calculating BrayCurtis dissimilarity¶
Using the example multispectral data from Table A.1, we can manually calculate the value of the BrayCurtis dissimilarity for timestep \(t\).