GeoMAD cloudfree composites¶
Date modified: March 2022
Service overview¶
Background¶
GeoMAD is the Digital Earth Africa (DE Africa) surface reflectance geomedian and triple Median Absolute Deviation data service. It is a cloudfree composite of satellite data compiled over specific timeframes. This service is ideal for:
Longerterm time series analysis
Cloudfree imagery
Statistical accuracy
GeoMAD has two main components:
Geomedian
Median Absolute Deviations (MADs)
The geomedian component combines measurements collected over the specified timeframe to produce one representative, multispectral measurement for every pixel 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.
For each pixel, invalid data is discarded, and the remaining observations are mathematically summarised using the geomedian statistic. Flyover coverage provided by collecting data over a period of time also helps scope intermittently cloudy areas.
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.
Calculating GeoMAD over different timeframes and sensors provides a range of insights to the environment. An annual timeframe allows better correction for cloud cover and reduces artifacts for comparison over multiple years. A shorter timeframe, for example, sixmonth or threemonth blocks, better captures seasonal variation within one year, and can be used to compare equivalent periods from different years. Landsat sensors allow full utility of the surface reflectance archive dating back to 1984, while more recent Sentinel2 data provides higherfrequency flyovers and better resolution.
The Digital Earth Africa GeoMAD service currently provides annual, sixmonth semiannual and rolling monthly (3month moving window) datasets, with separate services for Landsat and Sentinel2 sensors.
Jupyter Notebooks that demonstrate loading and using GeoMAD and Rolling GeoMAD datasets in the Sandbox are also available.
Specifications¶
The GeoMAD service spans all of Africa and is currently available in three timeframes:
Annual: observations from one calendar year are summarised into one set of measurements
Sentinel2 Annual GeoMAD
Landsat 8 Annual GeoMAD
Landsat 5 & 7 Annual GeoMAD
Landsat 8 & 9 Annual GeoMAD
Semiannual: observations are summarised for each half of the calendar year, giving one set of measurements for January–June, and one for July–December
Sentinel2 Semiannual GeoMAD
Rolling GeoMAD: observations are summarised over rolling threemonth periods starting on the first day of each calendar month, giving a set of new measurements each month e.g JanMar, FebApr, MarMay and so on.
Sentinel2 Rolling Monthly GeoMAD
GeoMAD metadata can also be viewed on Digital Earth Africa’s Metadata Explorer.
Table 1: Sentinel2 GeoMAD service specifications
Specification 


Service name 
Sentinel2 Annual GeoMAD 
Sentinel2 Semiannual GeoMAD 
Sentinel2 Rolling Monthly GeoMAD 
Service status 
Operational 
Operational 
Operational 
Number of bands 
14 
14 
14 
Cell size  X (metres) 
10 
10 
10 
Cell size  Y (metres) 
10 
10 
10 
Coordinate reference system (CRS) 
EPSG: 6933 
EPSG: 6933 
EPSG: 6933 
Temporal resolution 
Annual 
6monthly (JanJun, JulDec) 
Monthly (3month rolling window) 
Temporal range 
2017 – current 
2017 – current 
2019 – current 
Parent dataset 
Sentinel2 Level2A 
Sentinel2 Level2A 
Sentinel2 Level2A 
Update frequency 
Annual 
6 months 
Monthly 
Update latency 
2 months from end of previous year 
2 months from end of previous half year 
1 month from end of previous period 
Table 2: Landsat GeoMAD service specifications
Specification 


