GeoMAD cloud-free composites

Date modified: March 2022

Header image - A section of Olifantsrivier in South Africa, imaged using the 2019 Sentinel-2 Annual GeoMAD. Contains modified Copernicus Sentinel data 2019, processed by Digital Earth Africa.

Service overview

Background

GeoMAD is the Digital Earth Africa (DE Africa) surface reflectance geomedian and triple Median Absolute Deviation data service. It is a cloud-free composite of satellite data compiled over specific timeframes. This service is ideal for:

  • Longer-term time series analysis

  • Cloudless 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 true-colour 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 higher-order 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 semiannual timeframe, for example six-month blocks, better captures seasonal variation within one year, but can also be used to compare equivalent periods from different years. Likewise, Landsat sensors allows full utility of the surface reflectance archive dating back to 1984, while more recent Sentinel-2 data provides higher-frequency flyovers and better resolution.

The Digital Earth Africa GeoMAD service currently provides annual and six-month semiannual datasets, with separate services for Landsat and Sentinel-2 sensors.

A Jupyter Notebook which demonstrates loading and using GeoMAD datasets in the Sandbox is also available.

Specifications

The GeoMAD service spans all of Africa and is currently available in two timeframes:

  • Annual: observations from one calendar year are summarised into one set of measurements

    • Sentinel-2 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

    • Sentinel-2 Semiannual GeoMAD

GeoMAD metadata can also be viewed on Digital Earth Africa’s Metadata Explorer.

Table 1: Sentinel-2 GeoMAD service specifications

Specification

Service name

Sentinel-2 Annual GeoMAD

Sentinel-2 Semiannual GeoMAD

Service status

Operational

Provisional

Number of bands

14

14

Cell size - X (metres)

10

10

Cell size - Y (metres)

10

10

Coordinate reference system (CRS)

EPSG: 6933

EPSG: 6933

Temporal resolution

Annual

6-monthly (Jan-Jun, Jul-Dec)

Temporal range

2017 – 2020

2017 – 2020

Parent dataset

Sentinel-2 Level-2A

Sentinel-2 Level-2A

Update frequency

Annual

TBD

Update latency

2 months from end of previous year

TBD

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 Level-2 Collection 2

Landsat Level-2 Collection 2

Landsat Level-2 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

S-2 GeoMAD data 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.

Sentinel-2, 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.

LS5-7 GeoMAD data extent over time LS5-7 GeoMAD count summary

Table 3.1: Sentinel-2 GeoMAD measurements

Band ID

Description

Value range

Data type

No data value

B02

Geomedian B02 (Blue)

1 - 10000

uint16

0

B03

Geomedian B03 (Green)

1 - 10000

uint16

0

B04

Geomedian B04 (Red)

1 - 10000

uint16

0

B05

Geomedian B05 (Red edge 1)

1 - 10000

uint16

0

B06

Geomedian B06 (Red edge 2)

1 - 10000

uint16

0

B07

Geomedian B07 (Red edge 3)

1 - 10000

uint16

0

B08

Geomedian B08 (Near infrared (NIR) 1)

1 - 10000

uint16

0

B8A

Geomedian B8A (NIR 2)

1 - 10000

uint16

0

B11

Geomedian B11 (Short-wave infrared (SWIR) 1)

1 - 10000

uint16

0

B12

Geomedian B12 (SWIR 2)

1 - 10000

uint16

0

SMAD

Spectral Median Absolute Deviation

0 - 1

float32

NaN

EMAD

Euclidean Median Absolute Deviation

0 - 31623

float32

NaN

BCMAD

Bray-Curtis Median Absolute Deviation

0 - 1

float32

NaN

COUNT

Number of clear observations

1 - 65535

uint16

0

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)

1 - 10000

uint16

0

SR_B3

Geomedian SR_B3 (Green)

1 - 10000

uint16

0

SR_B4

Geomedian SR_B4 (Red)

1 - 10000

uint16

0

SR_B5