Service name 
Landsat 8 & 9 Annual GeoMAD 
Landsat 8 Annual GeoMAD 
Landsat 5 & 7 Annual GeoMAD 
Service status 
Operational 
Operational 
Operational 
Number of bands 
10 
10 
10 
Cell size  X (metres) 
30 
30 
30 
Cell size  Y (metres) 
30 
30 
30 
Coordinate reference system (CRS) 
EPSG: 6933 
EPSG: 6933 
EPSG: 6933 
Temporal resolution 
Annual 
Annual 
Annual 
Temporal range 
2021 – 
2013 – 2020 
1984 – 2012 
Parent dataset 
Landsat Level2 Collection 2 
Landsat Level2 Collection 2 
Landsat Level2 Collection 2 
Update frequency 
Annual 
Annual 
NA 
Update latency 
2 months from end of previous year 
2 months from end of previous year 
NA 
The Landsat 5 & 7 Annual GeoMAD service combines data from both the Landsat 5 and Landsat 7 sensors. It is produced for the 1984 – 2012 archive only. Similarly, the Landsat 8 & 9 Annual GeoMAD came into production after Landsat 9 data became available in late 2021.
Figure 1: GeoMAD geographic extent
Digital Earth Africa GeoMAD 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.
Sentinel2, Landsat 8, and Landsat 8 & 9 GeoMAD services are generated from parent datasets with frequent satellite flyovers consistent with the geographic extent shown in Figure 1. While the Landsat 5 & Landsat 7 Annual GeoMAD also covers the same area, there are some limitations associated with the sparsity of the Landsat archive over Africa, particularly before the launch of Landsat 7 in 1999.
This is visually displayed in the count
imagery in Figure 2. Northern and southern Africa have the best coverage over the archive period. Central and western Africa show lower data counts, and therefore less years where a geomedian and MADs have been calculated. It is recommended to consult the Digital Earth Africa Metadata Explorer for precise data availability of the Landsat 5 and 7 GeoMAD.
Figure 2: Landsat 5 & 7 GeoMAD data extent
Left: An animation showing data count over Africa for the years 1984 to 2012. Darker blue indicates higher counts, while lighter blue is lower counts. Areas with very low or zero counts in one year do not have enough data to calculate a geomedian or MADs for that year.
Right: A static image showing the number of years for which a geomedian has been calculated. In the 29 years between 1984 and 2012, northern and southern Africa show almost complete geomedian coverage, while central and western Africa have a lower geomedian count associated with less data acquisition over the time period.
Table 3.1: Sentinel2 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 



Table 3.2: Landsat 8 and Landsat 8 & 9 GeoMAD measurements
Band ID 
Description 
Value range 
Data type 
No data value 

SR_B2 
Geomedian SR_B2 (Blue) 



SR_B3 
Geomedian SR_B3 (Green) 



SR_B4 
Geomedian SR_B4 (Red) 



SR_B5 
Geomedian SR_B5 (NIR) 



SR_B6 
Geomedian SR_B6 (SWIR 1) 



SR_B7 
Geomedian SR_B7 (SWIR 2) 



SMAD 
Spectral Median Absolute Deviation 



EMAD 
Euclidean Median Absolute Deviation 



BCMAD 
BrayCurtis Median Absolute Deviation 



COUNT 
Number of clear observations 



Table 3.3: Landsat 5 & Landsat 7 GeoMAD measurements
Band ID 
Description 
Value range 
Data type 
No data value 

SR_B1 
Geomedian SR_B1 (Blue) 



SR_B2 
Geomedian SR_B2 (Green) 



SR_B3 
Geomedian SR_B3 (Red) 



SR_B4 
Geomedian SR_B4 (NIR) 



SR_B5 
Geomedian SR_B5 (SWIR 1) 



SR_B7 
Geomedian SR_B7 (SWIR 2) 



SMAD 
Spectral Median Absolute Deviation 



EMAD 
Euclidean Median Absolute Deviation 



BCMAD 
BrayCurtis Median Absolute Deviation 



COUNT 
Number of clear observations 



Bands can be subdivided as follows:
Geomedian  10 bands (Sentinel2), 6 bands (Landsat 5/7/8/9): The geomedian is calculated using the spectral bands of data collected during the specified time period. Surface reflectance values have been scaled between
1
and10000
to allow for more efficient data storage as unsigned 16bit integers (uint16
). Note parent datasets often contain 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 Annual GeoMAD band
B02
contains the annual geomedian of the Sentinel2B02
band.
Median Absolute Deviations (MADs)  3 bands: Deviations from the geomedian are quantified through median absolute deviation calculations. The GeoMAD service 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 for Sentinel2 and 20 for Landsat 8 annually, 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 time period are collated from the parent dataset. Cloudy pixels are identified and excluded. The geomedian and MADs calculations are then performed by the hdstats package.
Media and example images¶
Image 1: Animations over northern Africa. 1984 – 2012 Landsat 5 & 7 Annual geomedian, truecolour (RGB).
Left: Agricultural area expansion near Sadat City, Egypt.
Right: Urbanisation in Tangier, Morocco.
Credit: Contains Landsat Level2 Surface Reflectance Science Product courtesy of the U.S. Geological Survey, processed by Digital Earth Africa.
Image 2: Mangroves in GuineaBissau. 2019 Sentinel2 Annual geomedian, truecolour (RGB).
Credit: Contains modified Copernicus Sentinel data 2019, processed by Digital Earth Africa.
Image 3: Croplands in South Africa. 2019 Sentinel2 Annual triple MADs, plotted as RGB.
Credit: Contains modified Copernicus Sentinel data 2019, processed by Digital Earth Africa.
Image 4: Irrigated fields in Egypt. JanJun 2020 (right), JulDec 2020 (left). 2020 Sentinel2 Semiannual geomedian, truecolour (RGB).
Credit: Contains modified Copernicus Sentinel data 2020, processed by Digital Earth Africa.
Image 5: The Okavango Delta, Botswana. 2020 Landsat 8 Annual triple MADs, plotted as RGB.
Credit: Contains Landsat Level2 Surface Reflectance Science Product courtesy of the U.S. Geological Survey, processed by Digital Earth Africa.
Image 6: Animations 2020–2022 using Sentinel2 Rolling GeoMAD, truecolour (RGB)
The dates shown are the midpoints of the Rolling GeoMAD.
Credit: Contains modified Copernicus Sentinel data 2022, processed by Digital Earth Africa.
Farmland in northern Egypt
Wetland in Guinea Bissua
Known limitations¶
The Landsat 8 (& 9) GeoMAD has a known issue with data quality over marine regions. The GeoMAD algorithm uses pixel quality information from the input data to identify and mask pixels with poor quality obervations. Landsat 8 & 9 analysis ready satellite images over the ocean often contain negative surface reflectance values, and the GeoMAD masking procedures remove pixels where any negative values occur. Thus, in regions where pixels are persistently negative throughout the year, the GeoMAD product will contain a nodata value. An example of this can be seen in Image 7 below where a shallow marine system contains nodata values in the GeoMAD because the NIR band values in the input data are persistently negative.
Image 7: Shallow marine system, Mauritius. 2021 Landsat 8 Annual GeoMAD (left), versus Sentinel2 GeoMAD (right), plotted as RGB.
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), 62546264. 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 service 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 service 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)¶
GeoMAD datasets are accessible in the Digital Earth Africa ODC API, which is available through the Digital Earth Africa Sandbox.
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.1: Sentinel2 GeoMAD ODC product and band names
ODC product name: gm_s2_annual
, gm_s2_semiannual
, gm_s2_rolling
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 
Table 6.2: Landsat 8 GeoMAD ODC product and band names
ODC product name: gm_ls8_annual
Band name 
Alternative names 

SR_B2 
band_2, blue 
SR_B3 
band_3, green 
SR_B4 
band_4, red 
SR_B5 
band_5, nir 
SR_B6 
band_6, swir_1 
SR_B7 
band_7, swir_2 
SMAD 
smad, sdev, SDEV 
EMAD 
emad, edev, EDEV 
BCMAD 
bcmad, bcdev, BCDEV 
COUNT 
count 
Table 6.3: Landsat 5 & 7 GeoMAD ODC product and band names
ODC product name: gm_ls5_ls7_annual
Band name 
Alternative names 