Geomedian SR_B5 (NIR)

1 - 10000

uint16

0

SR_B6

Geomedian SR_B6 (SWIR 1)

1 - 10000

uint16

0

SR_B7

Geomedian SR_B7 (SWIR 2)

1 - 10000

uint16

0

SMAD

Spectral Median Absolute Deviation

0 - 1

float32

NaN

EMAD

Euclidean Median Absolute Deviation

0 - 31623

float32

NaN

BCMAD

Bray-Curtis Median Absolute Deviation

0 - 1

float32

NaN

COUNT

Number of clear observations

1 - 65535

uint16

0

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)

1 - 10000

uint16

0

SR_B2

Geomedian SR_B2 (Green)

1 - 10000

uint16

0

SR_B3

Geomedian SR_B3 (Red)

1 - 10000

uint16

0

SR_B4

Geomedian SR_B4 (NIR)

1 - 10000

uint16

0

SR_B5

Geomedian SR_B5 (SWIR 1)

1 - 10000

uint16

0

SR_B7

Geomedian SR_B7 (SWIR 2)

1 - 10000

uint16

0

SMAD

Spectral Median Absolute Deviation

0 - 1

float32

NaN

EMAD

Euclidean Median Absolute Deviation

0 - 31623

float32

NaN

BCMAD

Bray-Curtis Median Absolute Deviation

0 - 1

float32

NaN

COUNT

Number of clear observations

1 - 65535

uint16

0

Bands can be subdivided as follows:

  • Geomedian - 10 bands (Sentinel-2), 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 and 10000 to allow for more efficient data storage as unsigned 16-bit 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 Sentinel-2 Level-2A data. For example, the Annual GeoMAD band B02 contains the annual geomedian of the Sentinel-2 B02 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 Bray-Curtis 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 and 1. 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 Sentinel-2 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. Annual and Semiannual GeoMAD datasets for the period of 1984 – 2020 use hdstats version 0.2.

Media and example images

Image 1: Animations over northern Africa. 1984 – 2012 Landsat 5 & 7 Annual geomedian, true-colour (RGB).

Left: Agricultural area expansion near Sadat City, Egypt.

Right: Urbanisation in Tangier, Morocco.

Credit: Contains Landsat Level-2 Surface Reflectance Science Product courtesy of the U.S. Geological Survey, processed by Digital Earth Africa.

2d004efd64ab48aba95567c672b69590 3af79c7275ed45889dc3fc98460a7282

Image 2: Mangroves in Guinea-Bissau. 2019 Sentinel-2 Annual geomedian, true-colour (RGB).

Credit: Contains modified Copernicus Sentinel data 2019, processed by Digital Earth Africa.

Geomedian composite

Image 3: Croplands in South Africa. 2019 Sentinel-2 Annual triple MADs, plotted as RGB.

Credit: Contains modified Copernicus Sentinel data 2019, processed by Digital Earth Africa.

MADs

Image 4: Irrigated fields in Egypt. Jan-Jun 2020 (right), Jul-Dec 2020 (left). 2020 Sentinel-2 Semiannual geomedian, true-colour (RGB).

Credit: Contains modified Copernicus Sentinel data 2020, processed by Digital Earth Africa.

Semiannual geomedian

Image 5: The Okavango Delta, Botswana. 2020 Landsat 8 Annual triple MADs, plotted as RGB.

Credit: Contains Landsat Level-2 Surface Reflectance Science Product courtesy of the U.S. Geological Survey, processed by Digital Earth Africa.

Annual LS8 geomedian

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 no-data value. An example of this can be seen in Image 6 below where a shallow marine system contains no-data values in the GeoMAD because the NIR band values in the input data are persistently negative.

Image 6: Shallow marine system, Mauritius. 2021 Landsat 8 Annual GeoMAD (left), versus Sentinel-2 GeoMAD (right), plotted as RGB.

Maps portal link

9b162aad86874ce9a11f2398fffced2e

References