SR_B1 
band_1, blue 
SR_B2 
band_2, green 
SR_B3 
band_3, red 
SR_B4 
band_4, nir 
SR_B5 
band_5, swir_1 
SR_B7 
band_7, swir_2 
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 six months or 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 2, \(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:
:nbsphinxmath:`begin{align*} mathbf{x}^{(t)} = begin{bmatrix}
x^{(t)}_{red} \ x^{(t)}_{green}\ x^{(t)}_{blue}
end{bmatrix} end{align*}`
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:
\begin{align*} \mathbf{m}_\text{example} = \underset{\mathbf{x} \in \mathbb{R}^3}{\mathrm{argmin}} \sum^{20}_{t=1} \lVert \mathbf{x}  \mathbf{x}^{(t)} \rVert \end{align*}
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.
:nbsphinxmath:`begin{align*} mathbf{m}_text{example} =begin{bmatrix}
m_{red} \ m_{green}\ m_{blue}
end{bmatrix} end{align*}`
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:
\begin{align*} \mathbf{m} = \underset{\mathbf{x} \in \mathbb{R}^p}{\mathrm{argmin}} \sum^{N}_{t=1} \lVert \mathbf{x}  \mathbf{x}^{(t)} \rVert \end{align*}
In the Sentinel2 GeoMADs, the ten bands of the geomedian result in a 10dimensional space  in the Landsat GeoMADs, the six geomedian bands give a 6dimensional space. It is difficult to illustrate this example in six (or ten!) dimensions, but we can instead provide the equation, and the 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 for Sentinel2 annual services), while pixels in the centre of the path are only observed once per flyover, and have \(N\) closer to 70.
:nbsphinxmath:`begin{align*} mathbf{m}_text{GeoMAD} = underset{mathbf{x} in mathbb{R}^{10}}{mathrm{argmin}} sum^{N}_{t=1} lVert mathbf{x}  mathbf{x}^{(t)} rVert = begin{bmatrix}
m_{B02} \ m_{B03}\ m_{B04}\ m_{B05} \ m_{B06}\ m_{B07}\ m_{B08} \ m_{B8A}\ m_{B11}\ m_{B12}
end{bmatrix} end{align*}`
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 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.
To interact with the widget, log in to the Digital Earth Africa Sandbox and navigate to the Datasets folder > GeoMAD.
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 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 service computes three different MADs for use in data analysis. We are calculating, in three separate ways, the deviation between the 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 services, 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:
\begin{align*} &\text{Multispectral Euclidean distance for timestep }t \\ =& \sqrt{ \left( x^{(t)}_{\text{band 1}}  m_{\text{band 1}} \right)^2 + \left( x^{(t)}_{\text{band 2}}  m_{\text{band 2}} \right)^2 + \dots + \left( x^{(t)}_{\text{band p}}  m_{\text{band p}} \right)^2 }\\ =& \lVert \mathbf{x}^{(t)}  m \rVert_{\mathbb{R}^p} \end{align*}
Then EMAD for \(N\) timesteps is given by Roberts et al, 2018, as the median of the Euclidean distances from all the timesteps.
\begin{align*} \text{EMAD} = \text{median} \left( \left\{ \lVert \mathbf{x}^{(t)}  \mathbf{m} \rVert_{\mathbb{R}^p}, t = 1, \dots , N \right\} \right) \end{align*}
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.
\begin{align*} \text{EMAD}_\text{GeoMAD} = \text{median} \left( \left\{ \lVert \mathbf{x}^{(t)}  \mathbf{m} \rVert_{\mathbb{R}^{10}}, t = 1, \dots , N \right\} \right) \end{align*}
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 at maximum 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:
\begin{align*} \text{Cosine distance (two dimensions)} = 1  \frac{x_1 y_1 + x_2 y_2}{ \left( \sqrt{ \left( x_1\right) ^2 + \left( x_2\right) ^2 } \right) \left( \sqrt{ \left( y_1\right) ^2 + \left( y_2\right) ^2 } \right)} \end{align*}
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:
\begin{align*}\small &\text{Multispectral cosdist}\left( \mathbf{x}^{(t)}, m \right) \text{ for timestep } t \\ &= 1  \frac{ \mathbf{x}^{(t)} \cdot \mathbf{m} }{ \lVert \mathbf{x}^{(t)} \rVert \ \lVert \mathbf{m} \rVert} \ \text{ for } \mathbf{x}^{(t)}, \mathbf{m} \in \mathbb{R}_{p}\\ &= 1  \left( \frac{\left( x_{\text{band 1}}^{(t)} \right) \left(m_{\text{band 1}} \right) + \left( x_{\text{band 2}}^{(t)} \right) \left(m_{\text{band 2}} \right) + \cdots + \left( x_{\text{band p}}^{(t)} \right) \left(m_{\text{band p}} \right)}{ \left(\sqrt{\left( x_{\text{band 1}}^{(t)} \right)^2 + \cdots+ \left( x_{\text{band p}}^{(t)} \right)^2} \right) \left( \sqrt{\left( m_{\text{band 1}} \right)^2 + \cdots+ \left( m_{\text{band p}} \right)^2 } \right)} \right) \end{align*}
Then for \(N\) timesteps, SMAD is the median of the cosine distances.
\begin{align*} \text{SMAD} = \text{median} \left( \left\{ \text{cosdist}\left( \mathbf{x}^{(t)}, \mathbf{m} \right), t = 1, \dots , N \right\} \right) \end{align*}
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:
\begin{align*} \text{Singleband BrayCurtis dissimilarity at timestep }t = \frac{\left x_{\text{red}}^{(t)}  m_{\text{red}}\right}{ \left x_{\text{red}}^{(t)} + m_{\text{red}} \right } \end{align*}
It can be generalised to a multispectral dataset with \(p\) bands:
\begin{align*} &\text{Multispectral BrayCurtis dissimilarity for timestep }t\\ &= \frac{\left x_{\text{band 1}}^{(t)}  m_{\text{band 1}}\right + \left x_{\text{band 2}}^{(t)}  m_{\text{band 2}} \right + \dots + \left x_{\text{band p}}^{(t)}  m_{\text{band p}} \right }{ \left x_{\text{band 1}}^{(t)} + m_{\text{band 1}} \right + \left x_{\text{band 2}}^{(t)} + m_{\text{band 2}} \right + \dots + \left x_{\text{band p}}^{(t)} + m_{\text{band p}} \right} \end{align*}
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\).
\begin{align*} \text{BCMAD} = \text{median} \left( \left\{ \frac{\left \mathbf{x}^{(t)}  \mathbf{m} \right_{\mathbb{R}^p}}{\left \mathbf{x}^{(t)} + \mathbf{m} \right _{\mathbb{R}^p}}, t = 1, \dots , N \right\} \right) \end{align*}
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:
\begin{align*} &\text{Euclidean distance} \\ \tiny &= \sqrt{\left( x^{(t)}_{\text{band 1}}  m_{\text{band 1}} \right)^2 + \left( x^{(t)}_{\text{band 2}}  m_{\text{band 2}} \right)^2 + \dots + \left( x^{(t)}_{\text{band p}}  m_{\text{band p}} \right)^2 }\\ &= \sqrt{\left( x^{(t)}_{\text{red}}  m_{\text{red}} \right)^2 + \left( x^{(t)}_{\text{green}}  m_{\text{green}} \right)^2 + \left( x^{(t)}_{\text{blue}}  m_{\text{blue}} \right)^2 + \left( x^{(t)}_{\text{nir1}}  m_{\text{nir1}} \right)^2}\\ &= \sqrt{\left( 2176  2032 \right)^2 + \left( 1468  1406 \right)^2 + \left(1028  969 \right)^2 + \left( 3090  3078 \right)^2}\\ &= 167.9 \end{align*}
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\).
\begin{align*} &\text{Cosine distance} \\ &= 1  \left( \frac{\left( x_{\text{band 1}}^{(t)} \right) \left(m_{\text{band 1}} \right) + \left( x_{\text{band 2}}^{(t)} \right) \left(m_{\text{band 2}} \right) + \cdots + \left( x_{\text{band p}}^{(t)} \right) \left(m_{\text{band p}} \right)}{ \left(\sqrt{\left( x_{\text{band 1}}^{(t)} \right)^2 + \cdots+ \left( x_{\text{band p}}^{(t)} \right)^2} \right) \left( \sqrt{\left( m_{\text{band 1}} \right)^2 + \cdots+ \left( m_{\text{band p}} \right)^2 } \right)} \right) \\ &= 1 \\ & \ \frac{\left( x_{\text{blue}}^{(t)} \right) \left(m_{\text{blue}} \right) + \left( x_{\text{green}}^{(t)} \right) \left(m_{\text{green}} \right) + \left( x_{\text{red}}^{(t)} \right) \left(m_{\text{red}} \right) + \left( x_{\text{nir1}}^{(t)} \right) \left(m_{\text{nir1}} \right)}{\sqrt{\left( x_{\text{blue}}^{(t)} \right)^2 + \left( x_{\text{green}}^{(t)} \right)^2 + \left( x_{\text{red}}^{(t)} \right)^2+ \left( x_{\text{nir1}}^{(t)} \right)^2}\sqrt{\left( m_{\text{blue}} \right)^2 + \left( m_{\text{green}} \right)^2 + \left( m_{\text{red}} \right)^2 + \left(m_{\text{nir1}}\right)^2 } }\\ &= 1 \\ & \ \frac{\left( 1028 \right) \left(969 \right) + \left(1468 \right) \left(1406 \right) + \left(2176 \right) \left(2032\right) + \left( 3090 \right) \left(3078\right)}{\sqrt{\left( 1028 \right)^2 + \left( 1468\right)^2 + \left( 2176 \right)^2+ \left(3090\right)^2} \sqrt{\left( 969 \right)^2 + \left(1406 \right)^2 + \left( 2032\right)^2 + \left(3078 \right)^2 }}\\ &=0.0004176 \end{align*}
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\).
\begin{align*} &\text{BrayCurtis dissimilarity} \\ &= \frac{\left x_{\text{band 1}}^{(t)}  m_{\text{band 1}}\right + \left x_{\text{band 2}}^{(t)}  m_{\text{band 2}} \right + \dots + \left x_{\text{band p}}^{(t)}  m_{\text{band p}} \right }{ \left x_{\text{band 1}}^{(t)} + m_{\text{band 1}} \right + \left x_{\text{band 2}}^{(t)} + m_{\text{band 2}} \right + \dots + \left x_{\text{band p}}^{(t)} + m_{\text{band p}} \right} \\ &= \frac{\left x_{\text{nir1}}^{(t)}  m_{\text{nir1}}\right + \left x_{\text{green}}^{(t)}  m_{\text{green}} \right + \left x_{\text{red}}^{(t)}  m_{\text{red}} \right + \left x_{\text{blue}}^{(t)}  m_{\text{blue}} \right}{\left x_{\text{nir1}}^{(t)} + m_{\text{nir1}}\right + \left x_{\text{green}}^{(t)} + m_{\text{green}} \right + \left x_{\text{red}}^{(t)} + m_{\text{red}} \right + \left x_{\text{blue}}^{(t)} + m_{\text{blue}} \right} \\ &= \frac{\left 3090  3078\right + \left 1468  1406 \right + \left2176  2032 \right + \left 1028  969 \right}{\left 3090 + 3078\right + \left 1468 + 1406 \right + \left2176 + 2032 \right + \left 1028 + 969 \right}\\ &= 0.01817 \end{align*}