Roberts, D., Mueller, N., & Mcintyre, A. (2017). High-dimensional 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 high-dimensional statistics of time series. 8647-8650. https://doi.org/10.1109/IGARSS.2018.8518312.

License

CC BY Attribution 4.0 International License

Acknowledgements

The high-dimensional 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

arn:aws:s3:::deafrica-services

Product name

gm_s2_annual, gm_ls8_annual, gm_ls5_ls7_annual, gm_s2_semiannual

The bucket is in the AWS region af-south-1 (Cape Town). Additional region specifications can be applied as follows:

aws s3 ls --region=af-south-1 s3://deafrica-services/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

productname

Product name

gm_s2_annual

version

Product version

1.0.0

x

Tile number in the x direction.

x17

y

Tile number in the y direction.

y156

timeperiod

Year of data collection followed by period of time and time unit in the format YYYY--P<period><unit>. Time units are designated by capitalised letters, Y for years and M for months. Annual data is P1Y while 6-monthly data (such as the semiannuals) is P6M.

2019--P1Y

x_y_timeperiod_band.extension

File name. Combines x, y, timeperiod with band, using Band IDs, and file extensions. For most tiles, the file extension is .tif.

Tile numbering uses 96 km grid squares based on the default CRS for this product, EPSG:6933, with the origin set at the bottom-leftmost 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

DE Africa Services

Web Map Services (WMS) URL

https://ows.digitalearth.africa/wms?version=1.3.0

Web Coverage Service (WCS) URL

https://ows.digitalearth.africa/wcs?version=2.1.0

Layer name

gm_s2_annual, gm_ls8_annual, gm_ls5_ls7_annual, gm_s2_semiannual

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: Sentinel-2 GeoMAD ODC product and band names

ODC product name: gm_s2_annual, gm_s2_semiannual

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 case-sensitive.

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 scientifically-rigorous 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.

Geomedian composite

Figure 2 uses a three-dimensional 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 3-dimensional vector \(\mathbf{x}\). A single vector at timestep \(t\), \(\mathbf{x}^{(t)}\), looks like this:

\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 (straight-line 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 3-dimensional vector, with a value each for red, green, and blue.

\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 Sentinel-2 GeoMADs, the ten bands of the geomedian result in a 10-dimensional space - in the Landsat GeoMADs, the six geomedian bands give a 6-dimensional 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 Sentinel-2 annual services), while pixels in the centre of the path are only observed once per flyover, and have \(N\) closer to 70.

\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 3-D example. An example screenshot is shown below.

Figure 3: Geomedian visualisation widget.

Geomedian composite

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 second-order 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)

  • Bray-Curtis MAD, BCMAD (based on Bray-Curtis 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 straight-line distances between points. In our three-dimensional world, it may look like this:

Figure 4: Euclidean distance in three dimensions.

Euclidean

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.

Euclidean

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 Sentinel-2 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.

Cosine distance

In a general sense, cosine distance is related to the angle between the two points \(\theta\), while Euclidean distance is related to the straight-line 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 01.

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.

Bray-Curtis MAD (BCMAD)

The Bray-Curtis MAD (BCMAD) is calculated from the Bray-Curtis dissimilarity. The Bray-Curtis dissimilarity emphasises differences in each band between the measurement and the geomedian.

For a single band of satellite data, the Bray-Curtis 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{Single-band Bray-Curtis 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 Bray-Curtis 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 Bray-Curtis 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 Bray-Curtis 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 Bray-Curtis dissimilarity

Using the example multispectral data from Table A.1, we can manually calculate the value of the Bray-Curtis dissimilarity for timestep \(t\).

\begin{align*} &\text{Bray-Curtis 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| + \left|2176 - 2032 \right| + \left| 1028 - 969 \right|}{\left| 3090 + 3078\right| + \left| 1468 + 1406 \right| + \left|2176 + 2032 \right| + \left| 1028 + 969 \right|}\\ &= 0.01817 \end{align*}