diff --git a/tutorials/imad-tutorial-pt1/index.ipynb b/tutorials/imad-tutorial-pt1/index.ipynb new file mode 100644 index 000000000..53421a5fa --- /dev/null +++ b/tutorials/imad-tutorial-pt1/index.ipynb @@ -0,0 +1,835 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "8kdsGkYJXXKc" + }, + "outputs": [], + "source": [ + "#@title Copyright 2023 The Earth Engine Community Authors { display-mode: \"form\" }\n", + "#\n", + "# Licensed under the Apache License, Version 2.0 (the \"License\");\n", + "# you may not use this file except in compliance with the License.\n", + "# You may obtain a copy of the License at\n", + "#\n", + "# https://www.apache.org/licenses/LICENSE-2.0\n", + "#\n", + "# Unless required by applicable law or agreed to in writing, software\n", + "# distributed under the License is distributed on an \"AS IS\" BASIS,\n", + "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n", + "# See the License for the specific language governing permissions and\n", + "# limitations under the License." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "pJrMq0aNrF-u" + }, + "source": [ + "# Change Detection on Google Earth Engine\n", + "# Part 1 The MAD Transformation\n", + "Authors: mortcanty" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "4xFjivyNFe7U" + }, + "source": [ + "## Context\n", + "\n", + "With the great tools available in Google Earth Engine (GEE) it is easy to generate satellite image differences and time series animations which visualize changes on the Earth surface, sometimes quite dramatically. Such representations are very good as eye-openers to give a conceptual, overall impression of change. They may, however, be perceived differently from person to person and not allow more than subjective interpretation. With more objective and rigorous statistical methods one can often extract more and better change information on both pixel and patch/field levels.\n", + "\n", + "The Multivariate Alteration Detection (MAD) transformation was proposed and developed some time ago by [Allan Nielsen and his coworkers](https://www2.imm.dtu.dk/pubdb/pubs/1220-full.html) at the Technical University of Denmark and later [extended to an iteratively re-weighted version](https://www2.imm.dtu.dk/pubdb/pubs/4695-full.html) called IR-MAD or iMAD. The iMAD algorithm has since found widespread application for the generation of change information from visual/infrared imagery as well as for performing [relative radiometric normalization](http://www2.imm.dtu.dk/pubdb/pubs/5362-full.html) of image pairs.\n", + "\n", + "This tutorial is intended to familiarize GEE users with the iMAD method so that they may use it with confidence in their research wherever they think it might be appropriate. It is based on material taken from Chapter 9 of [Canty (2019)](https://www.taylorfrancis.com/books/image-analysis-classification-change-detection-remote-sensing-morton-john-canty/10.1201/9780429464348). In preparing the tutorial we have tried to take as much advantage as possible of the GEE platform to illustrate and demonstrate the theory interactively.\n", + "\n", + "## Outline\n", + "\n", + "The tutorial is split into three parts, of which this is the first:\n", + "\n", + " 1. The MAD Transformation\n", + " 2. The iMAD Algorithm\n", + " 3. Radiometric Normalization and Harmonization\n", + "\n", + "In Part 3, a convenient Jupyter widget interface for running the iMAD algorithm from a Docker container is also described.\n", + "\n", + "## Prerequisites\n", + "Little prior knowledge is required to follow the discussion, apart from some familiarity with the GEE API and a basic, intuitive understanding of the simple statistical concepts of probability distributions and their moments, in particular, mean, variance (var) and covariance (cov). Since we'll be working with multispectral imagery, in the course of the tutorial we'll frequently mention the so-called *dispersion matrix* or *variance-covariance matrix*, or simply *covariance matrix* for short. The covariance matrix of a set of measured variables $X_i,\\ i=1\\dots N$, such as the band-wise values of a multispectral image pixel, is defined as\n", + "\n", + "$$\n", + "\\Sigma_X\\ =\\ \\begin{pmatrix} \\sigma^2_1 \u0026 \\sigma_{12} \u0026 \\dots \u0026 \\sigma_{1N} \\cr\n", + " \\sigma_{21} \u0026 \\sigma^2_2\u0026 \\dots \u0026 \\sigma_{2N} \\cr\n", + " \\vdots \u0026 \\vdots \u0026 \\vdots\u0026 \\vdots \\cr\n", + " \\sigma_{N1} \u0026 \\sigma_{N2} \u0026 \\dots \u0026 \\sigma^2_N \\end{pmatrix},\n", + "$$\n", + "\n", + "where\n", + "\n", + "$$\n", + "\\sigma_{ij} = {\\rm cov}(X_i,X_j) = \\langle (X_i - \\bar X_i)(X_j - \\bar X_j)\\rangle, \\quad\n", + "\\sigma^2_i = {\\rm cov}(X_i,X_i) = {\\rm var}(X_i) = \\langle (X_i-\\bar X_i)^2\\rangle, \\quad \\ i,j = 1\\dots N.\n", + "$$\n", + "\n", + "Here, $\\bar X_i$ denotes mean value and $\\langle\\dots\\rangle$ is ensemble average. The covariance matrix is symmetric about its principal diagonal. The *correlation* between $X_i$ and $X_j$ is defined as the normalized covariance\n", + "\n", + "$$\n", + "\\rho_{ij} = {{\\rm cov}(X_i,X_j)\\over \\sqrt{{\\rm var}(X_i)}\\sqrt{{\\rm var}(X_j)}} = {\\sigma_{ij} \\over \\sigma_i\\sigma_j}.\n", + "$$\n", + "\n", + "Since we are in a Colab environment we will use the Python API, but it is so similar to the JavaScript version that this should pose no problems to anyone not acquainted with it.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "cClM-tAWSnuS" + }, + "source": [ + "## Preliminaries\n", + "\n", + "The cells below execute the necessary formalities for accessing Earth Engine from this Colab notebook. They also import the various Python add-ons needed, including the [geemap](https://geemap.org/) interactive map package, and define a few helper routines for displaying results." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "ZdyBUBN4F1IA" + }, + "outputs": [], + "source": [ + "import ee\n", + "\n", + "ee.Authenticate(auth_mode='notebook')\n", + "ee.Initialize()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "NqcghIcEIWnZ" + }, + "outputs": [], + "source": [ + "# Import other packages used in the tutorial\n", + "%matplotlib inline\n", + "import geemap\n", + "import numpy as np\n", + "import random, time\n", + "import matplotlib.pyplot as plt\n", + "from scipy.stats import norm, chi2\n", + "\n", + "from pprint import pprint # for pretty printing" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "6z-Jcr3n8MhD" + }, + "outputs": [], + "source": [ + "# Truncate a 1-D array to dec decimal places\n", + "def trunc(values, dec = 3):\n", + " return np.trunc(values*10**dec)/(10**dec)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "JXWO7Qug8MhE" + }, + "outputs": [], + "source": [ + "# Display an image in a one percent linear stretch\n", + "def display_ls(image, map, name, centered = False):\n", + " bns = image.bandNames().length().getInfo()\n", + " if bns == 3:\n", + " image = image.rename('B1', 'B2', 'B3')\n", + " pb_99 = ['B1_p99', 'B2_p99', 'B3_p99']\n", + " pb_1 = ['B1_p1', 'B2_p1', 'B3_p1']\n", + " img = ee.Image.rgb(image.select('B1'), image.select('B2'), image.select('B3'))\n", + " else:\n", + " image = image.rename('B1')\n", + " pb_99 = ['B1_p99']\n", + " pb_1 = ['B1_p1']\n", + " img = image.select('B1')\n", + " percentiles = image.reduceRegion(ee.Reducer.percentile([1, 99]), maxPixels=1e11)\n", + " mx = percentiles.values(pb_99)\n", + " if centered:\n", + " mn = ee.Array(mx).multiply(-1).toList()\n", + " else:\n", + " mn = percentiles.values(pb_1)\n", + " map.addLayer(img, {'min': mn, 'max': mx}, name)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "bn5m_-sOpkDj" + }, + "source": [ + "## Simple Differences\n", + "Let's begin by considering two $N$-band optical/infrared images of the same scene acquired with the same sensor at two different times, between which ground reflectance changes have occurred at\n", + "some locations but not everywhere. To see those changes we can \"flick back and forth\" between gray scale or RGB representations of the images on a computer display. Alternatively we can simply subtract one image from the other and, provided the intensities at the two time points have been calibrated or normalized in some way, examine the difference. Symbolically the difference image is\n", + "\n", + "$$\n", + "Z = X-Y\n", + "$$\n", + "\n", + "where\n", + "$$\n", + "X= \\begin{pmatrix}X_1\\cr X_2\\cr\\vdots\\cr X_N\\end{pmatrix},\\quad Y= \\begin{pmatrix}Y_1\\cr Y_2\\cr\\vdots\\cr Y_N\\end{pmatrix}\n", + "$$\n", + "\n", + "are 'typical' pixel vectors (more formally: _random vectors_) representing the first and second image, respectively." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "KX-pPWU9wO5I" + }, + "source": [ + "### Landkreis Olpe\n", + "\n", + "For example, here is an area of interest (aoi) covering a heavily forested administrative district (Landkreis Olpe) in North Rhine Westphalia, Germany. Due to severe drought in recent years large swaths of shallow-root coniferous trees have died and been cleared away, leaving deciduous trees for the most part untouched." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "cellView": "form", + "id": "Xs39zZmNSdru" + }, + "outputs": [], + "source": [ + "aoi = ee.FeatureCollection('projects/google/imad_tutorial_aoi').geometry()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "rwCBEe49d-xx" + }, + "source": [ + "We first collect two Sentinel-2 scenes which bracket some of the recent clear cutting (June, 2021 and June, 2022) and print out their timestamps. For demonstration purposes we choose top of atmosphere reflectance since the MAD transformation does not require absolute surface reflectances." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "JzDg-puTSybq" + }, + "outputs": [], + "source": [ + "def collect(aoi, t1a ,t1b, t2a, t2b):\n", + " try:\n", + " im1 = ee.Image( ee.ImageCollection(\"COPERNICUS/S2_HARMONIZED\")\n", + " .filterBounds(aoi)\n", + " .filterDate(ee.Date(t1a), ee.Date(t1b))\n", + " .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo'))\n", + " .sort('CLOUDY_PIXEL_PERCENTAGE')\n", + " .first()\n", + " .clip(aoi) )\n", + " im2 = ee.Image( ee.ImageCollection(\"COPERNICUS/S2_HARMONIZED\")\n", + " .filterBounds(aoi)\n", + " .filterDate(ee.Date(t2a), ee.Date(t2b))\n", + " .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo'))\n", + " .sort('CLOUDY_PIXEL_PERCENTAGE')\n", + " .first()\n", + " .clip(aoi) )\n", + " timestamp = im1.date().format('E MMM dd HH:mm:ss YYYY')\n", + " print(timestamp.getInfo())\n", + " timestamp = im2.date().format('E MMM dd HH:mm:ss YYYY')\n", + " print(timestamp.getInfo())\n", + " return (im1, im2)\n", + " except Exception as e:\n", + " print('Error: %s'%e)\n", + "\n", + "im1, im2 = collect(aoi, '2021-06-01', '2021-06-30', '2022-06-01', '2022-06-30')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "zekLrceuaNtK" + }, + "source": [ + "In order to view the images, we instantiate an interactive map located at the center of the aoi. Since the map requires lat/long coordinates we have to reverse the long/lat order returned from the *coordinates()* method of the *ee.Geometry.centroid()* object:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "kc_pwScrphK1" + }, + "outputs": [], + "source": [ + "# Interactive map\n", + "M1 = geemap.Map()\n", + "M1.centerObject(aoi, 11)\n", + "\n", + "M1" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "vYKXgtJ2FdYP" + }, + "source": [ + "Then we display RGB composites of the first 3 (visual) bands of the Sentinel images as well as their difference. We use the previously defined function *display_ls()* to display the images in a one percent linear stretch." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "W8_al6RlBpNq" + }, + "outputs": [], + "source": [ + "visirbands = ['B2','B3','B4','B8','B11','B12']\n", + "visbands = ['B2','B3','B4']\n", + "\n", + "diff = im1.subtract(im2).select(visbands)\n", + "display_ls(im1.select(visbands), M1, 'Image 1')\n", + "display_ls(im2.select(visbands), M1, 'Image 2')\n", + "display_ls(diff, M1, 'Difference')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "h2d-vg-ezwpe" + }, + "source": [ + "Small intensity differences (the gray-ish background in the difference image above) indicate no change, large\n", + "positive (bright) or negative (dark) colors indicate change. The most obvious changes are in fact due to the clear cutting. Decision thresholds can be set to define significant changes. The thresholds are usually expressed in terms of standard deviations from the mean difference value, which is taken to correspond to no change. The per-band variances of the difference images are simply\n", + "\n", + "$$\n", + "{\\rm var}(Z_i) = {\\rm var}(X_i - Y_i) = {\\rm var}(X_i) + {\\rm var}(Y_i) - 2\\cdot{\\rm cov}(X_i,Y_i),\\quad i=1\\dots N,\n", + "$$\n", + "\n", + "or, if the bands are uncorrelated, about twice as noisy as the individual image bands. Normally $X_i$ and $Y_i$ are positively correlated (${\\rm cov}(X_i,Y_i) \u003e 0$) so the variance of $Z_i$ is somewhat smaller. When the difference signatures in the spectral channels are combined so as to try to characterize the nature of the changes that have taken place, one speaks of _spectral change vector analysis_.\n", + "\n", + "While the recent clear cuts are fairly easily identified as (some but not all of) the dark areas in the simple difference image above, it is possible to derive a better and more informative change map. This involves taking greater advantage of the statistical properties of the images." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "rhMR6KLhkar5" + }, + "source": [ + "## The MAD Transformation\n", + "Let's make a linear combination of the intensities for all $N$ bands in the first image $X$, thus creating a scalar image\n", + "\n", + "$$\n", + "U = a_1X_1 + a_2X_2 + \\dots + a_NX_N = {\\bf a}^\\top X. \\tag{1}\n", + "$$\n", + "\n", + "The symbol ${\\bf a}^\\top$ denotes the transpose of the column vector ${\\bf a}=\\begin{pmatrix}a_1\\cr a_2\\cr\\vdots\\cr a_N\\end{pmatrix}$, in other words the row vector ${\\bf a}^\\top = (a_1,a_2 \\dots a_N)$.\n", + "\n", + "\n", + "The expression ${\\bf a}^\\top X$ is called an _inner vector product_. The vector of coefficients ${\\bf a}$ is as yet unspecified.\n", + "We'll do the same for the second image $Y$, forming the linear combination\n", + "\n", + "$$\n", + "V = b_1Y_1 + b_2Y_2 + \\dots + b_NY_N = {\\bf b}^\\top Y, \\tag{2}\n", + "$$\n", + "\n", + "and then look at the scalar difference $U-V$. Change information\n", + "is now contained, for the time being, in a single image rather than distributed among all $N$ bands.\n", + "\n", + "We have of course to choose the coefficient vectors ${\\bf a}$ and\n", + "${\\bf b}$ in some suitable way. In [Nielsen et al. (1998)](https://www2.imm.dtu.dk/pubdb/pubs/1220-full.html) it is suggested that they be determined by making the transformed images $U$ and $V$ _as similar as they can be_ before taking their difference. This sounds at first counter-intuitive: Why make the images similar when we want to see the dissimilarities? The crux is that, in making them match one another by taking suitable linear combinations, we ensure that pixels in the transformed images $U$ and $V$ at which *indeed no change has taken place* will be as similar as possible before they are compared (subtracted). The genuine changes will then, it is hoped, be all the more apparent in the difference image.\n", + "\n", + "This process of \"similar making\" can be accomplished in statistics with standard [*Canonical Correlation Analysis* (CCA)](https://en.wikipedia.org/wiki/Canonical_correlation), a technique first described by Harold Hotelling in 1936. Canonical Correlation Analysis of the images represented by Equations (1) and (2) *maximizes the correlation*\n", + "\n", + "$$\n", + "\\rho = {{\\rm cov}(U,V)\\over \\sqrt{{\\rm var}(U)}\\sqrt{{\\rm var}(V)}}, \\tag{3}\n", + "$$\n", + "\n", + "between the random variables $U$ and $V$, which is just what we want.\n", + "\n", + "One technical point: Arbitrary multiples of $U$ and $V$ would clearly have the same correlation (the multiples will cancel off in numerator and denominator), so a constraint must be chosen. A convenient one is\n", + "\n", + "$$\n", + "{\\rm var}(U)={\\rm var}(V)=1. \\tag{4}\n", + "$$\n", + "\n", + "For those more versed in linear algebra the details of CCA are given in [Canty (2019)](https://www.taylorfrancis.com/books/image-analysis-classification-change-detection-remote-sensing-morton-john-canty/10.1201/9780429464348) beginning on page 385. It all boils down to solving two so-called *generalized eigenvalue problems* for determining the transformation vectors ${\\bf a}$ and ${\\bf b}$ in Equations (1) and (2), respectively. There is not just one solution but rather $N$ solutions. They consist of the $N$ pairs of _eigenvectors_ $({\\bf a}^i, {\\bf b}^i)$,\n", + "\n", + "$$\n", + "{\\bf a}^i=\\begin{pmatrix}a_1\\cr a_2\\cr\\vdots\\cr a_N \\end{pmatrix}_i,\\quad {\\bf b}^i=\\begin{pmatrix}b_1\\cr b_2\\cr\\vdots\\cr b_N \\end{pmatrix}_i, \\quad i=1\\dots N\n", + "$$\n", + "\n", + "and, correspondingly, the $N$ pairs of *canonical variates* $(U_i, V_i),\\ i= 1\\dots N$. (Readers familiar with Principal Components Analysis (PCA) will recognize the similarity. When applying PCA to an N-band image, an ordinary (rather than a generalized) eigenvalue problem is solved to maximize variance, resulting in $N$ principal component bands.)\n", + "\n", + "Unlike the bands of the original images $X$ and $Y$, which are ordered by spectral wavelength, the canonical variates are ordered by similarity or correlation. The difference images\n", + "\n", + "$$\n", + "M_i = U_i - V_i = ({\\bf a}^i)^\\top X - ({\\bf b}^i)^\\top Y,\\quad i=1\\dots N, \\tag{5}\n", + "$$\n", + "\n", + "contain the change information and are called the *MAD variates*.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "ePBGinhxsSg_" + }, + "source": [ + "### Auxiliary functions" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "o0j6--qcsSg_" + }, + "source": [ + "To run the MAD transformation on GEE image pairs we need some helper routines. These are coded below. The first, called *covarw()*, returns the weighted covariance matrix of an $N$-band image as well as a weighted, centered (mean-zero) version of the image. Why it is important to weight the pixels before centering or calculating the covariance matrix will become apparent in Part 2 of the tutorial." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "xJXqZ3ZJsSg_" + }, + "outputs": [], + "source": [ + "def covarw(image, weights = None, scale = 20, maxPixels = 1e10):\n", + " '''Return the weighted centered image and its weighted covariance matrix'''\n", + " try:\n", + " geometry = image.geometry()\n", + " bandNames = image.bandNames()\n", + " N = bandNames.length()\n", + " if weights is None:\n", + " weights = image.constant(1)\n", + " weightsImage = image.multiply(ee.Image.constant(0)).add(weights)\n", + " means = image.addBands(weightsImage) \\\n", + " .reduceRegion(ee.Reducer.mean().repeat(N).splitWeights(),\n", + " scale = scale,\n", + " maxPixels = maxPixels) \\\n", + " .toArray() \\\n", + " .project([1])\n", + " centered = image.toArray().subtract(means)\n", + " B1 = centered.bandNames().get(0)\n", + " b1 = weights.bandNames().get(0)\n", + " nPixels = ee.Number(centered.reduceRegion(ee.Reducer.count(),\n", + " scale = scale,\n", + " maxPixels = maxPixels).get(B1))\n", + " sumWeights = ee.Number(weights.reduceRegion(ee.Reducer.sum(),\n", + " geometry = geometry,\n", + " scale = scale,\n", + " maxPixels = maxPixels).get(b1))\n", + " covw = centered.multiply(weights.sqrt()) \\\n", + " .toArray() \\\n", + " .reduceRegion(ee.Reducer.centeredCovariance(),\n", + " geometry = geometry,\n", + " scale = scale,\n", + " maxPixels = maxPixels) \\\n", + " .get('array')\n", + " covw = ee.Array(covw).multiply(nPixels).divide(sumWeights)\n", + " return (centered.arrayFlatten([bandNames]), covw)\n", + " except Exception as e:\n", + " print('Error: %s'%e)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "bEtccjvpsShA" + }, + "source": [ + "The second routine, called _corr()_, transforms a covariance matrix to the equivalent correlation matrix by dividing by the square roots of the variances, see Equation (3)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "olCi6RvJsShA" + }, + "outputs": [], + "source": [ + "def corr(cov):\n", + " '''Transform covariance matrix to correlation matrix'''\n", + " # diagonal matrix of inverse sigmas\n", + " sInv = cov.matrixDiagonal().sqrt().matrixToDiag().matrixInverse()\n", + " # transform\n", + " corr = sInv.matrixMultiply(cov).matrixMultiply(sInv).getInfo()\n", + " # truncate\n", + " return [list(map(trunc, corr[i])) for i in range(len(corr))]" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "G34S_Z2gQAwp" + }, + "source": [ + "For example the code below calculates the (unweighted) covariance matrix of the visual and infrared bands of the image _im1_ and prints out the corresponding correlation matrix." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "TmFXfPhYDSNi" + }, + "outputs": [], + "source": [ + "_, cov = covarw(im1.select(visirbands))\n", + "pprint(corr(cov))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "oWhTaIRXsShB" + }, + "source": [ + "Note that the spectral bands of the image are generally highly and positively correlated with one another.\n", + "\n", + "By stacking the images one atop the other we can use the helper functions to display the between image correlations\n", + "\n", + "$$\n", + "\\rm{corr}(X_i,Y_i), \\quad i= 1\\dots N.\n", + "$$\n", + "\n", + "They can be found in the diagonal of the upper left $6\\times 6$ submatrix of the correlation matrix for the stacked image." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "vXJkxrZmsShB" + }, + "outputs": [], + "source": [ + "im12 = im1.select(visirbands).addBands(im2.select(visirbands))\n", + "_, covar = covarw(im12)\n", + "correl = np.array(corr(covar))\n", + "print(np.diag(correl[:6, 6:]))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "iJC-EC_PAXo9" + }, + "source": [ + "It is these between image correlations that the MAD algorithm tries to maximize.\n", + "\n", + "The third and last auxiliary routine, *geneiv()*, is the core of CCA and the MAD transformation. It solves the *generalized eigenproblem*,\n", + "\n", + "$$\n", + "CX = \\lambda BX\n", + "$$\n", + "\n", + "for two $N\\times N$ matrices $C$ and $B$, returning the $N$ solutions, or eigenvectors, $X_i$ and the $N$ eigenvalues $\\lambda_i, \\ i=1\\dots N$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "cellView": "form", + "id": "lK3AMH709Sg5" + }, + "outputs": [], + "source": [ + "def geneiv(C,B):\n", + " '''Return the eigenvalues and eigenvectors of the generalized eigenproblem\n", + " C*X = lambda*B*X'''\n", + " try:\n", + " C = ee.Array(C)\n", + " B = ee.Array(B)\n", + " # Li = choldc(B)^-1\n", + " Li = ee.Array(B.matrixCholeskyDecomposition().get('L')).matrixInverse()\n", + " # solve symmetric, ordinary eigenproblem Li*C*Li^T*x = lambda*x\n", + " Xa = Li.matrixMultiply(C) \\\n", + " .matrixMultiply(Li.matrixTranspose()) \\\n", + " .eigen()\n", + " # eigenvalues as a row vector\n", + " lambdas = Xa.slice(1, 0, 1).matrixTranspose()\n", + " # eigenvectors as columns\n", + " X = Xa.slice(1, 1).matrixTranspose()\n", + " # generalized eigenvectors as columns, Li^T*X\n", + " eigenvecs = Li.matrixTranspose().matrixMultiply(X)\n", + " return (lambdas, eigenvecs)\n", + " except Exception as e:\n", + " print('Error: %s'%e)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "qksWAxsrIV4g" + }, + "source": [ + "The next cell codes the MAD transformation itself in the funcion *mad_run()*, taking as input two multiband images and returning the _canonical variates_\n", + "\n", + "$$\n", + "U_i, \\ V_i, \\quad i=1\\dots N,\n", + "$$\n", + "\n", + "the _MAD variates_\n", + "\n", + "$$\n", + "M_i = U_i - V_i, \\quad i=1\\dots N,\n", + "$$\n", + "\n", + "as well as the sum of the squares of the standardized MAD variates,\n", + "\n", + "$$\n", + "Z = \\sum_{i=1}^N\\left({M_i\\over \\sigma_{M_i}}\\right)^2.\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "-4I6Ig_zFICj" + }, + "outputs": [], + "source": [ + "#@title The MAD transformation\n", + "def mad_run(image1, image2, scale = 20):\n", + " '''The MAD transformation of two multiband images'''\n", + " try:\n", + " image = image1.addBands(image2)\n", + " region = image.geometry()\n", + " nBands = image.bandNames().length().divide(2)\n", + " centeredImage,covarArray = covarw(image,scale=scale)\n", + " bNames = centeredImage.bandNames()\n", + " bNames1 = bNames.slice(0,nBands)\n", + " bNames2 = bNames.slice(nBands)\n", + " centeredImage1 = centeredImage.select(bNames1)\n", + " centeredImage2 = centeredImage.select(bNames2)\n", + " s11 = covarArray.slice(0,0,nBands).slice(1,0,nBands)\n", + " s22 = covarArray.slice(0,nBands).slice(1,nBands)\n", + " s12 = covarArray.slice(0,0,nBands).slice(1,nBands)\n", + " s21 = covarArray.slice(0,nBands).slice(1,0,nBands)\n", + " c1 = s12.matrixMultiply(s22.matrixInverse()).matrixMultiply(s21)\n", + " b1 = s11\n", + " c2 = s21.matrixMultiply(s11.matrixInverse()).matrixMultiply(s12)\n", + " b2 = s22\n", + " # solution of generalized eigenproblems\n", + " lambdas, A = geneiv(c1,b1)\n", + " _, B = geneiv(c2,b2)\n", + " rhos = lambdas.sqrt().project(ee.List([1]))\n", + " # MAD variances\n", + " sigma2s = rhos.subtract(1).multiply(-2).toList()\n", + " sigma2s = ee.Image.constant(sigma2s)\n", + " # ensure sum of correlations between X and U is positive\n", + " tmp = s11.matrixDiagonal().sqrt()\n", + " ones = tmp.multiply(0).add(1)\n", + " tmp = ones.divide(tmp).matrixToDiag()\n", + " s = tmp.matrixMultiply(s11).matrixMultiply(A).reduce(ee.Reducer.sum(),[0]).transpose()\n", + " A = A.matrixMultiply(s.divide(s.abs()).matrixToDiag())\n", + " # ensure positive correlation between U and V\n", + " tmp = A.transpose().matrixMultiply(s12).matrixMultiply(B).matrixDiagonal()\n", + " tmp = tmp.divide(tmp.abs()).matrixToDiag()\n", + " B = B.matrixMultiply(tmp)\n", + " # canonical and MAD variates\n", + " centeredImage1Array = centeredImage1.toArray().toArray(1)\n", + " centeredImage2Array = centeredImage2.toArray().toArray(1)\n", + " U = ee.Image(A.transpose()).matrixMultiply(centeredImage1Array) \\\n", + " .arrayProject([0]) \\\n", + " .arrayFlatten([bNames2])\n", + " V = ee.Image(B.transpose()).matrixMultiply(centeredImage2Array) \\\n", + " .arrayProject([0]) \\\n", + " .arrayFlatten([bNames2])\n", + " MAD = U.subtract(V)\n", + " # chi square image\n", + " Z = MAD.pow(2) \\\n", + " .divide(sigma2s) \\\n", + " .reduce(ee.Reducer.sum()) \\\n", + " .clip(region)\n", + " return (U, V, MAD, Z)\n", + " except Exception as e:\n", + " print('Error: %s'%e)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "jLuzqZ4lR1Ck" + }, + "source": [ + "After setting up the MAD transformation on the six visual and infrared bands of the Landkreis Olpe images," + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "zTK-8PV5NvtO", + "scrolled": false + }, + "outputs": [], + "source": [ + "U, V, MAD, Z = mad_run(im1.select(visirbands), im2.select(visirbands))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "vbeJ6WF_sShC" + }, + "source": [ + "we display an RGB composite of the first 3 MAD variates, together with the original images and their simple difference:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "zNAGaOr8sShC" + }, + "outputs": [], + "source": [ + "M2 = geemap.Map()\n", + "M2.centerObject(aoi, 11)\n", + "display_ls(im1.select(visbands), M2, 'Image 1')\n", + "display_ls(im2.select(visbands), M2, 'Image 2')\n", + "display_ls(diff, M2, 'Difference')\n", + "display_ls(MAD.select(0, 1, 2), M2, 'MAD Image', True)\n", + "\n", + "M2" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "ALAZOS0EdM24" + }, + "source": [ + "The richness in colors in the MAD image is a consequence of the decoupling (orthogonalization) of the MAD components, as will be explained below.\n", + "\n", + "Pretty colors notwithstanding, when compared with the simple difference image the result is _mittelprächtig_ (German slang for \"not so hot\"), and certainly no easier to interpret. However we're not finished.\n", + "\n", + "The canonical variates have very nice properties indeed. They are _all mutually uncorrelated_ except for the pairs $(U_i,V_i)$, and these are ordered by decreasing positive correlation:\n", + "\n", + "$$\n", + "\\rho_i = {\\rm cov}(U_i, V_i),\\quad i=1\\dots N, \\tag{6}\n", + "$$\n", + "\n", + "with $\\rho_1 \\ge \\rho_2 \\ge\\dots \\ge\\rho_N$.\n", + "\n", + "Stacking the images $U$ and $V$, these facts can be read from the $12\\times 12$ correlation matrix:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "fwQjSqbPdPUQ" + }, + "outputs": [], + "source": [ + "_, covar = covarw(U.addBands(V))\n", + "correl = np.array(corr(covar))\n", + "pprint(correl)\n", + "print('rho =', np.diag(correl[:6,6:]))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "aIWj8a7j8MhO" + }, + "source": [ + "The MAD variates themselves are consequently also mutually uncorrelated, their covariances being given by\n", + "\n", + "$$\n", + "{\\rm cov}(M_i,M_j) = {\\rm cov}(U_i-V_i,U_j-V_j)= 0,\\quad i\\ne j=1\\dots N, \\tag{7}\n", + "$$\n", + "\n", + "and their variances by\n", + "\n", + "$$\n", + "\\sigma_{M_i}^2={\\rm var}(U_i-V_i)= {\\rm var}(U_i)+{\\rm var}(V_i) -2{\\rm cov}(U_i,V_i) = 2(1-\\rho_i),\\quad i=1\\dots N, \\tag{8}\n", + "$$\n", + "\n", + "where the last equality follows from the constraint Equation (4). Let's check." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "kXOtyyIThIOw" + }, + "outputs": [], + "source": [ + "# display MAD covariance matrix\n", + "_, covar = covarw(MAD)\n", + "covar = covar.getInfo()\n", + "[list(map(trunc,covar[i])) for i in range(len(covar))]" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "HOCSrVuo6jcL" + }, + "source": [ + "\n", + "\n", + "The first MAD variate has minimum variance in its pixel\n", + "intensities. The second MAD variate has minimum variance subject to the condition that its pixel intensities are statistically\n", + "uncorrelated with those in the first variate; the third has\n", + "minimum spread subject to being uncorrelated with the first two,\n", + "and so on. Depending on the type of change present, any of the\n", + "components may exhibit significant change information, although it tends to be concentrated in the first few MAD variates. One of the nicest aspects of the MAD transformation is that it can sort different categories of change into different, uncorrelated image bands.\n", + "\n", + "However, the result, as we saw, needs some improvement. This is the subject of Part 2." + ] + } + ], + "metadata": { + "colab": { + "provenance": [] + }, + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/tutorials/imad-tutorial-pt1/index.md b/tutorials/imad-tutorial-pt1/index.md deleted file mode 100644 index 60582b668..000000000 --- a/tutorials/imad-tutorial-pt1/index.md +++ /dev/null @@ -1,549 +0,0 @@ ---- -title: Change Detection in GEE - Part 1 The MAD Transformation -description: Iteratively re-weighted Multivariate Alteration Detection. -author: mortcanty -tags: change-detection, imad -date_published: 2023-02-07 ---- - - -{% include "earth-engine/_widgets/_notebook_buttons.html" with github_path="mortcanty/eeimad/blob/main/src/gee_change_detection_pt_1.ipynb" %} - -## Context - -With the great tools available in Google Earth Engine (GEE) it is easy to generate satellite image differences and time series animations which visualize changes on the Earth surface, sometimes quite dramatically. Such representations are very good as eye-openers to give a conceptual, overall impression of change. They may, however, be perceived differently from person to person and not allow more than subjective interpretation. With more objective and rigorous statistical methods one can often extract more and better change information on both pixel and patch/field levels. - -The Multivariate Alteration Detection (MAD) transformation was proposed and developed some time ago by [Allan Nielsen and his coworkers](https://www2.imm.dtu.dk/pubdb/pubs/1220-full.html) at the Technical University of Denmark and later [extended to an iteratively re-weighted version](https://www2.imm.dtu.dk/pubdb/pubs/4695-full.html) called IR-MAD or iMAD. The iMAD algorithm has since found widespread application for the generation of change information from visual/infrared imagery as well as for performing [relative radiometric normalization](http://www2.imm.dtu.dk/pubdb/pubs/5362-full.html) of image pairs. - -This tutorial is intended to familiarize GEE users with the iMAD method so that they may use it with confidence in their research wherever they think it might be appropriate. It is based on material taken from Chapter 9 of [Canty (2019)](https://www.taylorfrancis.com/books/image-analysis-classification-change-detection-remote-sensing-morton-john-canty/10.1201/9780429464348). In preparing the tutorial we have tried to take as much advantage as possible of the GEE platform to illustrate and demonstrate the theory interactively. - -## Outline - -The tutorial is split into three parts, of which this is the first: - - 1. The MAD Transformation - 2. The iMAD Algorithm - 3. Radiometric Normalization and Harmonization - -In Part 3, a convenient Jupyter widget interface for running the iMAD algorithm from a Docker container is also described. - -## Prerequisites -Little prior knowledge is required to follow the discussion, apart from some familiarity with the GEE API and a basic, intuitive understanding of the simple statistical concepts of probability distributions and their moments, in particular, mean, variance (var) and covariance (cov). Since we'll be working with multispectral imagery, in the course of the tutorial we'll frequently mention the so-called *dispersion matrix* or *variance-covariance matrix*, or simply *covariance matrix* for short. The covariance matrix of a set of measured variables $X_i,\ i=1\dots N$, such as the band-wise values of a multispectral image pixel, is defined as - -$$ -\Sigma_X\ =\ \begin{pmatrix} \sigma^2_1 & \sigma_{12} & \dots & \sigma_{1N} \cr - \sigma_{21} & \sigma^2_2& \dots & \sigma_{2N} \cr - \vdots & \vdots & \vdots& \vdots \cr - \sigma_{N1} & \sigma_{N2} & \dots & \sigma^2_N \end{pmatrix}, -$$ - -where - -$$ -\sigma_{ij} = {\rm cov}(X_i,X_j) = \langle (X_i - \bar X_i)(X_j - \bar X_j)\rangle, \quad -\sigma^2_i = {\rm cov}(X_i,X_i) = {\rm var}(X_i) = \langle (X_i-\bar X_i)^2\rangle, \quad \ i,j = 1\dots N. -$$ - -Here, $\bar X_i$ denotes mean value and $\langle\dots\rangle$ is ensemble average. The covariance matrix is symmetric about its principal diagonal. The *correlation* between $X_i$ and $X_j$ is defined as the normalized covariance - -$$ -\rho_{ij} = {{\rm cov}(X_i,X_j)\over \sqrt{{\rm var}(X_i)}\sqrt{{\rm var}(X_j)}} = {\sigma_{ij} \over \sigma_i\sigma_j}. -$$ - -Since we are in a Colab environment we will use the Python API, but it is so similar to the JavaScript version that this should pose no problems to anyone not acquainted with it. - - -## Preliminaries - -The cells below execute the necessary formalities for accessing Earth Engine from this Colab notebook. They also import the various Python add-ons needed, including the [earthengine-jupyter](https://github.com/google/earthengine-jupyter) interactive map package, and define a few helper routines for displaying results. - - -```python -import ee - -ee.Authenticate(auth_mode='notebook') -ee.Initialize() -``` - - -```python -# Import other packages used in the tutorial -%matplotlib inline -import geemap -import numpy as np -import random, time -import matplotlib.pyplot as plt -from scipy.stats import norm, chi2 - -from pprint import pprint # for pretty printing -``` - - -```python -# Truncate a 1-D array to dec decimal places -def trunc(values, dec = 3): - return np.trunc(values*10**dec)/(10**dec) -``` - - -```python -# Display an image in a one percent linear stretch -def display_ls(image, map, name, centered = False): - bns = image.bandNames().length().getInfo() - if bns == 3: - image = image.rename('B1', 'B2', 'B3') - pb_99 = ['B1_p99', 'B2_p99', 'B3_p99'] - pb_1 = ['B1_p1', 'B2_p1', 'B3_p1'] - img = ee.Image.rgb(image.select('B1'), image.select('B2'), image.select('B3')) - else: - image = image.rename('B1') - pb_99 = ['B1_p99'] - pb_1 = ['B1_p1'] - img = image.select('B1') - percentiles = image.reduceRegion(ee.Reducer.percentile([1, 99]), maxPixels=1e11) - mx = percentiles.values(pb_99) - if centered: - mn = ee.Array(mx).multiply(-1).toList() - else: - mn = percentiles.values(pb_1) - map.addLayer(img, {'min': mn, 'max': mx}, name) -``` - -## Simple Differences -Let's begin by considering two $N$-band optical/infrared images of the same scene acquired with the same sensor at two different times, between which ground reflectance changes have occurred at -some locations but not everywhere. To see those changes we can "flick back and forth" between gray scale or RGB representations of the images on a computer display. Alternatively we can simply subtract one image from the other and, provided the intensities at the two time points have been calibrated or normalized in some way, examine the difference. Symbolically the difference image is - -$$ -Z = X-Y -$$ - -where - -$$ -X= \begin{pmatrix}X_1\cr X_2\cr\vdots\cr X_N\end{pmatrix},\quad Y= \begin{pmatrix}Y_1\cr Y_2\cr\vdots\cr Y_N\end{pmatrix} -$$ - -are 'typical' pixel vectors (more formally: _random vectors_) representing the first and second image, respectively. - -### Landkreis Olpe - -For example, here is an area of interest (aoi) covering a heavily forested administrative district (Landkreis Olpe) in North Rhine Westphalia, Germany. Due to severe drought in recent years large swaths of shallow-root coniferous trees have died and been cleared away, leaving deciduous trees for the most part untouched. - - -```python -aoi = ee.FeatureCollection('projects/google/imad_tutorial_aoi').geometry() -``` - -We first collect two Sentinel-2 scenes which bracket some of the recent clear cutting (June, 2021 and June, 2022) and print out their timestamps. For demonstration purposes we choose top of atmosphere reflectance since the MAD transformation does not require absolute surface reflectances. - - -```python -def collect(aoi, t1a ,t1b, t2a, t2b): - try: - im1 = ee.Image( ee.ImageCollection("COPERNICUS/S2_HARMONIZED") - .filterBounds(aoi) - .filterDate(ee.Date(t1a), ee.Date(t1b)) - .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo')) - .sort('CLOUDY_PIXEL_PERCENTAGE') - .first() - .clip(aoi) ) - im2 = ee.Image( ee.ImageCollection("COPERNICUS/S2_HARMONIZED") - .filterBounds(aoi) - .filterDate(ee.Date(t2a), ee.Date(t2b)) - .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo')) - .sort('CLOUDY_PIXEL_PERCENTAGE') - .first() - .clip(aoi) ) - timestamp = im1.date().format('E MMM dd HH:mm:ss YYYY') - print(timestamp.getInfo()) - timestamp = im2.date().format('E MMM dd HH:mm:ss YYYY') - print(timestamp.getInfo()) - return (im1, im2) - except Exception as e: - print('Error: %s'%e) - -im1, im2 = collect(aoi, '2021-06-01', '2021-06-30', '2022-06-01', '2022-06-30') -``` - -In order to view the images, we instantiate an interactive map located at the center of the aoi. Since the map requires lat/long coordinates we have to reverse the long/lat order returned from the *coordinates()* method of the *ee.Geometry.centroid()* object: - - -```python -# Interactive map -M1 = geemap.Map() -M1.centerObject(aoi, 11) - -M1 -``` - -Then we display RGB composites of the first 3 (visual) bands of the Sentinel images as well as their difference. We use the previously defined function *display_ls()* to display the images in a one percent linear stretch. - - -```python -visirbands = ['B2','B3','B4','B8','B11','B12'] -visbands = ['B2','B3','B4'] - -diff = im1.subtract(im2).select(visbands) -display_ls(im1.select(visbands), M1, 'Image 1') -display_ls(im2.select(visbands), M1, 'Image 2') -display_ls(diff, M1, 'Difference') -``` - -Small intensity differences (the gray-ish background in the difference image above) indicate no change, large -positive (bright) or negative (dark) colors indicate change. The most obvious changes are in fact due to the clear cutting. Decision thresholds can be set to define significant changes. The thresholds are usually expressed in terms of standard deviations from the mean difference value, which is taken to correspond to no change. The per-band variances of the difference images are simply - -$$ -{\rm var}(Z_i) = {\rm var}(X_i - Y_i) = {\rm var}(X_i) + {\rm var}(Y_i) - 2\cdot{\rm cov}(X_i,Y_i),\quad i=1\dots N, -$$ - -or, if the bands are uncorrelated, about twice as noisy as the individual image bands. Normally $X_i$ and $Y_i$ are positively correlated $({\rm cov}(X_i,Y_i) > 0)$ so the variance of $Z_i$ is somewhat smaller. When the difference signatures in the spectral channels are combined so as to try to characterize the nature of the changes that have taken place, one speaks of _spectral change vector analysis_. - -While the recent clear cuts are fairly easily identified as (some but not all of) the dark areas in the simple difference image above, it is possible to derive a better and more informative change map. This involves taking greater advantage of the statistical properties of the images. - -## The MAD Transformation -Let's make a linear combination of the intensities for all $N$ bands in the first image $X$, thus creating a scalar image - -$$ -U = a_1X_1 + a_2X_2 + \dots + a_NX_N = {\bf a}^\top X. \tag{1} -$$ - -The symbol ${\bf a}^\top$ denotes the transpose of the column vector - -$$ -{\bf a}=\begin{pmatrix}a_1\cr a_2\cr\vdots\cr a_N\end{pmatrix} -$$ - -in other words the row vector ${\bf a}^\top = (a_1,a_2 \dots a_N)$. - - -The expression ${\bf a}^\top X$ is called an _inner vector product_. The vector of coefficients ${\bf a}$ is as yet unspecified. -We'll do the same for the second image $Y$, forming the linear combination - -$$ -V = b_1Y_1 + b_2Y_2 + \dots + b_NY_N = {\bf b}^\top Y, \tag{2} -$$ - -and then look at the scalar difference $U-V$. Change information -is now contained, for the time being, in a single image rather than distributed among all $N$ bands. - -We have of course to choose the coefficient vectors ${\bf a}$ and -${\bf b}$ in some suitable way. In [Nielsen et al. (1998)](https://www2.imm.dtu.dk/pubdb/pubs/1220-full.html) it is suggested that they be determined by making the transformed images $U$ and $V$ _as similar as they can be_ before taking their difference. This sounds at first counter-intuitive: Why make the images similar when we want to see the dissimilarities? The crux is that, in making them match one another by taking suitable linear combinations, we ensure that pixels in the transformed images $U$ and $V$ at which *indeed no change has taken place* will be as similar as possible before they are compared (subtracted). The genuine changes will then, it is hoped, be all the more apparent in the difference image. - -This process of "similar making" can be accomplished in statistics with standard [*Canonical Correlation Analysis* (CCA)](https://en.wikipedia.org/wiki/Canonical_correlation), a technique first described by Harold Hotelling in 1936. Canonical Correlation Analysis of the images represented by Equations (1) and (2) *maximizes the correlation* - -$$ -\rho = {{\rm cov}(U,V)\over \sqrt{{\rm var}(U)}\sqrt{{\rm var}(V)}}, \tag{3} -$$ - -between the random variables $U$ and $V$, which is just what we want. - -One technical point: Arbitrary multiples of $U$ and $V$ would clearly have the same correlation (the multiples will cancel off in numerator and denominator), so a constraint must be chosen. A convenient one is - -$$ -{\rm var}(U)={\rm var}(V)=1. \tag{4} -$$ - -For those more versed in linear algebra the details of CCA are given in [Canty (2019)](https://www.taylorfrancis.com/books/image-analysis-classification-change-detection-remote-sensing-morton-john-canty/10.1201/9780429464348) beginning on page 385. It all boils down to solving two so-called *generalized eigenvalue problems* for determining the transformation vectors ${\bf a}$ and ${\bf b}$ in Equations (1) and (2), respectively. There is not just one solution but rather $N$ solutions. They consist of the $N$ pairs of _eigenvectors_ $({\bf a}^i, {\bf b}^i)$, - -$$ -{\bf a}^i=\begin{pmatrix}a_1\cr a_2\cr\vdots\cr a_N \end{pmatrix}_i,\quad {\bf b}^i=\begin{pmatrix}b_1\cr b_2\cr\vdots\cr b_N \end{pmatrix}_i, \quad i=1\dots N -$$ - -and, correspondingly, the $N$ pairs of *canonical variates* $(U_i, V_i),\ i= 1\dots N$. (Readers familiar with Principal Components Analysis (PCA) will recognize the similarity. When applying PCA to an N-band image, an ordinary (rather than a generalized) eigenvalue problem is solved to maximize variance, resulting in $N$ principal component bands.) - -Unlike the bands of the original images $X$ and $Y$, which are ordered by spectral wavelength, the canonical variates are ordered by similarity or correlation. The difference images - -$$ -M_i = U_i - V_i = ({\bf a}^i)^\top X - ({\bf b}^i)^\top Y,\quad i=1\dots N, \tag{5} -$$ - -contain the change information and are called the *MAD variates*. - - -### Auxiliary functions - -To run the MAD transformation on GEE image pairs we need some helper routines. These are coded below. The first, called *covarw()*, returns the weighted covariance matrix of an $N$-band image as well as a weighted, centered (mean-zero) version of the image. Why it is important to weight the pixels before centering or calculating the covariance matrix will become apparent in Part 2 of the tutorial. - - -```python -def covarw(image, weights = None, scale = 20, maxPixels = 1e10): - '''Return the weighted centered image and its weighted covariance matrix''' - try: - geometry = image.geometry() - bandNames = image.bandNames() - N = bandNames.length() - if weights is None: - weights = image.constant(1) - weightsImage = image.multiply(ee.Image.constant(0)).add(weights) - means = image.addBands(weightsImage) \ - .reduceRegion(ee.Reducer.mean().repeat(N).splitWeights(), - scale = scale, - maxPixels = maxPixels) \ - .toArray() \ - .project([1]) - centered = image.toArray().subtract(means) - B1 = centered.bandNames().get(0) - b1 = weights.bandNames().get(0) - nPixels = ee.Number(centered.reduceRegion(ee.Reducer.count(), - scale = scale, - maxPixels = maxPixels).get(B1)) - sumWeights = ee.Number(weights.reduceRegion(ee.Reducer.sum(), - geometry = geometry, - scale = scale, - maxPixels = maxPixels).get(b1)) - covw = centered.multiply(weights.sqrt()) \ - .toArray() \ - .reduceRegion(ee.Reducer.centeredCovariance(), - geometry = geometry, - scale = scale, - maxPixels = maxPixels) \ - .get('array') - covw = ee.Array(covw).multiply(nPixels).divide(sumWeights) - return (centered.arrayFlatten([bandNames]), covw) - except Exception as e: - print('Error: %s'%e) -``` - -The second routine, called _corr()_, transforms a covariance matrix to the equivalent correlation matrix by dividing by the square roots of the variances, see Equation (3). - - -```python -def corr(cov): - '''Transform covariance matrix to correlation matrix''' - # diagonal matrix of inverse sigmas - sInv = cov.matrixDiagonal().sqrt().matrixToDiag().matrixInverse() - # transform - corr = sInv.matrixMultiply(cov).matrixMultiply(sInv).getInfo() - # truncate - return [list(map(trunc, corr[i])) for i in range(len(corr))] -``` - -For example the code below calculates the (unweighted) covariance matrix of the visual and infrared bands of the image _im1_ and prints out the corresponding correlation matrix. - - -```python -_, cov = covarw(im1.select(visirbands)) -pprint(corr(cov)) -``` - -Note that the spectral bands of the image are generally highly and positively correlated with one another. - -By stacking the images one atop the other we can use the helper functions to display the between image correlations - -$$ -\rm{corr}(X_i,Y_i), \quad i= 1\dots N. -$$ - -They can be found in the diagonal of the upper left $6\times 6$ submatrix of the correlation matrix for the stacked image. - - -```python -im12 = im1.select(visirbands).addBands(im2.select(visirbands)) -_, covar = covarw(im12) -correl = np.array(corr(covar)) -print(np.diag(correl[:6, 6:])) -``` - -It is these between image correlations that the MAD algorithm tries to maximize. - -The third and last auxiliary routine, *geneiv()*, is the core of CCA and the MAD transformation. It solves the *generalized eigenproblem*, - -$$ -CX = \lambda BX -$$ - -for two $N\times N$ matrices $C$ and $B$, returning the $N$ solutions, or eigenvectors, $X_i$ and the $N$ eigenvalues $\lambda_i, \ i=1\dots N$. - - -```python -def geneiv(C,B): - '''Return the eigenvalues and eigenvectors of the generalized eigenproblem - C*X = lambda*B*X''' - try: - C = ee.Array(C) - B = ee.Array(B) - # Li = choldc(B)^-1 - Li = ee.Array(B.matrixCholeskyDecomposition().get('L')).matrixInverse() - # solve symmetric, ordinary eigenproblem Li*C*Li^T*x = lambda*x - Xa = Li.matrixMultiply(C) \ - .matrixMultiply(Li.matrixTranspose()) \ - .eigen() - # eigenvalues as a row vector - lambdas = Xa.slice(1, 0, 1).matrixTranspose() - # eigenvectors as columns - X = Xa.slice(1, 1).matrixTranspose() - # generalized eigenvectors as columns, Li^T*X - eigenvecs = Li.matrixTranspose().matrixMultiply(X) - return (lambdas, eigenvecs) - except Exception as e: - print('Error: %s'%e) -``` - -The next cell codes the MAD transformation itself in the function *mad_run()*, taking as input two multiband images and returning the _canonical variates_ - -$$ -U_i, \ V_i, \quad i=1\dots N, -$$ - -the _MAD variates_ - -$$ -M_i = U_i - V_i, \quad i=1\dots N, -$$ - -as well as the sum of the squares of the standardized MAD variates, - -$$ -Z = \sum_{i=1}^N\left({M_i\over \sigma_{M_i}}\right)^2. -$$ - - -```python -#@title The MAD transformation -def mad_run(image1, image2, scale = 20): - '''The MAD transformation of two multiband images''' - try: - image = image1.addBands(image2) - region = image.geometry() - nBands = image.bandNames().length().divide(2) - centeredImage,covarArray = covarw(image,scale=scale) - bNames = centeredImage.bandNames() - bNames1 = bNames.slice(0,nBands) - bNames2 = bNames.slice(nBands) - centeredImage1 = centeredImage.select(bNames1) - centeredImage2 = centeredImage.select(bNames2) - s11 = covarArray.slice(0,0,nBands).slice(1,0,nBands) - s22 = covarArray.slice(0,nBands).slice(1,nBands) - s12 = covarArray.slice(0,0,nBands).slice(1,nBands) - s21 = covarArray.slice(0,nBands).slice(1,0,nBands) - c1 = s12.matrixMultiply(s22.matrixInverse()).matrixMultiply(s21) - b1 = s11 - c2 = s21.matrixMultiply(s11.matrixInverse()).matrixMultiply(s12) - b2 = s22 - # solution of generalized eigenproblems - lambdas, A = geneiv(c1,b1) - _, B = geneiv(c2,b2) - rhos = lambdas.sqrt().project(ee.List([1])) - # MAD variances - sigma2s = rhos.subtract(1).multiply(-2).toList() - sigma2s = ee.Image.constant(sigma2s) - # ensure sum of correlations between X and U is positive - tmp = s11.matrixDiagonal().sqrt() - ones = tmp.multiply(0).add(1) - tmp = ones.divide(tmp).matrixToDiag() - s = tmp.matrixMultiply(s11).matrixMultiply(A).reduce(ee.Reducer.sum(),[0]).transpose() - A = A.matrixMultiply(s.divide(s.abs()).matrixToDiag()) - # ensure positive correlation between U and V - tmp = A.transpose().matrixMultiply(s12).matrixMultiply(B).matrixDiagonal() - tmp = tmp.divide(tmp.abs()).matrixToDiag() - B = B.matrixMultiply(tmp) - # canonical and MAD variates - centeredImage1Array = centeredImage1.toArray().toArray(1) - centeredImage2Array = centeredImage2.toArray().toArray(1) - U = ee.Image(A.transpose()).matrixMultiply(centeredImage1Array) \ - .arrayProject([0]) \ - .arrayFlatten([bNames2]) - V = ee.Image(B.transpose()).matrixMultiply(centeredImage2Array) \ - .arrayProject([0]) \ - .arrayFlatten([bNames2]) - MAD = U.subtract(V) - # chi square image - Z = MAD.pow(2) \ - .divide(sigma2s) \ - .reduce(ee.Reducer.sum()) \ - .clip(region) - return (U, V, MAD, Z) - except Exception as e: - print('Error: %s'%e) - -``` - -After setting up the MAD transformation on the six visual and infrared bands of the Landkreis Olpe images, - - -```python -U, V, MAD, Z = mad_run(im1.select(visirbands), im2.select(visirbands)) -``` - -we display an RGB composite of the first 3 MAD variates, together with the original images and their simple difference: - - -```python -M2 = geemap.Map() -M2.centerObject(aoi, 11) -display_ls(im1.select(visbands), M2, 'Image 1') -display_ls(im2.select(visbands), M2, 'Image 2') -display_ls(diff, M2, 'Difference') -display_ls(MAD.select(0, 1, 2), M2, 'MAD Image', True) - -M2 -``` - -The richness in colors in the MAD image is a consequence of the decoupling (orthogonalization) of the MAD components, as will be explained below. - -Pretty colors notwithstanding, when compared with the simple difference image the result is _mittelprächtig_ (German slang for "not so hot"), and certainly no easier to interpret. However we're not finished. - -The canonical variates have very nice properties indeed. They are _all mutually uncorrelated_ except for the pairs $(U_i,V_i)$, and these are ordered by decreasing positive correlation: - -$$ -\rho_i = {\rm cov}(U_i, V_i),\quad i=1\dots N, \tag{6} -$$ - -with $\rho_1 \ge \rho_2 \ge\dots \ge\rho_N$. - -Stacking the images $U$ and $V$, these facts can be read from the $12\times 12$ correlation matrix: - - -```python -_, covar = covarw(U.addBands(V)) -correl = np.array(corr(covar)) -pprint(correl) -print('rho =', np.diag(correl[:6,6:])) -``` - -The MAD variates themselves are consequently also mutually uncorrelated, their covariances being given by - -$$ -{\rm cov}(M_i,M_j) = {\rm cov}(U_i-V_i,U_j-V_j)= 0,\quad i\ne j=1\dots N, \tag{7} -$$ - -and their variances by - -$$ -\sigma_{M_i}^2={\rm var}(U_i-V_i)= {\rm var}(U_i)+{\rm var}(V_i) -2{\rm cov}(U_i,V_i) = 2(1-\rho_i),\quad i=1\dots N, \tag{8} -$$ - -where the last equality follows from the constraint Equation (4). Let's check. - - -```python -# display MAD covariance matrix -_, covar = covarw(MAD) -covar = covar.getInfo() -[list(map(trunc,covar[i])) for i in range(len(covar))] -``` - - - -The first MAD variate has minimum variance in its pixel -intensities. The second MAD variate has minimum variance subject to the condition that its pixel intensities are statistically -uncorrelated with those in the first variate; the third has -minimum spread subject to being uncorrelated with the first two, -and so on. Depending on the type of change present, any of the -components may exhibit significant change information, although it tends to be concentrated in the first few MAD variates. One of the nicest aspects of the MAD transformation is that it can sort different categories of change into different, uncorrelated image bands. - -However, the result, as we saw, needs some improvement. This is the subject of Part 2. diff --git a/tutorials/imad-tutorial-pt2/index.ipynb b/tutorials/imad-tutorial-pt2/index.ipynb new file mode 100644 index 000000000..5130079ca --- /dev/null +++ b/tutorials/imad-tutorial-pt2/index.ipynb @@ -0,0 +1,1114 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "9064f134", + "metadata": { + "id": "9064f134" + }, + "outputs": [], + "source": [ + "#@title Copyright 2023 The Earth Engine Community Authors { display-mode: \"form\" }\n", + "#\n", + "# Licensed under the Apache License, Version 2.0 (the \"License\");\n", + "# you may not use this file except in compliance with the License.\n", + "# You may obtain a copy of the License at\n", + "#\n", + "# https://www.apache.org/licenses/LICENSE-2.0\n", + "#\n", + "# Unless required by applicable law or agreed to in writing, software\n", + "# distributed under the License is distributed on an \"AS IS\" BASIS,\n", + "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n", + "# See the License for the specific language governing permissions and\n", + "# limitations under the License." + ] + }, + { + "cell_type": "markdown", + "id": "31afb58c", + "metadata": { + "id": "31afb58c" + }, + "source": [ + "# Change Detection on Google Earth Engine\n", + "# Part 2 The iMAD Algorithm\n", + "Authors: mortcanty" + ] + }, + { + "cell_type": "markdown", + "id": "4040ed83", + "metadata": { + "id": "4040ed83" + }, + "source": [ + "## Context\n", + "\n", + "In Part 1 of this tutorial, a statistical approach to detecting changes in pairs of\n", + "multispectral remote sensing images, called Multivariate Alteration Detection (MAD), was\n", + "explained. The MAD change images were obtained by first maximizing the correlations\n", + "between the original image bands and then subtracting one from the other. The resulting\n", + "difference bands, called MAD variates, contained the change information. They were shown\n", + "to be ordered by increasing variance and to be mutually uncorrelated. However, they were\n", + "also seen to be not easily interpretable.\n", + "\n", + "Now, in Part 2, we introduce an iteration scheme that performs the MAD transformation\n", + "exclusively on those pixels that mark areas in the images which have not physically\n", + "changed. This establishes a well-defined background of invariant pixels against which\n", + "to discriminate changes.\n" + ] + }, + { + "cell_type": "markdown", + "id": "912fe8bf", + "metadata": { + "id": "912fe8bf" + }, + "source": [ + "## Preliminaries" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0ca8fe9a", + "metadata": { + "id": "0ca8fe9a" + }, + "outputs": [], + "source": [ + "# Enter your own export to assets path name here -----------\n", + "EXPORT_PATH = 'projects/YOUR_GEE_PROJECT_NAME/assets/imad/'\n", + "# ------------------------------------------------" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "511ab83e", + "metadata": { + "id": "511ab83e", + "tags": [] + }, + "outputs": [], + "source": [ + "import ee\n", + "\n", + "ee.Authenticate(auth_mode='notebook')\n", + "ee.Initialize()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f35057db", + "metadata": { + "id": "f35057db" + }, + "outputs": [], + "source": [ + "# Import other packages used in the tutorial\n", + "%matplotlib inline\n", + "import geemap\n", + "import numpy as np\n", + "import random, time\n", + "import matplotlib.pyplot as plt\n", + "from scipy.stats import norm, chi2\n", + "\n", + "from pprint import pprint # for pretty printing" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47df0458", + "metadata": { + "id": "47df0458" + }, + "outputs": [], + "source": [ + "#@title Routines from Part 1\n", + "\n", + "def trunc(values, dec = 3):\n", + " '''Truncate a 1-D array to dec decimal places.'''\n", + " return np.trunc(values*10**dec)/(10**dec)\n", + "\n", + "# Display an image in a one percent linear stretch.\n", + "def display_ls(image, map, name, centered = False):\n", + " bns = image.bandNames().length().getInfo()\n", + " if bns == 3:\n", + " image = image.rename('B1', 'B2', 'B3')\n", + " pb_99 = ['B1_p99', 'B2_p99', 'B3_p99']\n", + " pb_1 = ['B1_p1', 'B2_p1', 'B3_p1']\n", + " img = ee.Image.rgb(image.select('B1'), image.select('B2'), image.select('B3'))\n", + " else:\n", + " image = image.rename('B1')\n", + " pb_99 = ['B1_p99']\n", + " pb_1 = ['B1_p1']\n", + " img = image.select('B1')\n", + " percentiles = image.reduceRegion(ee.Reducer.percentile([1, 99]), maxPixels=1e11)\n", + " mx = percentiles.values(pb_99)\n", + " if centered:\n", + " mn = ee.Array(mx).multiply(-1).toList()\n", + " else:\n", + " mn = percentiles.values(pb_1)\n", + " map.addLayer(img, {'min': mn, 'max': mx}, name)\n", + "\n", + "def collect(aoi, t1a ,t1b, t2a, t2b):\n", + " try:\n", + " im1 = ee.Image( ee.ImageCollection(\"COPERNICUS/S2_HARMONIZED\")\n", + " .filterBounds(aoi)\n", + " .filterDate(ee.Date(t1a), ee.Date(t1b))\n", + " .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo'))\n", + " .sort('CLOUDY_PIXEL_PERCENTAGE')\n", + " .first()\n", + " .clip(aoi) )\n", + " im2 = ee.Image( ee.ImageCollection(\"COPERNICUS/S2_HARMONIZED\")\n", + " .filterBounds(aoi)\n", + " .filterDate(ee.Date(t2a), ee.Date(t2b))\n", + " .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo'))\n", + " .sort('CLOUDY_PIXEL_PERCENTAGE')\n", + " .first()\n", + " .clip(aoi) )\n", + " timestamp = im1.date().format('E MMM dd HH:mm:ss YYYY')\n", + " print(timestamp.getInfo())\n", + " timestamp = im2.date().format('E MMM dd HH:mm:ss YYYY')\n", + " print(timestamp.getInfo())\n", + " return (im1, im2)\n", + " except Exception as e:\n", + " print('Error: %s'%e)\n", + "\n", + "def covarw(image, weights=None, scale=20, maxPixels=1e10):\n", + " '''Return the centered image and its weighted covariance matrix.'''\n", + " try:\n", + " geometry = image.geometry()\n", + " bandNames = image.bandNames()\n", + " N = bandNames.length()\n", + " if weights is None:\n", + " weights = image.constant(1)\n", + " weightsImage = image.multiply(ee.Image.constant(0)).add(weights)\n", + " means = image.addBands(weightsImage) \\\n", + " .reduceRegion(ee.Reducer.mean().repeat(N).splitWeights(),\n", + " scale = scale,\n", + " maxPixels = maxPixels) \\\n", + " .toArray() \\\n", + " .project([1])\n", + " centered = image.toArray().subtract(means)\n", + " B1 = centered.bandNames().get(0)\n", + " b1 = weights.bandNames().get(0)\n", + " nPixels = ee.Number(centered.reduceRegion(ee.Reducer.count(),\n", + " scale=scale,\n", + " maxPixels=maxPixels).get(B1))\n", + " sumWeights = ee.Number(weights.reduceRegion(ee.Reducer.sum(),\n", + " geometry=geometry,\n", + " scale=scale,\n", + " maxPixels=maxPixels).get(b1))\n", + " covw = centered.multiply(weights.sqrt()) \\\n", + " .toArray() \\\n", + " .reduceRegion(ee.Reducer.centeredCovariance(),\n", + " geometry=geometry,\n", + " scale=scale,\n", + " maxPixels=maxPixels) \\\n", + " .get('array')\n", + " covw = ee.Array(covw).multiply(nPixels).divide(sumWeights)\n", + " return (centered.arrayFlatten([bandNames]), covw)\n", + " except Exception as e:\n", + " print('Error: %s'%e)\n", + "\n", + "def corr(cov):\n", + " '''Transfrom covariance matrix to correlation matrix.'''\n", + " # Diagonal matrix of inverse sigmas.\n", + " sInv = cov.matrixDiagonal().sqrt().matrixToDiag().matrixInverse()\n", + " # Transform.\n", + " corr = sInv.matrixMultiply(cov).matrixMultiply(sInv).getInfo()\n", + " # Truncate.\n", + " return [list(map(trunc, corr[i])) for i in range(len(corr))]\n", + "\n", + "def geneiv(C,B):\n", + " '''Return the eignvalues and eigenvectors of the generalized eigenproblem\n", + " C*X = lambda*B*X'''\n", + " try:\n", + " C = ee.Array(C)\n", + " B = ee.Array(B)\n", + " # Li = choldc(B)^-1\n", + " Li = ee.Array(B.matrixCholeskyDecomposition().get('L')).matrixInverse()\n", + " # Solve symmetric, ordinary eigenproblem Li*C*Li^T*x = lambda*x\n", + " Xa = Li.matrixMultiply(C) \\\n", + " .matrixMultiply(Li.matrixTranspose()) \\\n", + " .eigen()\n", + " # Eigenvalues as a row vector.\n", + " lambdas = Xa.slice(1, 0, 1).matrixTranspose()\n", + " # Eigenvectors as columns.\n", + " X = Xa.slice(1, 1).matrixTranspose()\n", + " # Generalized eigenvectors as columns, Li^T*X\n", + " eigenvecs = Li.matrixTranspose().matrixMultiply(X)\n", + " return (lambdas, eigenvecs)\n", + " except Exception as e:\n", + " print('Error: %s'%e)\n", + "\n", + "def mad_run(image1, image2, scale=20):\n", + " '''The MAD transformation of two multiband images.'''\n", + " try:\n", + " image = image1.addBands(image2)\n", + " nBands = image.bandNames().length().divide(2)\n", + " centeredImage,covarArray = covarw(image,scale=scale)\n", + " bNames = centeredImage.bandNames()\n", + " bNames1 = bNames.slice(0,nBands)\n", + " bNames2 = bNames.slice(nBands)\n", + " centeredImage1 = centeredImage.select(bNames1)\n", + " centeredImage2 = centeredImage.select(bNames2)\n", + " s11 = covarArray.slice(0, 0, nBands).slice(1, 0, nBands)\n", + " s22 = covarArray.slice(0, nBands).slice(1, nBands)\n", + " s12 = covarArray.slice(0, 0, nBands).slice(1, nBands)\n", + " s21 = covarArray.slice(0, nBands).slice(1, 0, nBands)\n", + " c1 = s12.matrixMultiply(s22.matrixInverse()).matrixMultiply(s21)\n", + " b1 = s11\n", + " c2 = s21.matrixMultiply(s11.matrixInverse()).matrixMultiply(s12)\n", + " b2 = s22\n", + " # Solution of generalized eigenproblems.\n", + " lambdas, A = geneiv(c1, b1)\n", + " _, B = geneiv(c2, b2)\n", + " rhos = lambdas.sqrt().project(ee.List([1]))\n", + " # MAD variances.\n", + " sigma2s = rhos.subtract(1).multiply(-2).toList()\n", + " sigma2s = ee.Image.constant(sigma2s)\n", + " # Ensure sum of positive correlations between X and U is positive.\n", + " tmp = s11.matrixDiagonal().sqrt()\n", + " ones = tmp.multiply(0).add(1)\n", + " tmp = ones.divide(tmp).matrixToDiag()\n", + " s = tmp.matrixMultiply(s11).matrixMultiply(A).reduce(ee.Reducer.sum(),[0]).transpose()\n", + " A = A.matrixMultiply(s.divide(s.abs()).matrixToDiag())\n", + " # Ensure positive correlation.\n", + " tmp = A.transpose().matrixMultiply(s12).matrixMultiply(B).matrixDiagonal()\n", + " tmp = tmp.divide(tmp.abs()).matrixToDiag()\n", + " B = B.matrixMultiply(tmp)\n", + " # Canonical and MAD variates as images.\n", + " centeredImage1Array = centeredImage1.toArray().toArray(1)\n", + " centeredImage2Array = centeredImage2.toArray().toArray(1)\n", + " U = ee.Image(A.transpose()).matrixMultiply(centeredImage1Array) \\\n", + " .arrayProject([0]) \\\n", + " .arrayFlatten([bNames2])\n", + " V = ee.Image(B.transpose()).matrixMultiply(centeredImage2Array) \\\n", + " .arrayProject([0]) \\\n", + " .arrayFlatten([bNames2])\n", + " MAD = U.subtract(V)\n", + " # Chi-square image.\n", + " Z = MAD.pow(2) \\\n", + " .divide(sigma2s) \\\n", + " .reduce(ee.Reducer.sum())\n", + " return (U, V, MAD, Z)\n", + " except Exception as e:\n", + " print('Error: %s'%e)" + ] + }, + { + "cell_type": "markdown", + "id": "e1bfac16", + "metadata": { + "id": "e1bfac16" + }, + "source": [ + "## Iterative re-weighting\n", + "\n", + "Let's consider two images of the same scene acquired at different times under similar\n", + "conditions, similar to the Sentinel-2 images from Part 1, but for which no ground\n", + "reflectance changes have occurred. In this case, the only differences between the images\n", + "will be due to random effects such as instrument noise and atmospheric fluctuations.\n", + "Therefore, we can expect that the histogram of any difference component we generate will\n", + "be nearly Gaussian. Specifically, the MAD variates, which we have seen to be uncorrelated,\n", + "should follow a multivariate, zero-mean normal distribution with a diagonal covariance\n", + "matrix.\n", + "\n", + "$$\n", + "\\Sigma_M = \\pmatrix{\\sigma^2_{M_1} \u00260 \u0026\\cdots \u00260 \\cr\n", + " 0 \u0026 \\sigma^2_{M_2} \u0026\\cdots \u00260 \\cr\n", + " \\vdots \u0026\\vdots \u0026\\cdots \u00260 \\cr\n", + " 0 \u0026 0 \u0026\\cdots \u0026\\sigma^2_{M_N}}.\n", + "$$\n", + "\n", + "Change observations would deviate more or less strongly from such a\n", + "distribution. We might therefore expect an improvement in the sensitivity\n", + "of the MAD transformation *if we can establish a better background of no\n", + "change against which to detect change.* This can be done in an iteration\n", + "scheme ([Nielsen 2007](https://www2.imm.dtu.dk/pubdb/pubs/4695-full.html))\n", + "in which, when calculating the statistics for each successive iteration\n", + "of the MAD transformation, observations are weighted in some appropriate\n", + "fashion.\n", + "\n", + "Recall that the variable $Z$ represents the sum of the squares of the standardized MAD variates,\n", + "\n", + "$$\n", + "Z = \\sum_{i=1}^N\\left({M_i\\over \\sigma_{M_i}}\\right)^2,\n", + "$$\n", + "\n", + "where $\\sigma^2_{M_i}$ is given by Equation (8) in the first Tutorial,\n", + "\n", + "$$\n", + "\\sigma_{M_i}^2={\\rm var}(U_i-V_i) = 2(1-\\rho_i),\\quad i=1\\dots N,\n", + "$$\n", + "\n", + "and $\\rho_i = {\\rm cov}(U_i,V_i)$. Then, since the no-change observations are expected\n", + "to be normally distributed and uncorrelated, basic statistical theory tells us\n", + "that the values of $Z$ corresponding to no-change observations should be *chi-square\n", + "distributed* with $N$ degrees of freedom. Let's check to what extent this is true for\n", + "the MAD variates that we have determined so far for the Landkreis Olpe scene." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4be3d71f", + "metadata": { + "id": "4be3d71f" + }, + "outputs": [], + "source": [ + "# Landkreis Olpe.\n", + "aoi = ee.FeatureCollection('projects/google/imad_tutorial_aoi').geometry()\n", + "\n", + "visirbands = ['B2', 'B3', 'B4', 'B8', 'B11', 'B12']\n", + "visbands = ['B2', 'B3', 'B4']\n", + "rededgebands = ['B5', 'B6', 'B7', 'B8A']\n", + "\n", + "# Collect the two Sentinel-2 images.\n", + "im1, im2 = collect(aoi, '2021-06-01', '2021-06-30', '2022-06-01', '2022-06-30')\n", + "\n", + "# Re-run MAD.\n", + "U, V, MAD, Z = mad_run(im1.select(visirbands), im2.select(visirbands), scale=20)\n", + "\n", + "# Plot histogram of Z.\n", + "hist = Z.reduceRegion(ee.Reducer.fixedHistogram(0, 50, 500), aoi, scale=20).get('sum').getInfo()\n", + "a = np.array(hist)\n", + "x = a[:, 0] # array of bucket edge positions\n", + "y = a[:, 1]/np.sum(a[:, 1]) # normalized array of bucket contents\n", + "plt.plot(x, y, '.', label='data')\n", + "# The chi-square distribution with 6 degrees of freedom.\n", + "plt.plot(x, chi2.pdf(x, 6)/10, '-r', label='Chi-square')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "9adf096b", + "metadata": { + "id": "9adf096b" + }, + "source": [ + "Clearly not the case at all. Which is to be expected, since there are many change pixels in\n", + "the scene and we have made no attempt to discriminate them.\n", + "\n", + "In fact, it is easy to show that $Z$ is a *likelihood ratio test statistic* for change, see\n", + "the discussion of statistical hypothesis testing in the\n", + "[SAR Tutorial](https://developers.google.com/earth-engine/tutorials/community/detecting-changes-in-sentinel-1-imagery-pt-2)\n", + "and the discussion on pp.390-391 in\n", + "[Canty (2019)](https://www.taylorfrancis.com/books/image-analysis-classification-change-detection-remote-sensing-morton-john-canty/10.1201/9780429464348).\n", + "Under the hypothesis that no change has occurred, the test statistic $Z$ will follow, as we\n", + "said, a chi-square distribution. The so-called $p$-*value* is a measure of the extent to which\n", + "this is true. For an observation $z$ of the test statistic, the $p$-value is the probability\n", + "that any sample drawn from the chi-square distribution could be as large as $z$ or larger.\n", + "This is given by\n", + "\n", + "$$\n", + "p(z) = 1-P_{\\chi^2;N}(z),\\quad 0 \u003c p(z) \u003c 1,\n", + "$$\n", + "\n", + "where $P_{\\chi^2;N}(z)$ is the cumulative chi-square probability distribution, i.e., the area\n", + "under the chi-square distribution up to the value $z$, and $p(z)$ is its complement.\n", + "All $p$-values are equally likely if no change has occurred at that pixel\n", + "location$^\\star$, but **change will always be associated with small $p$-values**.\n", + "Therefore in order to eliminate the change observations from the MAD transformation,\n", + "the $p$-value itself can be used to weight each pixel before re-sampling the images to\n", + "determine the statistics for the next iteration. (This was the motivation for coding\n", + "a *weighted* covariance matrix routine in Part 1 earlier). The influence of the\n", + "change observations on the MAD transformation is thereby gradually reduced. Iteration\n", + "continues until some stopping criterion is met, such as lack of significant change in\n", + "the canonical correlations $\\rho_i$. The whole procedure constitutes the *iMAD algorithm*.\n", + "It is implemented in the GEE Python API in the following cell:\n", + "\n", + "$\\star$ Thus the $p$-value is *not* a no-change probability, a common misapprehension!\n", + "See again the [SAR Tutorial](https://developers.google.com/earth-engine/tutorials/community/detecting-changes-in-sentinel-1-imagery-pt-2)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cbcdaad1", + "metadata": { + "cellView": "form", + "id": "cbcdaad1" + }, + "outputs": [], + "source": [ + "#@title The iMAD code\n", + "def chi2cdf(Z,df):\n", + " '''Chi-square cumulative distribution function with df degrees of freedom.'''\n", + " return ee.Image(Z.divide(2)).gammainc(ee.Number(df).divide(2))\n", + "\n", + "def imad(current,prev):\n", + " '''Iterator function for iMAD.'''\n", + " done = ee.Number(ee.Dictionary(prev).get('done'))\n", + " return ee.Algorithms.If(done, prev, imad1(current, prev))\n", + "\n", + "def imad1(current,prev):\n", + " '''Iteratively re-weighted MAD.'''\n", + " image = ee.Image(ee.Dictionary(prev).get('image'))\n", + " Z = ee.Image(ee.Dictionary(prev).get('Z'))\n", + " allrhos = ee.List(ee.Dictionary(prev).get('allrhos'))\n", + " nBands = image.bandNames().length().divide(2)\n", + " weights = chi2cdf(Z,nBands).subtract(1).multiply(-1)\n", + " scale = ee.Dictionary(prev).getNumber('scale')\n", + " niter = ee.Dictionary(prev).getNumber('niter')\n", + " # Weighted stacked image and weighted covariance matrix.\n", + " centeredImage, covarArray = covarw(image, weights, scale)\n", + " bNames = centeredImage.bandNames()\n", + " bNames1 = bNames.slice(0, nBands)\n", + " bNames2 = bNames.slice(nBands)\n", + " centeredImage1 = centeredImage.select(bNames1)\n", + " centeredImage2 = centeredImage.select(bNames2)\n", + " s11 = covarArray.slice(0, 0, nBands).slice(1, 0, nBands)\n", + " s22 = covarArray.slice(0, nBands).slice(1, nBands)\n", + " s12 = covarArray.slice(0, 0, nBands).slice(1, nBands)\n", + " s21 = covarArray.slice(0, nBands).slice(1, 0, nBands)\n", + " c1 = s12.matrixMultiply(s22.matrixInverse()).matrixMultiply(s21)\n", + " b1 = s11\n", + " c2 = s21.matrixMultiply(s11.matrixInverse()).matrixMultiply(s12)\n", + " b2 = s22\n", + " # Solution of generalized eigenproblems.\n", + " lambdas, A = geneiv(c1, b1)\n", + " _, B = geneiv(c2, b2)\n", + " rhos = lambdas.sqrt().project(ee.List([1]))\n", + " # Test for convergence.\n", + " lastrhos = ee.Array(allrhos.get(-1))\n", + " done = rhos.subtract(lastrhos) \\\n", + " .abs() \\\n", + " .reduce(ee.Reducer.max(), ee.List([0])) \\\n", + " .lt(ee.Number(0.0001)) \\\n", + " .toList() \\\n", + " .get(0)\n", + " allrhos = allrhos.cat([rhos.toList()])\n", + " # MAD variances.\n", + " sigma2s = rhos.subtract(1).multiply(-2).toList()\n", + " sigma2s = ee.Image.constant(sigma2s)\n", + " # Ensure sum of positive correlations between X and U is positive.\n", + " tmp = s11.matrixDiagonal().sqrt()\n", + " ones = tmp.multiply(0).add(1)\n", + " tmp = ones.divide(tmp).matrixToDiag()\n", + " s = tmp.matrixMultiply(s11).matrixMultiply(A).reduce(ee.Reducer.sum(), [0]).transpose()\n", + " A = A.matrixMultiply(s.divide(s.abs()).matrixToDiag())\n", + " # Ensure positive correlation.\n", + " tmp = A.transpose().matrixMultiply(s12).matrixMultiply(B).matrixDiagonal()\n", + " tmp = tmp.divide(tmp.abs()).matrixToDiag()\n", + " B = B.matrixMultiply(tmp)\n", + " # Canonical and MAD variates.\n", + " centeredImage1Array = centeredImage1.toArray().toArray(1)\n", + " centeredImage2Array = centeredImage2.toArray().toArray(1)\n", + " U = ee.Image(A.transpose()).matrixMultiply(centeredImage1Array) \\\n", + " .arrayProject([0]) \\\n", + " .arrayFlatten([bNames1])\n", + " V = ee.Image(B.transpose()).matrixMultiply(centeredImage2Array) \\\n", + " .arrayProject([0]) \\\n", + " .arrayFlatten([bNames2])\n", + " iMAD = U.subtract(V)\n", + " # Chi-square image.\n", + " Z = iMAD.pow(2) \\\n", + " .divide(sigma2s) \\\n", + " .reduce(ee.Reducer.sum())\n", + " return ee.Dictionary({'done': done, 'scale': scale, 'niter': niter.add(1),\n", + " 'image': image, 'allrhos': allrhos, 'Z': Z, 'iMAD': iMAD})" + ] + }, + { + "cell_type": "markdown", + "id": "50da8871", + "metadata": { + "id": "50da8871" + }, + "source": [ + "The following code cell is a routine to run the iMAD algorithm as an export task,\n", + "avoiding memory and time limitations in the active runtime. The asset will be exported\n", + "to the location specified by the `EXPORT_PATH` variable defined earlier. It requires\n", + "about 130 MB of space and can take 15 to 20 minutes to complete." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a49bba67", + "metadata": { + "id": "a49bba67" + }, + "outputs": [], + "source": [ + "#@title Run iMAD algorithm as export task\n", + "def run_imad(aoi, image1, image2, assetFN, scale=20, maxiter=100):\n", + " try:\n", + " N = image1.bandNames().length().getInfo()\n", + " imadnames = ['iMAD'+str(i+1) for i in range(N)]\n", + " imadnames.append('Z')\n", + " # Maximum iterations.\n", + " inputlist = ee.List.sequence(1, maxiter)\n", + " first = ee.Dictionary({'done':0,\n", + " 'scale': scale,\n", + " 'niter': ee.Number(0),\n", + " 'image': image1.addBands(image2),\n", + " 'allrhos': [ee.List.sequence(1, N)],\n", + " 'Z': ee.Image.constant(0),\n", + " 'iMAD': ee.Image.constant(0)})\n", + " # Iteration.\n", + " result = ee.Dictionary(inputlist.iterate(imad, first))\n", + " # Retrieve results.\n", + " iMAD = ee.Image(result.get('iMAD')).clip(aoi)\n", + " rhos = ee.String.encodeJSON(ee.List(result.get('allrhos')).get(-1))\n", + " Z = ee.Image(result.get('Z'))\n", + " niter = result.getNumber('niter')\n", + " # Export iMAD and Z as a singe image, including rhos and number of iterations in properties.\n", + " iMAD_export = ee.Image.cat(iMAD, Z).rename(imadnames).set('rhos', rhos, 'niter', niter)\n", + " assetId = EXPORT_PATH + assetFN\n", + " assexport = ee.batch.Export.image.toAsset(iMAD_export,\n", + " description='assetExportTask',\n", + " assetId=assetId, scale=scale, maxPixels=1e10)\n", + " assexport.start()\n", + " print('Exporting iMAD to %s\\n task id: %s'%(assetId, str(assexport.id)))\n", + " except Exception as e:\n", + " print('Error: %s'%e)" + ] + }, + { + "cell_type": "markdown", + "id": "0c50b992", + "metadata": { + "id": "0c50b992" + }, + "source": [ + "and here we run it on our two images:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d3616dd", + "metadata": { + "id": "9d3616dd" + }, + "outputs": [], + "source": [ + "run_imad(aoi, im1.select(visirbands), im2.select(visirbands), 'LandkreisOlpe')" + ] + }, + { + "cell_type": "markdown", + "id": "21061d67", + "metadata": { + "id": "21061d67" + }, + "source": [ + "After the export finishes, the number of iterations and the final canonical correlations\n", + "can be read from properties of the exported image:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a67a360d", + "metadata": { + "id": "a67a360d" + }, + "outputs": [], + "source": [ + "im_imad = ee.Image(EXPORT_PATH + 'LandkreisOlpe').select(0, 1, 2, 3, 4, 5)\n", + "im_z = ee.Image(EXPORT_PATH + 'LandkreisOlpe').select(6).rename('Z')\n", + "niter = im_imad.get('niter').getInfo()\n", + "rhos = ee.List(im_imad.get('rhos')).getInfo()\n", + "print('iteratons: %i'%niter)\n", + "print('canonical correlations: %s'%rhos)" + ] + }, + { + "cell_type": "markdown", + "id": "9b71f923", + "metadata": { + "id": "9b71f923" + }, + "source": [ + "We got convergence after 28 iterations, and the correlations are very close to one\n", + "for the first canonical variates. It might now be interesting to check if the test\n", + "statistic $Z$ has the expected chi-square distribution when evaluated for the no\n", + "change pixels. To to eliminate the changes at the 10% significance level we set a\n", + "lower threshold of $\\alpha = 0.1$ on the $p$-values (recall: small p-values signify change)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "754fd67b", + "metadata": { + "id": "754fd67b" + }, + "outputs": [], + "source": [ + "scale = 20\n", + "# p-values image.\n", + "pval = chi2cdf(im_z, 6).subtract(1).multiply(-1).rename('pval')\n", + "# No-change mask (use p-values greater than 0.1).\n", + "noChangeMask = pval.gt(0.1)\n", + "hist = im_z.updateMask(noChangeMask).reduceRegion(ee.Reducer \\\n", + " .fixedHistogram(0, 50, 500), aoi, scale=scale, maxPixels=1e11) \\\n", + " .get('Z').getInfo()\n", + "a = np.array(hist)\n", + "x = a[:, 0] # array of bucket edge positions\n", + "y = a[:, 1]/np.sum(a[:, 1]) # normalized array of bucket contents\n", + "plt.plot(x, y, '.', label = 'data')\n", + "plt.plot(x, chi2.pdf(x, 6)/10, '-r', label='chi2(6)')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "8e13d7bd", + "metadata": { + "id": "8e13d7bd" + }, + "source": [ + "Agreement is not perfect, but the plot is certainly closer to the ideal chi-square distribution\n", + "after iteration than after the single MAD transformation. So let's display the *iMAD* image:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0e5a6d31", + "metadata": { + "id": "0e5a6d31" + }, + "outputs": [], + "source": [ + "M1 = geemap.Map()\n", + "M1.centerObject(aoi, 11)\n", + "display_ls(im1.select(visbands), M1, 'im1')\n", + "display_ls(im2.select(visbands), M1, 'im2')\n", + "display_ls(im_imad.select('iMAD1', 'iMAD2', 'iMAD3'), M1, 'iMAD123', True)\n", + "\n", + "M1" + ] + }, + { + "cell_type": "markdown", + "id": "aebc78ac", + "metadata": { + "id": "aebc78ac" + }, + "source": [ + "Gray pixels point to no change, while the wide range of color in the iMAD variates\n", + "indicates a good discrimination of the types of change occuring.\n", + "\n", + "**Aside:** We are of course primarily interested in extracting the changes in the iMAD\n", + "image, especially those which mark clear cutting, and we'll come back to them in a moment.\n", + "However, now that what we think the no change pixels have been isolated, we could\n", + "also perform a regression analysis on them to determine how well the radiometric parameters\n", + "of the two Sentinel-2 acquisitions compare. If surface rather than top of atmosphere (TOA)\n", + "reflectance images had been used for the example, we would expect a good match, i.e., a\n", + "slope close to one and an intercept close to zero at all spectral wavelengths. In general,\n", + "for uncalibrated images, this will not be the case. In that event, the regression\n", + "coefficients can be applied to normalize one image (the target, say) to the other\n", + "(the reference). This might be desirable for tracing features such as NDVI indices over\n", + "a time series of acquisitions when the images have not been reduced to surface reflectances,\n", + "see e.g., [Gan et al. (2021)](https://ieeexplore.ieee.org/document/9392311), or indeed for\n", + "*harmonizing* the data from two different sensors of the same family such as Landsat 7 with\n", + "Landsat 8. These topics will be the subject of Part 3.\n", + "\n", + "But now let's look in more detail at the changes in the Landkreis Olpe scene." + ] + }, + { + "cell_type": "markdown", + "id": "99616ab1", + "metadata": { + "id": "99616ab1" + }, + "source": [ + "### Clustering\n", + "\n", + "To better interpret the change image, we can attempt an unsupervised classification.\n", + "We'll see if we can get away with an ordinary K-means clusterer and a simple Euclidean distance\n", + "measure for the complete iMAD image. We choose the number of clusters as 4 and leave all 12(!)\n", + "other input parameters of the *wekaKmeans()* clusterer at their default values. We will also\n", + "first standardize the iMAD image by dividing by the square root of the variances of the no-change\n", + "pixels, $\\sigma_i = \\sqrt{2(1-\\rho_i)},\\ i=1\\dots 6$. This will favour a more compact\n", + "no-change cluster." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2b57369e", + "metadata": { + "id": "2b57369e" + }, + "outputs": [], + "source": [ + "# Standardize to no change sigmas.\n", + "sigma2s = ee.Image.constant([2*(1-x) for x in eval(rhos)])\n", + "im_imadstd = im_imad.divide(sigma2s.sqrt())\n", + "# Collect training data.\n", + "training = im_imadstd.sample(region=aoi, scale=scale, numPixels=50000)\n", + "# Train the clusterer.\n", + "clusterer = ee.Clusterer.wekaKMeans(4).train(training)\n", + "# Classify the standardized imad image.\n", + "result = im_imadstd.cluster(clusterer)" + ] + }, + { + "cell_type": "markdown", + "id": "22554e72", + "metadata": { + "id": "22554e72" + }, + "source": [ + "Here we display the four clusters overlayed onto the two Sentinel 2 images:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a0ba56eb", + "metadata": { + "id": "a0ba56eb" + }, + "outputs": [], + "source": [ + "M2 = geemap.Map()\n", + "M2.centerObject(aoi, 13)\n", + "display_ls(im1.select(visbands), M2, 'im1')\n", + "display_ls(im2.select(visbands), M2, 'im2')\n", + "cluster0 = result.updateMask(result.eq(0))\n", + "cluster1 = result.updateMask(result.eq(1))\n", + "cluster2 = result.updateMask(result.eq(2))\n", + "cluster3 = result.updateMask(result.eq(3))\n", + "palette = ['red', 'yellow', 'blue', 'black']\n", + "vis_params = {'min': 0, 'max': 3, 'palette': palette}\n", + "M2.addLayer(cluster0, vis_params, 'new clearcuts')\n", + "M2.addLayer(cluster1, vis_params, 'agriculture')\n", + "M2.addLayer(cluster2, vis_params, 'prior clearcuts')\n", + "M2.addLayer(cluster3, vis_params, 'no change')\n", + "\n", + "M2" + ] + }, + { + "cell_type": "markdown", + "id": "4d500eb8", + "metadata": { + "id": "4d500eb8" + }, + "source": [ + "### Interpretation\n", + "\n", + "**Cluster 0** (colored red in the preceding map) appears to classify the clear\n", + "cuts occurring over the observation period quite well.\n", + "\n", + "**Cluster 1** (yellow) marks changes in the agricultural fields and pastures.\n", + "\n", + "**Cluster 2** (blue) is more ambiguous but can be mainly associated with changes\n", + "in previously cleared forest (seasonal or new growth) as well as with some changes\n", + "in agricultural fields and in built up areas.\n", + "\n", + "**Cluster 3** (black) is no change." + ] + }, + { + "cell_type": "markdown", + "id": "a98f12ec", + "metadata": { + "id": "a98f12ec" + }, + "source": [ + "### Comparison with Dynamic World\n", + "\n", + "Google's recently released [Dynamic World](https://developers.google.com/earth-engine/tutorials/community/introduction-to-dynamic-world-pt-1)\n", + "dataset contains near real-time land use land cover predictions created from Sentinel-2\n", + "imagery for nine land use land cover classes including forest (trees class). It is interesting\n", + "to compare the loss in forest cover as determined from our iMAD/Cluster pipeline with the\n", + "Dynamic World tree map for the comparable time period. In the code snippet below, we gather\n", + "an image collection covering our observation period and simply mosaic them. The *mosaic()*\n", + "method composites the overlapping images according to their order in the collection (last on top),\n", + "which is what we want because the changes in tree cover are occurring over the entire period.\n", + "\n", + "Generally the agreement is excellent, although the iMAD change map registers a number of\n", + "small-area clear cuts missed in the Dynamic World map:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2ab45c3", + "metadata": { + "id": "b2ab45c3" + }, + "outputs": [], + "source": [ + "dyn = ee.ImageCollection('GOOGLE/DYNAMICWORLD/V1') \\\n", + " .filterDate('2021-06-01', '2022-06-30') \\\n", + " .filterBounds(aoi) \\\n", + " .select('label').mosaic()\n", + "# 'trees' class = class 1\n", + "dw = dyn.clip(aoi).updateMask(dyn.eq(1))\n", + "\n", + "M3 = geemap.Map()\n", + "M3.centerObject(aoi, 13)\n", + "display_ls(im1.select(visbands), M3, 'im1')\n", + "display_ls(im2.select(visbands), M3, 'im2')\n", + "M3.addLayer(dw, {'min': 0, 'max': 1, 'palette': ['black', 'green']}, 'dynamic world')\n", + "M3.addLayer(cluster0, vis_params, 'new clearcuts')\n", + "\n", + "M3" + ] + }, + { + "cell_type": "markdown", + "id": "13b0e6f2", + "metadata": { + "id": "13b0e6f2" + }, + "source": [ + "### Simple difference revisited\n", + "\n", + "In fact, K-means clustering of the simple difference image also gives a fairly good discrimination\n", + "of the clear cuts. This is because the NIR band is especially sensitive to all vegetation changes,\n", + "and is also only weakly correlated with the other 5 bands. However, close inspection indicates\n", + "that there are many more false positives, especially in agricultural fields, as well as in the reservoir." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1513b02f", + "metadata": { + "id": "1513b02f" + }, + "outputs": [], + "source": [ + "M4 = geemap.Map()\n", + "M4.centerObject(aoi, 13)\n", + "diff = im1.subtract(im2).select(visirbands)\n", + "training = diff.sample(region=aoi, scale=20, numPixels=50000)\n", + "clusterer = ee.Clusterer.wekaKMeans(4).train(training)\n", + "result1 = diff.cluster(clusterer)\n", + "cluster0d = result1.updateMask(result1.eq(0))\n", + "\n", + "display_ls(im1.select(visbands), M4, 'im1')\n", + "display_ls(im2.select(visbands), M4, 'im2')\n", + "M4.addLayer(cluster0d, {'min': 0, 'max': 3,\n", + " 'palette': ['orange', 'yellow', 'blue', 'black']}, 'clearcuts (diff)')\n", + "M4.addLayer(cluster0, vis_params, 'clearcuts (iMAD)')\n", + "\n", + "M4" + ] + }, + { + "cell_type": "markdown", + "id": "1d66e00a", + "metadata": { + "id": "1d66e00a" + }, + "source": [ + "### Deforestation quantified\n", + "\n", + "From the clear cuts class number 0, and using the *pixelArea()* function, we can extract the\n", + "total area cleared between June, 2021 and June, 2022 within the Landkreis Olpe, whereby we\n", + "exclude small areas covering less than 0.2 hectare:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0559159b", + "metadata": { + "id": "0559159b" + }, + "outputs": [], + "source": [ + "# Minimum contiguous area requirement (0.2 hectare).\n", + "contArea = cluster0.connectedPixelCount().selfMask()\n", + "# 0.2 hectare = 5 pixels @ 400 sq. meters.\n", + "mp = contArea.gte(ee.Number(5)).selfMask()\n", + "# Clear cuts in hectares.\n", + "pixelArea = mp.multiply(ee.Image.pixelArea()).divide(10000)\n", + "clearcutArea = pixelArea.reduceRegion(\n", + " reducer=ee.Reducer.sum(),\n", + " geometry=aoi,\n", + " scale=scale,\n", + " maxPixels=1e11)\n", + "ccA = clearcutArea.get('cluster').getInfo()\n", + "print(ccA, 'hectare')" + ] + }, + { + "cell_type": "markdown", + "id": "5177d831", + "metadata": { + "id": "5177d831" + }, + "source": [ + "The most recent [commercial forest inventory](https://www.it.nrw/itnrw) recorded for\n", + "the Landkreis Olpe as of December, 2019 was 40,178 hectare, so we have determined that,\n", + "allowing for further decreases in 2020 and the first half of 2021, more than 9.3% of\n", + "woodland was lost to drought/clearing within the time period measured.\n", + "\n", + "Finally, repeating the calculation with the clear cuts determined with the simple difference image:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef9eb1ea", + "metadata": { + "id": "ef9eb1ea" + }, + "outputs": [], + "source": [ + "contArea = cluster0d.connectedPixelCount().selfMask()\n", + "mp = contArea.gte(ee.Number(5)).selfMask()\n", + "pixelArea = mp.multiply(ee.Image.pixelArea()).divide(10000)\n", + "clearcutArea = pixelArea.reduceRegion(\n", + " reducer=ee.Reducer.sum(),\n", + " geometry=aoi,\n", + " scale=scale,\n", + " maxPixels=1e11)\n", + "ccA = clearcutArea.get('cluster').getInfo()\n", + "print(ccA, 'hectare')" + ] + }, + { + "cell_type": "markdown", + "id": "003dc804", + "metadata": { + "id": "003dc804" + }, + "source": [ + "thus overestimating the loss from clear cutting by about one third." + ] + }, + { + "cell_type": "markdown", + "id": "7a086c34", + "metadata": { + "id": "7a086c34" + }, + "source": [ + "## Summary\n", + "\n", + "In Part 2 of this tutorial, we have generalized the MAD transformation to an iterative scheme,\n", + "the iMAD algorithm. Then, we illustrated change detection with the algorithm by focussing a\n", + "particular application, namely detection of clear cutting of coniferous trees destroyed by\n", + "drought in an administrative district in Germany.\n", + "\n", + "While simple image comparison or differencing can be useful, the statistical transformations\n", + "implicit in the iMAD algorithm offer a more powerful means of analyzing and categorizing changes\n", + "in bitemporal image data. In Part 3, we will examine the use of iMAD for image calibration tasks,\n", + "giving some examples of relative radiometric normalization over an image sequence, as well as\n", + "harmonization of reflectances from different sensors." + ] + }, + { + "cell_type": "markdown", + "id": "6d6bb69e", + "metadata": { + "id": "6d6bb69e" + }, + "source": [ + "## Exercises\n", + "\n", + "1. Try another parameter set, or one of the other clusterers in the GEE arsenal to see if you can\n", + "improve on the above classification.\n", + "\n", + "2. In the discussion up till now we have not included the Sentinel-2 red edge bands. Repeat the\n", + "analysis with all 10 visual/infrared bands:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "07c3c37b", + "metadata": { + "id": "07c3c37b" + }, + "outputs": [], + "source": [ + "visirbands + rededgebands" + ] + }, + { + "cell_type": "markdown", + "id": "26cd6668", + "metadata": { + "id": "26cd6668" + }, + "source": [ + "3. Determine with the aid of a cloud-free S2 image from summer, 2020 the forest cover loss in the\n", + "district for the 2-year period ending June, 2022.\n", + "\n", + "\n", + "4. Urban and suburban sprawl accompany the growth of many large cities. If they are located in at\n", + "least partly forested areas, deforestation due to new housing and infrastructure development can\n", + "be very rapid and widespread. An extreme example is the city of Houston, Texas, where massive\n", + "encroachment on the surrounding coutryside is a recognized problem. Below is an area of interest\n", + "comprising Montgomery County, which encompasses heavily wooded areas to the north of the city,\n", + "and two cloud-free Sentinel-2 images from July, 2021 and June, 2022. Repeat the analysis with\n", + "the iMAD/Cluster pipeline to determine the loss of woodland within the County over that time\n", + "period. (Hint: Since a variety of mad-made changes occur in the scene, the interpretation of\n", + "unsupervised classification of the change image is critical.)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bfcd4420", + "metadata": { + "id": "bfcd4420" + }, + "outputs": [], + "source": [ + "# TIGER------: US Census Counties from the GEE Data Archive.\n", + "counties = ee.FeatureCollection('TIGER/2016/Counties')\n", + "filtered = counties.filter(ee.Filter.eq('NAMELSAD', 'Montgomery County'))\n", + "aois = filtered.geometry()\n", + "# There are many Montgomery Counties in USA, we want the 12th in the list.\n", + "aoi = ee.Geometry(aois.geometries().get(12))\n", + "im1, im2 = collect(aoi, '2021-07-01', '2021-07-30', '2022-06-01', '2022-06-30')\n", + "M5 = geemap.Map()\n", + "M5.centerObject(aoi, 10)\n", + "display_ls(im1.select(visbands), M5, 'Image 1')\n", + "display_ls(im2.select(visbands), M5, 'Image 2')\n", + "\n", + "M5" + ] + } + ], + "metadata": { + "colab": { + "provenance": [] + }, + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tutorials/imad-tutorial-pt2/index.md b/tutorials/imad-tutorial-pt2/index.md deleted file mode 100644 index 979e29b86..000000000 --- a/tutorials/imad-tutorial-pt2/index.md +++ /dev/null @@ -1,753 +0,0 @@ ---- -title: Change Detection in GEE - Part 2 The iMAD Algorithm -description: Iteratively re-weighted Multivariate Alteration Detection. -author: mortcanty -tags: change-detection, imad -date_published: 2023-03-24 ---- - - -{% include "earth-engine/_widgets/_notebook_buttons.html" with github_path="mortcanty/eeimad/blob/main/src/gee_change_detection_pt_2.ipynb" %} - -## Context - -In Part 1 of this tutorial, a statistical approach to detecting changes in pairs of -multispectral remote sensing images, called Multivariate Alteration Detection (MAD), was -explained. The MAD change images were obtained by first maximizing the correlations -between the original image bands and then subtracting one from the other. The resulting -difference bands, called MAD variates, contained the change information. They were shown -to be ordered by increasing variance and to be mutually uncorrelated. However, they were -also seen to be not easily interpretable. - -Now, in Part 2, we introduce an iteration scheme that performs the MAD transformation -exclusively on those pixels that mark areas in the images which have not physically -changed. This establishes a well-defined background of invariant pixels against which -to discriminate changes. - -## Preliminaries - -```python -# Enter your own export to assets path name here ----------- -EXPORT_PATH = 'projects/YOUR_GEE_PROJECT_NAME/assets/imad/' -# ------------------------------------------------ -``` - -```python -import ee - -ee.Authenticate(auth_mode='notebook') -ee.Initialize() -``` - -```python -# Import other packages used in the tutorial. -%matplotlib inline -import geemap -import numpy as np -import random, time -import matplotlib.pyplot as plt -from scipy.stats import norm, chi2 - -from pprint import pprint # for pretty printing -``` - -```python -# Routines from Part 1 - -def trunc(values, dec = 3): - '''Truncate a 1-D array to dec decimal places.''' - return np.trunc(values*10**dec)/(10**dec) - -# Display an image in a one percent linear stretch. -def display_ls(image, map, name, centered = False): - bns = image.bandNames().length().getInfo() - if bns == 3: - image = image.rename('B1', 'B2', 'B3') - pb_99 = ['B1_p99', 'B2_p99', 'B3_p99'] - pb_1 = ['B1_p1', 'B2_p1', 'B3_p1'] - img = ee.Image.rgb(image.select('B1'), image.select('B2'), image.select('B3')) - else: - image = image.rename('B1') - pb_99 = ['B1_p99'] - pb_1 = ['B1_p1'] - img = image.select('B1') - percentiles = image.reduceRegion(ee.Reducer.percentile([1, 99]), maxPixels=1e11) - mx = percentiles.values(pb_99) - if centered: - mn = ee.Array(mx).multiply(-1).toList() - else: - mn = percentiles.values(pb_1) - map.addLayer(img, {'min': mn, 'max': mx}, name) - -def collect(aoi, t1a ,t1b, t2a, t2b): - try: - im1 = ee.Image( ee.ImageCollection("COPERNICUS/S2_HARMONIZED") - .filterBounds(aoi) - .filterDate(ee.Date(t1a), ee.Date(t1b)) - .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo')) - .sort('CLOUDY_PIXEL_PERCENTAGE') - .first() - .clip(aoi) ) - im2 = ee.Image( ee.ImageCollection("COPERNICUS/S2_HARMONIZED") - .filterBounds(aoi) - .filterDate(ee.Date(t2a), ee.Date(t2b)) - .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo')) - .sort('CLOUDY_PIXEL_PERCENTAGE') - .first() - .clip(aoi) ) - timestamp = im1.date().format('E MMM dd HH:mm:ss YYYY') - print(timestamp.getInfo()) - timestamp = im2.date().format('E MMM dd HH:mm:ss YYYY') - print(timestamp.getInfo()) - return (im1, im2) - except Exception as e: - print('Error: %s'%e) - -def covarw(image, weights=None, scale=20, maxPixels=1e10): - '''Return the centered image and its weighted covariance matrix.''' - try: - geometry = image.geometry() - bandNames = image.bandNames() - N = bandNames.length() - if weights is None: - weights = image.constant(1) - weightsImage = image.multiply(ee.Image.constant(0)).add(weights) - means = image.addBands(weightsImage) \ - .reduceRegion(ee.Reducer.mean().repeat(N).splitWeights(), - scale = scale, - maxPixels = maxPixels) \ - .toArray() \ - .project([1]) - centered = image.toArray().subtract(means) - B1 = centered.bandNames().get(0) - b1 = weights.bandNames().get(0) - nPixels = ee.Number(centered.reduceRegion(ee.Reducer.count(), - scale=scale, - maxPixels=maxPixels).get(B1)) - sumWeights = ee.Number(weights.reduceRegion(ee.Reducer.sum(), - geometry=geometry, - scale=scale, - maxPixels=maxPixels).get(b1)) - covw = centered.multiply(weights.sqrt()) \ - .toArray() \ - .reduceRegion(ee.Reducer.centeredCovariance(), - geometry=geometry, - scale=scale, - maxPixels=maxPixels) \ - .get('array') - covw = ee.Array(covw).multiply(nPixels).divide(sumWeights) - return (centered.arrayFlatten([bandNames]), covw) - except Exception as e: - print('Error: %s'%e) - -def corr(cov): - '''Transfrom covariance matrix to correlation matrix.''' - # Diagonal matrix of inverse sigmas. - sInv = cov.matrixDiagonal().sqrt().matrixToDiag().matrixInverse() - # Transform. - corr = sInv.matrixMultiply(cov).matrixMultiply(sInv).getInfo() - # Truncate. - return [list(map(trunc, corr[i])) for i in range(len(corr))] - -def geneiv(C,B): - '''Return the eignvalues and eigenvectors of the generalized eigenproblem - C*X = lambda*B*X''' - try: - C = ee.Array(C) - B = ee.Array(B) - # Li = choldc(B)^-1 - Li = ee.Array(B.matrixCholeskyDecomposition().get('L')).matrixInverse() - # Solve symmetric, ordinary eigenproblem Li*C*Li^T*x = lambda*x - Xa = Li.matrixMultiply(C) \ - .matrixMultiply(Li.matrixTranspose()) \ - .eigen() - # Eigenvalues as a row vector. - lambdas = Xa.slice(1, 0, 1).matrixTranspose() - # Eigenvectors as columns. - X = Xa.slice(1, 1).matrixTranspose() - # Generalized eigenvectors as columns, Li^T*X - eigenvecs = Li.matrixTranspose().matrixMultiply(X) - return (lambdas, eigenvecs) - except Exception as e: - print('Error: %s'%e) - -def mad_run(image1, image2, scale=20): - '''The MAD transformation of two multiband images.''' - try: - image = image1.addBands(image2) - nBands = image.bandNames().length().divide(2) - centeredImage,covarArray = covarw(image,scale=scale) - bNames = centeredImage.bandNames() - bNames1 = bNames.slice(0,nBands) - bNames2 = bNames.slice(nBands) - centeredImage1 = centeredImage.select(bNames1) - centeredImage2 = centeredImage.select(bNames2) - s11 = covarArray.slice(0, 0, nBands).slice(1, 0, nBands) - s22 = covarArray.slice(0, nBands).slice(1, nBands) - s12 = covarArray.slice(0, 0, nBands).slice(1, nBands) - s21 = covarArray.slice(0, nBands).slice(1, 0, nBands) - c1 = s12.matrixMultiply(s22.matrixInverse()).matrixMultiply(s21) - b1 = s11 - c2 = s21.matrixMultiply(s11.matrixInverse()).matrixMultiply(s12) - b2 = s22 - # Solution of generalized eigenproblems. - lambdas, A = geneiv(c1, b1) - _, B = geneiv(c2, b2) - rhos = lambdas.sqrt().project(ee.List([1])) - # MAD variances. - sigma2s = rhos.subtract(1).multiply(-2).toList() - sigma2s = ee.Image.constant(sigma2s) - # Ensure sum of positive correlations between X and U is positive. - tmp = s11.matrixDiagonal().sqrt() - ones = tmp.multiply(0).add(1) - tmp = ones.divide(tmp).matrixToDiag() - s = tmp.matrixMultiply(s11).matrixMultiply(A).reduce(ee.Reducer.sum(),[0]).transpose() - A = A.matrixMultiply(s.divide(s.abs()).matrixToDiag()) - # Ensure positive correlation. - tmp = A.transpose().matrixMultiply(s12).matrixMultiply(B).matrixDiagonal() - tmp = tmp.divide(tmp.abs()).matrixToDiag() - B = B.matrixMultiply(tmp) - # Canonical and MAD variates as images. - centeredImage1Array = centeredImage1.toArray().toArray(1) - centeredImage2Array = centeredImage2.toArray().toArray(1) - U = ee.Image(A.transpose()).matrixMultiply(centeredImage1Array) \ - .arrayProject([0]) \ - .arrayFlatten([bNames2]) - V = ee.Image(B.transpose()).matrixMultiply(centeredImage2Array) \ - .arrayProject([0]) \ - .arrayFlatten([bNames2]) - MAD = U.subtract(V) - # Chi-square image. - Z = MAD.pow(2) \ - .divide(sigma2s) \ - .reduce(ee.Reducer.sum()) - return (U, V, MAD, Z) - except Exception as e: - print('Error: %s'%e) -``` - -## Iterative re-weighting - -Let's consider two images of the same scene acquired at different times under similar -conditions, similar to the Sentinel-2 images from Part 1, but for which no ground -reflectance changes have occurred. In this case, the only differences between the images -will be due to random effects such as instrument noise and atmospheric fluctuations. -Therefore, we can expect that the histogram of any difference component we generate will -be nearly Gaussian. Specifically, the MAD variates, which we have seen to be uncorrelated, -should follow a multivariate, zero-mean normal distribution with a diagonal covariance -matrix. - -$$ -\Sigma_M = \pmatrix{\sigma^2_{M_1} &0 &\cdots &0 \cr - 0 & \sigma^2_{M_2} &\cdots &0 \cr - \vdots &\vdots &\cdots &0 \cr - 0 & 0 &\cdots &\sigma^2_{M_N}}. -$$ - -Change observations would deviate more or less strongly from such a -distribution. We might therefore expect an improvement in the sensitivity -of the MAD transformation *if we can establish a better background of no -change against which to detect change.* This can be done in an iteration -scheme ([Nielsen 2007](https://www2.imm.dtu.dk/pubdb/pubs/4695-full.html)) -in which, when calculating the statistics for each successive iteration -of the MAD transformation, observations are weighted in some appropriate -fashion. - -Recall that the variable $Z$ represents the sum of the squares of the standardized MAD variates, - -$$ -Z = \sum_{i=1}^N\left({M_i\over \sigma_{M_i}}\right)^2, -$$ - -where $\sigma^2_{M_i}$ is given by Equation (8) in the first Tutorial, - -$$ -\sigma_{M_i}^2={\rm var}(U_i-V_i) = 2(1-\rho_i),\quad i=1\dots N, -$$ - -and $\rho_i = {\rm cov}(U_i,V_i)$. Then, since the no-change observations are expected -to be normally distributed and uncorrelated, basic statistical theory tells us -that the values of $Z$ corresponding to no-change observations should be *chi-square -distributed* with $N$ degrees of freedom. Let's check to what extent this is true for -the MAD variates that we have determined so far for the Landkreis Olpe scene. - -```python -# Landkreis Olpe. -aoi = ee.FeatureCollection('projects/google/imad_tutorial_aoi').geometry() - -visirbands = ['B2', 'B3', 'B4', 'B8', 'B11', 'B12'] -visbands = ['B2', 'B3', 'B4'] -rededgebands = ['B5', 'B6', 'B7', 'B8A'] - -# Collect the two Sentinel-2 images. -im1, im2 = collect(aoi, '2021-06-01', '2021-06-30', '2022-06-01', '2022-06-30') - -# Re-run MAD. -U, V, MAD, Z = mad_run(im1.select(visirbands), im2.select(visirbands), scale=20) - -# Plot histogram of Z. -hist = Z.reduceRegion(ee.Reducer.fixedHistogram(0, 50, 500), aoi, scale=20).get('sum').getInfo() -a = np.array(hist) -x = a[:, 0] # array of bucket edge positions -y = a[:, 1]/np.sum(a[:, 1]) # normalized array of bucket contents -plt.plot(x, y, '.', label='data') -# The chi-square distribution with 6 degrees of freedom. -plt.plot(x, chi2.pdf(x, 6)/10, '-r', label='Chi-square') -plt.legend() -plt.show() -``` - -Clearly not the case at all. Which is to be expected, since there are many change pixels in -the scene and we have made no attempt to discriminate them. - -In fact, it is easy to show that $Z$ is a *likelihood ratio test statistic* for change, see -the discussion of statistical hypothesis testing in the -[SAR Tutorial](https://developers.google.com/earth-engine/tutorials/community/detecting-changes-in-sentinel-1-imagery-pt-2) -and the discussion on pp.390-391 in -[Canty (2019)](https://www.taylorfrancis.com/books/image-analysis-classification-change-detection-remote-sensing-morton-john-canty/10.1201/9780429464348). -Under the hypothesis that no change has occurred, the test statistic $Z$ will follow, as we -said, a chi-square distribution. The so-called $p$-*value* is a measure of the extent to which -this is true. For an observation $z$ of the test statistic, the $p$-value is the probability -that any sample drawn from the chi-square distribution could be as large as $z$ or larger. -This is given by - -$$ -p(z) = 1-P_{\chi^2;N}(z),\quad 0 < p(z) < 1, -$$ - -where $P_{\chi^2;N}(z)$ is the cumulative chi-square probability distribution, i.e., the area -under the chi-square distribution up to the value $z$, and $p(z)$ is its complement. -All $p$-values are equally likely if no change has occurred at that pixel -location$^\star$, but **change will always be associated with small $p$-values**. -Therefore in order to eliminate the change observations from the MAD transformation, -the $p$-value itself can be used to weight each pixel before re-sampling the images to -determine the statistics for the next iteration. (This was the motivation for coding -a *weighted* covariance matrix routine in Part 1 earlier). The influence of the -change observations on the MAD transformation is thereby gradually reduced. Iteration -continues until some stopping criterion is met, such as lack of significant change in -the canonical correlations $\rho_i$. The whole procedure constitutes the *iMAD algorithm*. -It is implemented in the GEE Python API in the following cell: - -$\star$ Thus the $p$-value is *not* a no-change probability, a common misapprehension! -See again the [SAR Tutorial](https://developers.google.com/earth-engine/tutorials/community/detecting-changes-in-sentinel-1-imagery-pt-2). - -```python -# The iMAD code - -def chi2cdf(Z,df): - '''Chi-square cumulative distribution function with df degrees of freedom.''' - return ee.Image(Z.divide(2)).gammainc(ee.Number(df).divide(2)) - -def imad(current,prev): - '''Iterator function for iMAD.''' - done = ee.Number(ee.Dictionary(prev).get('done')) - return ee.Algorithms.If(done, prev, imad1(current, prev)) - -def imad1(current,prev): - '''Iteratively re-weighted MAD.''' - image = ee.Image(ee.Dictionary(prev).get('image')) - Z = ee.Image(ee.Dictionary(prev).get('Z')) - allrhos = ee.List(ee.Dictionary(prev).get('allrhos')) - nBands = image.bandNames().length().divide(2) - weights = chi2cdf(Z,nBands).subtract(1).multiply(-1) - scale = ee.Dictionary(prev).getNumber('scale') - niter = ee.Dictionary(prev).getNumber('niter') - # Weighted stacked image and weighted covariance matrix. - centeredImage, covarArray = covarw(image, weights, scale) - bNames = centeredImage.bandNames() - bNames1 = bNames.slice(0, nBands) - bNames2 = bNames.slice(nBands) - centeredImage1 = centeredImage.select(bNames1) - centeredImage2 = centeredImage.select(bNames2) - s11 = covarArray.slice(0, 0, nBands).slice(1, 0, nBands) - s22 = covarArray.slice(0, nBands).slice(1, nBands) - s12 = covarArray.slice(0, 0, nBands).slice(1, nBands) - s21 = covarArray.slice(0, nBands).slice(1, 0, nBands) - c1 = s12.matrixMultiply(s22.matrixInverse()).matrixMultiply(s21) - b1 = s11 - c2 = s21.matrixMultiply(s11.matrixInverse()).matrixMultiply(s12) - b2 = s22 - # Solution of generalized eigenproblems. - lambdas, A = geneiv(c1, b1) - _, B = geneiv(c2, b2) - rhos = lambdas.sqrt().project(ee.List([1])) - # Test for convergence. - lastrhos = ee.Array(allrhos.get(-1)) - done = rhos.subtract(lastrhos) \ - .abs() \ - .reduce(ee.Reducer.max(), ee.List([0])) \ - .lt(ee.Number(0.0001)) \ - .toList() \ - .get(0) - allrhos = allrhos.cat([rhos.toList()]) - # MAD variances. - sigma2s = rhos.subtract(1).multiply(-2).toList() - sigma2s = ee.Image.constant(sigma2s) - # Ensure sum of positive correlations between X and U is positive. - tmp = s11.matrixDiagonal().sqrt() - ones = tmp.multiply(0).add(1) - tmp = ones.divide(tmp).matrixToDiag() - s = tmp.matrixMultiply(s11).matrixMultiply(A).reduce(ee.Reducer.sum(), [0]).transpose() - A = A.matrixMultiply(s.divide(s.abs()).matrixToDiag()) - # Ensure positive correlation. - tmp = A.transpose().matrixMultiply(s12).matrixMultiply(B).matrixDiagonal() - tmp = tmp.divide(tmp.abs()).matrixToDiag() - B = B.matrixMultiply(tmp) - # Canonical and MAD variates. - centeredImage1Array = centeredImage1.toArray().toArray(1) - centeredImage2Array = centeredImage2.toArray().toArray(1) - U = ee.Image(A.transpose()).matrixMultiply(centeredImage1Array) \ - .arrayProject([0]) \ - .arrayFlatten([bNames1]) - V = ee.Image(B.transpose()).matrixMultiply(centeredImage2Array) \ - .arrayProject([0]) \ - .arrayFlatten([bNames2]) - iMAD = U.subtract(V) - # Chi-square image. - Z = iMAD.pow(2) \ - .divide(sigma2s) \ - .reduce(ee.Reducer.sum()) - return ee.Dictionary({'done': done, 'scale': scale, 'niter': niter.add(1), - 'image': image, 'allrhos': allrhos, 'Z': Z, 'iMAD': iMAD}) -``` - -The following code cell is a routine to run the iMAD algorithm as an export task, -avoiding memory and time limitations in the active runtime. The asset will be exported -to the location specified by the `EXPORT_PATH` variable defined earlier. It requires -about 130 MB of space and can take 15 to 20 minutes to complete. - -```python -# Run iMAD algorithm as export task - -def run_imad(aoi, image1, image2, assetFN, scale=20, maxiter=100): - try: - N = image1.bandNames().length().getInfo() - imadnames = ['iMAD'+str(i+1) for i in range(N)] - imadnames.append('Z') - # Maximum iterations. - inputlist = ee.List.sequence(1, maxiter) - first = ee.Dictionary({'done':0, - 'scale': scale, - 'niter': ee.Number(0), - 'image': image1.addBands(image2), - 'allrhos': [ee.List.sequence(1, N)], - 'Z': ee.Image.constant(0), - 'iMAD': ee.Image.constant(0)}) - # Iteration. - result = ee.Dictionary(inputlist.iterate(imad, first)) - # Retrieve results. - iMAD = ee.Image(result.get('iMAD')).clip(aoi) - rhos = ee.String.encodeJSON(ee.List(result.get('allrhos')).get(-1)) - Z = ee.Image(result.get('Z')) - niter = result.getNumber('niter') - # Export iMAD and Z as a single image, including rhos and number of iterations in properties. - iMAD_export = ee.Image.cat(iMAD, Z).rename(imadnames).set('rhos', rhos, 'niter', niter) - assetId = EXPORT_PATH + assetFN - assexport = ee.batch.Export.image.toAsset(iMAD_export, - description='assetExportTask', - assetId=assetId, scale=scale, maxPixels=1e10) - assexport.start() - print('Exporting iMAD to %s\n task id: %s'%(assetId, str(assexport.id))) - except Exception as e: - print('Error: %s'%e) -``` - -and here we run it on our two images: - -```python -run_imad(aoi, im1.select(visirbands), im2.select(visirbands), 'LandkreisOlpe') -``` - -After the export finishes, the number of iterations and the final canonical correlations -can be read from properties of the exported image: - -```python -im_imad = ee.Image(EXPORT_PATH + 'LandkreisOlpe').select(0, 1, 2, 3, 4, 5) -im_z = ee.Image(EXPORT_PATH + 'LandkreisOlpe').select(6).rename('Z') -niter = im_imad.get('niter').getInfo() -rhos = ee.List(im_imad.get('rhos')).getInfo() -print('iteratons: %i'%niter) -print('canonical correlations: %s'%rhos) -``` - -We got convergence after 28 iterations, and the correlations are very close to one -for the first canonical variates. It might now be interesting to check if the test -statistic $Z$ has the expected chi-square distribution when evaluated for the no -change pixels. To to eliminate the changes at the 10% significance level we set a -lower threshold of $\alpha = 0.1$ on the $p$-values (recall: small p-values signify change). - -```python -scale = 20 -# p-values image. -pval = chi2cdf(im_z, 6).subtract(1).multiply(-1).rename('pval') -# No-change mask (use p-values greater than 0.1). -noChangeMask = pval.gt(0.1) -hist = im_z.updateMask(noChangeMask).reduceRegion(ee.Reducer \ - .fixedHistogram(0, 50, 500), aoi, scale=scale, maxPixels=1e11) \ - .get('Z').getInfo() -a = np.array(hist) -x = a[:, 0] # array of bucket edge positions -y = a[:, 1]/np.sum(a[:, 1]) # normalized array of bucket contents -plt.plot(x, y, '.', label = 'data') -plt.plot(x, chi2.pdf(x, 6)/10, '-r', label='chi2(6)') -plt.legend() -plt.show() -``` - -Agreement is not perfect, but the plot is certainly closer to the ideal chi-square distribution -after iteration than after the single MAD transformation. So let's display the *iMAD* image: - -```python -M1 = geemap.Map() -M1.centerObject(aoi, 11) -display_ls(im1.select(visbands), M1, 'im1') -display_ls(im2.select(visbands), M1, 'im2') -display_ls(im_imad.select('iMAD1', 'iMAD2', 'iMAD3'), M1, 'iMAD123', True) - -M1 -``` - -Gray pixels point to no change, while the wide range of color in the iMAD variates -indicates a good discrimination of the types of change occurring. - -**Aside:** We are of course primarily interested in extracting the changes in the iMAD -image, especially those which mark clear cutting, and we'll come back to them in a moment. -However, now that what we think the no change pixels have been isolated, we could -also perform a regression analysis on them to determine how well the radiometric parameters -of the two Sentinel-2 acquisitions compare. If surface rather than top of atmosphere (TOA) -reflectance images had been used for the example, we would expect a good match, i.e., a -slope close to one and an intercept close to zero at all spectral wavelengths. In general, -for uncalibrated images, this will not be the case. In that event, the regression -coefficients can be applied to normalize one image (the target, say) to the other -(the reference). This might be desirable for tracing features such as NDVI indices over -a time series of acquisitions when the images have not been reduced to surface reflectances, -see e.g., [Gan et al. (2021)](https://ieeexplore.ieee.org/document/9392311), or indeed for -*harmonizing* the data from two different sensors of the same family such as Landsat 7 with -Landsat 8. These topics will be the subject of Part 3. - -But now let's look in more detail at the changes in the Landkreis Olpe scene. - -### Clustering - -To better interpret the change image, we can attempt an unsupervised classification. -We'll see if we can get away with an ordinary K-means clusterer and a simple Euclidean distance -measure for the complete iMAD image. We choose the number of clusters as 4 and leave all 12(!) -other input parameters of the *wekaKmeans()* clusterer at their default values. We will also -first standardize the iMAD image by dividing by the square root of the variances of the no-change -pixels, $\sigma_i = \sqrt{2(1-\rho_i)},\ i=1\dots 6$. This will favour a more compact -no-change cluster. - -```python -# Standardize to no change sigmas. -sigma2s = ee.Image.constant([2*(1-x) for x in eval(rhos)]) -im_imadstd = im_imad.divide(sigma2s.sqrt()) -# Collect training data. -training = im_imadstd.sample(region=aoi, scale=scale, numPixels=50000) -# Train the clusterer. -clusterer = ee.Clusterer.wekaKMeans(4).train(training) -# Classify the standardized imad image. -result = im_imadstd.cluster(clusterer) -``` - -Here we display the four clusters overlaid onto the two Sentinel 2 images: - -```python -M2 = geemap.Map() -M2.centerObject(aoi, 13) -display_ls(im1.select(visbands), M2, 'im1') -display_ls(im2.select(visbands), M2, 'im2') -cluster0 = result.updateMask(result.eq(0)) -cluster1 = result.updateMask(result.eq(1)) -cluster2 = result.updateMask(result.eq(2)) -cluster3 = result.updateMask(result.eq(3)) -palette = ['red', 'yellow', 'blue', 'black'] -vis_params = {'min': 0, 'max': 3, 'palette': palette} -M2.addLayer(cluster0, vis_params, 'new clearcuts') -M2.addLayer(cluster1, vis_params, 'agriculture') -M2.addLayer(cluster2, vis_params, 'prior clearcuts') -M2.addLayer(cluster3, vis_params, 'no change') - -M2 -``` - -### Interpretation - -**Cluster 0** (colored red in the preceding map) appears to classify the clear -cuts occurring over the observation period quite well. - -**Cluster 1** (yellow) marks changes in the agricultural fields and pastures. - -**Cluster 2** (blue) is more ambiguous but can be mainly associated with changes -in previously cleared forest (seasonal or new growth) as well as with some changes -in agricultural fields and in built up areas. - -**Cluster 3** (black) is no change. - -### Comparison with Dynamic World - -Google's recently released [Dynamic World](https://developers.google.com/earth-engine/tutorials/community/introduction-to-dynamic-world-pt-1) -dataset contains near real-time land use land cover predictions created from Sentinel-2 -imagery for nine land use land cover classes including forest (trees class). It is interesting -to compare the loss in forest cover as determined from our iMAD/Cluster pipeline with the -Dynamic World tree map for the comparable time period. In the code snippet below, we gather -an image collection covering our observation period and simply mosaic them. The *mosaic()* -method composites the overlapping images according to their order in the collection (last on top), -which is what we want because the changes in tree cover are occurring over the entire period. - -Generally the agreement is excellent, although the iMAD change map registers a number of -small-area clear cuts missed in the Dynamic World map: - -```python -dyn = ee.ImageCollection('GOOGLE/DYNAMICWORLD/V1') \ - .filterDate('2021-06-01', '2022-06-30') \ - .filterBounds(aoi) \ - .select('label').mosaic() -# 'trees' class = class 1 -dw = dyn.clip(aoi).updateMask(dyn.eq(1)) - -M3 = geemap.Map() -M3.centerObject(aoi, 13) -display_ls(im1.select(visbands), M3, 'im1') -display_ls(im2.select(visbands), M3, 'im2') -M3.addLayer(dw, {'min': 0, 'max': 1, 'palette': ['black', 'green']}, 'dynamic world') -M3.addLayer(cluster0, vis_params, 'new clearcuts') - -M3 -``` - -### Simple difference revisited - -In fact, K-means clustering of the simple difference image also gives a fairly good discrimination -of the clear cuts. This is because the NIR band is especially sensitive to all vegetation changes, -and is also only weakly correlated with the other 5 bands. However, close inspection indicates -that there are many more false positives, especially in agricultural fields, as well as in the reservoir. - -```python -M4 = geemap.Map() -M4.centerObject(aoi, 13) -diff = im1.subtract(im2).select(visirbands) -training = diff.sample(region=aoi, scale=20, numPixels=50000) -clusterer = ee.Clusterer.wekaKMeans(4).train(training) -result1 = diff.cluster(clusterer) -cluster0d = result1.updateMask(result1.eq(0)) - -display_ls(im1.select(visbands), M4, 'im1') -display_ls(im2.select(visbands), M4, 'im2') -M4.addLayer(cluster0d, {'min': 0, 'max': 3, - 'palette': ['orange', 'yellow', 'blue', 'black']}, 'clearcuts (diff)') -M4.addLayer(cluster0, vis_params, 'clearcuts (iMAD)') - -M4 -``` - -### Deforestation quantified - -From the clear cuts class number 0, and using the *pixelArea()* function, we can extract the -total area cleared between June, 2021 and June, 2022 within the Landkreis Olpe, whereby we -exclude small areas covering less than 0.2 hectare: - -```python -# Minimum contiguous area requirement (0.2 hectare). -contArea = cluster0.connectedPixelCount().selfMask() -# 0.2 hectare = 5 pixels @ 400 sq. meters. -mp = contArea.gte(ee.Number(5)).selfMask() -# Clear cuts in hectares. -pixelArea = mp.multiply(ee.Image.pixelArea()).divide(10000) -clearcutArea = pixelArea.reduceRegion( - reducer=ee.Reducer.sum(), - geometry=aoi, - scale=scale, - maxPixels=1e11) -ccA = clearcutArea.get('cluster').getInfo() -print(ccA, 'hectare') -``` - -The most recent [commercial forest inventory](https://www.it.nrw/itnrw) recorded for -the Landkreis Olpe as of December, 2019 was 40,178 hectare, so we have determined that, -allowing for further decreases in 2020 and the first half of 2021, more than 9.3% of -woodland was lost to drought/clearing within the time period measured. - -Finally, repeating the calculation with the clear cuts determined with the simple difference image: - -```python -contArea = cluster0d.connectedPixelCount().selfMask() -mp = contArea.gte(ee.Number(5)).selfMask() -pixelArea = mp.multiply(ee.Image.pixelArea()).divide(10000) -clearcutArea = pixelArea.reduceRegion( - reducer=ee.Reducer.sum(), - geometry=aoi, - scale=scale, - maxPixels=1e11) -ccA = clearcutArea.get('cluster').getInfo() -print(ccA, 'hectare') -``` - -thus overestimating the loss from clear cutting by about one third. - -## Summary - -In Part 2 of this tutorial, we have generalized the MAD transformation to an iterative scheme, -the iMAD algorithm. Then, we illustrated change detection with the algorithm by focussing a -particular application, namely detection of clear cutting of coniferous trees destroyed by -drought in an administrative district in Germany. - -While simple image comparison or differencing can be useful, the statistical transformations -implicit in the iMAD algorithm offer a more powerful means of analyzing and categorizing changes -in bitemporal image data. In Part 3, we will examine the use of iMAD for image calibration tasks, -giving some examples of relative radiometric normalization over an image sequence, as well as -harmonization of reflectances from different sensors. - -## Exercises - -1. Try another parameter set, or one of the other clusterers in the GEE arsenal to see if you can -improve on the above classification. - -2. In the discussion up till now we have not included the Sentinel-2 red edge bands. Repeat the -analysis with all 10 visual/infrared bands: - -```python -visirbands + rededgebands -``` - -3. Determine with the aid of a cloud-free S2 image from summer, 2020 the forest cover loss in the -district for the 2-year period ending June, 2022. - -4. Urban and suburban sprawl accompany the growth of many large cities. If they are located in at -least partly forested areas, deforestation due to new housing and infrastructure development can -be very rapid and widespread. An extreme example is the city of Houston, Texas, where massive -encroachment on the surrounding coutryside is a recognized problem. Below is an area of interest -comprising Montgomery County, which encompasses heavily wooded areas to the north of the city, -and two cloud-free Sentinel-2 images from July, 2021 and June, 2021. Repeat the analysis with -the iMAD/Cluster pipeline to determine the loss of woodland within the County over that time -period. (Hint: Since a variety of mad-made changes occur in the scene, the interpretation of -unsupervised classification of the change image is critical.) - -```python -# TIGER------: US Census Counties from the GEE Data Archive. -counties = ee.FeatureCollection('TIGER/2016/Counties') -filtered = counties.filter(ee.Filter.eq('NAMELSAD', 'Montgomery County')) -aois = filtered.geometry() -# There are many Montgomery Counties in USA, we want the 12th in the list. -aoi = ee.Geometry(aois.geometries().get(12)) -im1, im2 = collect(aoi, '2021-07-01', '2021-07-30', '2022-06-01', '2022-06-30') -M5 = geemap.Map() -M5.centerObject(aoi, 10) -display_ls(im1.select(visbands), M5, 'Image 1') -display_ls(im2.select(visbands), M5, 'Image 2') - -M5 -``` \ No newline at end of file diff --git a/tutorials/imad-tutorial-pt3/index.ipynb b/tutorials/imad-tutorial-pt3/index.ipynb new file mode 100644 index 000000000..e7a1c5e72 --- /dev/null +++ b/tutorials/imad-tutorial-pt3/index.ipynb @@ -0,0 +1,1039 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "zandql6BLXYn" + }, + "outputs": [], + "source": [ + "#@title Copyright 2023 The Earth Engine Community Authors { display-mode: \"form\" }\n", + "#\n", + "# Licensed under the Apache License, Version 2.0 (the \"License\");\n", + "# you may not use this file except in compliance with the License.\n", + "# You may obtain a copy of the License at\n", + "#\n", + "# https://www.apache.org/licenses/LICENSE-2.0\n", + "#\n", + "# Unless required by applicable law or agreed to in writing, software\n", + "# distributed under the License is distributed on an \"AS IS\" BASIS,\n", + "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n", + "# See the License for the specific language governing permissions and\n", + "# limitations under the License." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "pJrMq0aNrF-u" + }, + "source": [ + "# Change Detection in GEE - Part 3 Relative Radiometric Normalization and Harmonization\n", + " Authors: mortcanty" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "xsCfTJVeLXYp" + }, + "source": [ + "Test of pixel intensities across different image acquisitions requires that the received signals have similar radiometric scales.\n", + "As we saw in [Part 2](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt2) of the present tutorial, provided the iMAD algorithm converges satisfactorily, it will isolate the pixels\n", + "which have not significantly changed in the two input images presented to it. A regression analysis on these no-change\n", + "pixels can then determine how well the radiometric parameters of the two acquisitions compare.\n", + "If, for instance, images that were corrected to absolute surface reflectance are examined in this way, we would expect a good match, i.e., a\n", + "slope close to one and an intercept close to zero at all spectral wavelengths. In general,\n", + "for uncalibrated or poorly corrected images, this will not be the case. Then, the regression\n", + "coefficients can be applied to normalize one image (the target, say) to the other\n", + "(the reference). This is a prerequisite, for example, for tracing features such as NDVI indices over\n", + "a time series when the images have not been reduced to surface reflectances,\n", + "see e.g., [Gan et al. (2021)](https://ieeexplore.ieee.org/document/9392311), or indeed for\n", + "*harmonizing* the data from two different sensors of the same family, such as Landsat 8 with\n", + "Landsat 9. In this final part of the GEE Change Detection Tutorial we'll have a closer look at this approach, which we'll refer to here as *relative radiometric normalization*. For more detailed treatment, see also [Canty and Nielsen (2008)](https://www.sciencedirect.com/science/article/abs/pii/S0034425707003495), [Canty et al. (2004)](https://www.sciencedirect.com/science/article/abs/pii/S0034425704001208), and [Schroeder et al. (2006)](https://www.sciencedirect.com/science/article/abs/pii/S0034425706001179). A different (but similar) method involving identification of so-called *pseudoinvariant features*\n", + "is explained in [another GEE community tutorial](https://developers.google.com/earth-engine/tutorials/community/pseudo-invariant-feature-matching) in this series." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "A_dVj6pBLXYp" + }, + "source": [ + "## Preliminaries" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "x9Fw9beDLXYp", + "tags": [] + }, + "outputs": [], + "source": [ + "# Enter your own export to assets path name here -----------\n", + "EXPORT_PATH = 'projects/YOUR_GEE_PROJECT_NAME/assets/imad/'\n", + "# ------------------------------------------------" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "4w4PeArZLXYp", + "tags": [] + }, + "outputs": [], + "source": [ + "import ee\n", + "\n", + "ee.Authenticate(auth_mode='notebook')\n", + "ee.Initialize()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "1NQTCj6TLXYq", + "tags": [] + }, + "outputs": [], + "source": [ + "# Import other packages used in the tutorial\n", + "%matplotlib inline\n", + "import geemap\n", + "import numpy as np\n", + "import random, time\n", + "import matplotlib.pyplot as plt\n", + "from scipy.stats import norm, chi2\n", + "\n", + "from pprint import pprint # for pretty printing" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "LPxfmiKVLXYq" + }, + "outputs": [], + "source": [ + "#@title Routines from Parts 1 and 2\n", + "\n", + "def trunc(values, dec = 3):\n", + " '''Truncate a 1-D array to dec decimal places.'''\n", + " return np.trunc(values*10**dec)/(10**dec)\n", + "\n", + "# Display an image in a one percent linear stretch.\n", + "def display_ls(image, map, name, centered = False):\n", + " bns = image.bandNames().length().getInfo()\n", + " if bns == 3:\n", + " image = image.rename('B1', 'B2', 'B3')\n", + " pb_99 = ['B1_p99', 'B2_p99', 'B3_p99']\n", + " pb_1 = ['B1_p1', 'B2_p1', 'B3_p1']\n", + " img = ee.Image.rgb(\n", + " image.select('B1'), image.select('B2'), image.select('B3'))\n", + " else:\n", + " image = image.rename('B1')\n", + " pb_99 = ['B1_p99']\n", + " pb_1 = ['B1_p1']\n", + " img = image.select('B1')\n", + " percentiles = image.reduceRegion(\n", + " ee.Reducer.percentile([1, 99]), maxPixels=1e11)\n", + " mx = percentiles.values(pb_99)\n", + " if centered:\n", + " mn = ee.Array(mx).multiply(-1).toList()\n", + " else:\n", + " mn = percentiles.values(pb_1)\n", + " map.addLayer(img, {'min': mn, 'max': mx}, name)\n", + "\n", + "def covarw(image, weights=None, scale=20, maxPixels=1e10):\n", + " '''Return the centered image and its weighted covariance matrix.'''\n", + " try:\n", + " geometry = image.geometry()\n", + " bandNames = image.bandNames()\n", + " N = bandNames.length()\n", + " if weights is None:\n", + " weights = image.constant(1)\n", + " weightsImage = image.multiply(ee.Image.constant(0)).add(weights)\n", + " means = image.addBands(weightsImage) \\\n", + " .reduceRegion(ee.Reducer.mean().repeat(N).splitWeights(),\n", + " scale=scale,\n", + " maxPixels=maxPixels) \\\n", + " .toArray() \\\n", + " .project([1])\n", + " centered = image.toArray().subtract(means)\n", + " B1 = centered.bandNames().get(0)\n", + " b1 = weights.bandNames().get(0)\n", + " nPixels = ee.Number(centered.reduceRegion(ee.Reducer.count(),\n", + " scale=scale,\n", + " maxPixels=maxPixels).get(B1))\n", + " sumWeights = ee.Number(weights.reduceRegion(ee.Reducer.sum(),\n", + " geometry=geometry,\n", + " scale=scale,\n", + " maxPixels=maxPixels).get(b1))\n", + " covw = centered.multiply(weights.sqrt()) \\\n", + " .toArray() \\\n", + " .reduceRegion(ee.Reducer.centeredCovariance(),\n", + " geometry=geometry,\n", + " scale=scale,\n", + " maxPixels=maxPixels) \\\n", + " .get('array')\n", + " covw = ee.Array(covw).multiply(nPixels).divide(sumWeights)\n", + " return (centered.arrayFlatten([bandNames]), covw)\n", + " except Exception as e:\n", + " print('Error: %s'%e)\n", + "\n", + "def corr(cov):\n", + " '''Transform covariance matrix to correlation matrix.'''\n", + " # Diagonal matrix of inverse sigmas.\n", + " sInv = cov.matrixDiagonal().sqrt().matrixToDiag().matrixInverse()\n", + " # Transform.\n", + " corr = sInv.matrixMultiply(cov).matrixMultiply(sInv).getInfo()\n", + " # Truncate.\n", + " return [list(map(trunc, corr[i])) for i in range(len(corr))]\n", + "\n", + "def geneiv(C,B):\n", + " '''Return the eigenvalues and eigenvectors of the generalized eigenproblem\n", + " C*X = lambda*B*X'''\n", + " try:\n", + " C = ee.Array(C)\n", + " B = ee.Array(B)\n", + " # Li = choldc(B)^-1\n", + " Li = ee.Array(B.matrixCholeskyDecomposition().get('L')).matrixInverse()\n", + " # Solve symmetric, ordinary eigenproblem Li*C*Li^T*x = lambda*x\n", + " Xa = Li.matrixMultiply(C) \\\n", + " .matrixMultiply(Li.matrixTranspose()) \\\n", + " .eigen()\n", + " # Eigenvalues as a row vector.\n", + " lambdas = Xa.slice(1, 0, 1).matrixTranspose()\n", + " # Eigenvectors as columns.\n", + " X = Xa.slice(1, 1).matrixTranspose()\n", + " # Generalized eigenvectors as columns, Li^T*X\n", + " eigenvecs = Li.matrixTranspose().matrixMultiply(X)\n", + " return (lambdas, eigenvecs)\n", + " except Exception as e:\n", + " print('Error: %s'%e)\n", + "\n", + "# Collect a Sentinel-2 image pair.\n", + "def collect(aoi, t1a ,t1b, t2a, t2b, coll1, coll2 = None):\n", + " try:\n", + " if coll2 == None:\n", + " coll2 = coll1\n", + " im1 = ee.Image( ee.ImageCollection(coll1)\n", + " .filterBounds(aoi)\n", + " .filterDate(ee.Date(t1a), ee.Date(t1b))\n", + " .filter(ee.Filter.contains(rightValue=aoi, leftField='.geo'))\n", + " .sort('CLOUDY_PIXEL_PERCENTAGE')\n", + " .sort('CLOUD_COVER_LAND')\n", + " .first()\n", + " .clip(aoi) )\n", + " im2 = ee.Image( ee.ImageCollection(coll2)\n", + " .filterBounds(aoi)\n", + " .filterDate(ee.Date(t2a), ee.Date(t2b))\n", + " .filter(ee.Filter.contains(rightValue=aoi, leftField='.geo'))\n", + " .sort('CLOUDY_PIXEL_PERCENTAGE')\n", + " .sort('CLOUD_COVER_LAND')\n", + " .first()\n", + " .clip(aoi) )\n", + " timestamp = im1.date().format('E MMM dd HH:mm:ss YYYY')\n", + " print(timestamp.getInfo())\n", + " timestamp = im2.date().format('E MMM dd HH:mm:ss YYYY')\n", + " print(timestamp.getInfo())\n", + " return (im1, im2)\n", + " except Exception as e:\n", + " print('Error: %s'%e)\n", + "\n", + "#@title The iMAD code\n", + "def chi2cdf(Z,df):\n", + " '''Chi-square cumulative distribution function with df degrees of freedom.'''\n", + " return ee.Image(Z.divide(2)).gammainc(ee.Number(df).divide(2))\n", + "\n", + "def imad(current,prev):\n", + " '''Iterator function for iMAD.'''\n", + " done = ee.Number(ee.Dictionary(prev).get('done'))\n", + " return ee.Algorithms.If(done, prev, imad1(current, prev))\n", + "\n", + "def imad1(current,prev):\n", + " '''Iteratively re-weighted MAD.'''\n", + " image = ee.Image(ee.Dictionary(prev).get('image'))\n", + " Z = ee.Image(ee.Dictionary(prev).get('Z'))\n", + " allrhos = ee.List(ee.Dictionary(prev).get('allrhos'))\n", + " nBands = image.bandNames().length().divide(2)\n", + " weights = chi2cdf(Z,nBands).subtract(1).multiply(-1)\n", + " scale = ee.Dictionary(prev).getNumber('scale')\n", + " niter = ee.Dictionary(prev).getNumber('niter')\n", + " # Weighted stacked image and weighted covariance matrix.\n", + " centeredImage, covarArray = covarw(image, weights, scale)\n", + " bNames = centeredImage.bandNames()\n", + " bNames1 = bNames.slice(0, nBands)\n", + " bNames2 = bNames.slice(nBands)\n", + " centeredImage1 = centeredImage.select(bNames1)\n", + " centeredImage2 = centeredImage.select(bNames2)\n", + " s11 = covarArray.slice(0, 0, nBands).slice(1, 0, nBands)\n", + " s22 = covarArray.slice(0, nBands).slice(1, nBands)\n", + " s12 = covarArray.slice(0, 0, nBands).slice(1, nBands)\n", + " s21 = covarArray.slice(0, nBands).slice(1, 0, nBands)\n", + " c1 = s12.matrixMultiply(s22.matrixInverse()).matrixMultiply(s21)\n", + " b1 = s11\n", + " c2 = s21.matrixMultiply(s11.matrixInverse()).matrixMultiply(s12)\n", + " b2 = s22\n", + " # Solution of generalized eigenproblems.\n", + " lambdas, A = geneiv(c1, b1)\n", + " _, B = geneiv(c2, b2)\n", + " rhos = lambdas.sqrt().project(ee.List([1]))\n", + " # Test for convergence.\n", + " lastrhos = ee.Array(allrhos.get(-1))\n", + " done = rhos.subtract(lastrhos) \\\n", + " .abs() \\\n", + " .reduce(ee.Reducer.max(), ee.List([0])) \\\n", + " .lt(ee.Number(0.0001)) \\\n", + " .toList() \\\n", + " .get(0)\n", + " allrhos = allrhos.cat([rhos.toList()])\n", + " # MAD variances.\n", + " sigma2s = rhos.subtract(1).multiply(-2).toList()\n", + " sigma2s = ee.Image.constant(sigma2s)\n", + " # Ensure sum of positive correlations between X and U is positive.\n", + " tmp = s11.matrixDiagonal().sqrt()\n", + " ones = tmp.multiply(0).add(1)\n", + " tmp = ones.divide(tmp).matrixToDiag()\n", + " s = tmp.matrixMultiply(s11).matrixMultiply(A).reduce(ee.Reducer.sum(), [0]).transpose()\n", + " A = A.matrixMultiply(s.divide(s.abs()).matrixToDiag())\n", + " # Ensure positive correlation.\n", + " tmp = A.transpose().matrixMultiply(s12).matrixMultiply(B).matrixDiagonal()\n", + " tmp = tmp.divide(tmp.abs()).matrixToDiag()\n", + " B = B.matrixMultiply(tmp)\n", + " # Canonical and MAD variates.\n", + " centeredImage1Array = centeredImage1.toArray().toArray(1)\n", + " centeredImage2Array = centeredImage2.toArray().toArray(1)\n", + " U = ee.Image(A.transpose()).matrixMultiply(centeredImage1Array) \\\n", + " .arrayProject([0]) \\\n", + " .arrayFlatten([bNames1])\n", + " V = ee.Image(B.transpose()).matrixMultiply(centeredImage2Array) \\\n", + " .arrayProject([0]) \\\n", + " .arrayFlatten([bNames2])\n", + " iMAD = U.subtract(V)\n", + " # Chi-square image.\n", + " Z = iMAD.pow(2) \\\n", + " .divide(sigma2s) \\\n", + " .reduce(ee.Reducer.sum())\n", + " return ee.Dictionary({'done': done, 'scale': scale, 'niter': niter.add(1),\n", + " 'image': image, 'allrhos': allrhos, 'Z': Z, 'iMAD': iMAD})\n", + "\n", + "#@title Run iMAD algorithm as export task\n", + "def run_imad(aoi, image1, image2, assetFN, scale=10, maxiter=100):\n", + " try:\n", + " N = image1.bandNames().length().getInfo()\n", + " imadnames = ['iMAD'+str(i+1) for i in range(N)]\n", + " imadnames.append('Z')\n", + " # Maximum iterations.\n", + " inputlist = ee.List.sequence(1, maxiter)\n", + " first = ee.Dictionary({'done':0,\n", + " 'scale': scale,\n", + " 'niter': ee.Number(0),\n", + " 'image': image1.addBands(image2),\n", + " 'allrhos': [ee.List.sequence(1, N)],\n", + " 'Z': ee.Image.constant(0),\n", + " 'iMAD': ee.Image.constant(0)})\n", + " # Iteration.\n", + " result = ee.Dictionary(inputlist.iterate(imad, first))\n", + " # Retrieve results.\n", + " iMAD = ee.Image(result.get('iMAD')).clip(aoi)\n", + " rhos = ee.String.encodeJSON(ee.List(result.get('allrhos')).get(-1))\n", + " Z = ee.Image(result.get('Z'))\n", + " niter = result.getNumber('niter')\n", + " # Export iMAD and Z as a singe image, including rhos and number of iterations in properties.\n", + " iMAD_export = ee.Image.cat(iMAD, Z).rename(imadnames).set('rhos', rhos, 'niter', niter)\n", + " assetId = EXPORT_PATH + assetFN\n", + " assexport = ee.batch.Export.image.toAsset(iMAD_export,\n", + " description='assetExportTask',\n", + " assetId=assetId, scale=scale, maxPixels=1e10)\n", + " assexport.start()\n", + " print('Exporting iMAD to %s\\n task id: %s'%(assetId, str(assexport.id)))\n", + " except Exception as e:\n", + " print('Error: %s'%e)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "1EKu26_NLXYr" + }, + "source": [ + "### Relative radiometric normalization\n", + "\n", + "We will illustrate the idea with the German administrative district scene used in the first two parts of this tutorial ([Part 1](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt1), [Part 2](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt2)), the Landkreis Olpe. Once we have isolated what we think to be the no-change pixels for two different Sentinel-2 acquisitions in the aoi, we can perform a regression analysis on them to determine how well their radiometric parameters compare *relative* to the identified no-change pixels in the two scenes. Clearly, both variables involved have measurement uncertainty associated with them. In fact, which acquisition is termed reference and which is termed target data is arbitrary. Therefore it is preferable to use *orthogonal linear regression* rather than ordinary linear regression. The former method treats the data symmetrically whereas the latter assumes uncertainty only in the dependent variable.\n", + "\n", + "Orthogonal linear regression on two variables is explained in detail in the Appendix of [Canty et al. (2004)](https://www.sciencedirect.com/science/article/abs/pii/S0034425704001208), but can be briefly summarized as follows: let the no-change observations in a given spectral band of the reference and target images be $X_i$ and $Y_i, \\ i=1\\dots n$, respectively. Then define the $n\\times 2$ *data matrix*\n", + "\n", + "$$\n", + "d = \\begin{pmatrix} X_1 \u0026 Y_1\\cr X_2 \u0026 Y_2 \\cr \\vdots \u0026 \\vdots \\cr X_n \u0026 Y_n\\end{pmatrix},\n", + "$$\n", + "\n", + "and calculate from it the mean values $\\bar X,\\ \\bar Y$ and the covariance matrix\n", + "\n", + "$$\n", + "\\Sigma =\\begin{pmatrix} \\sigma^2_X \u0026 \\sigma_{XY} \\cr \\sigma_{YX}\\ \u0026 \\sigma^2_Y\\end{pmatrix}.\n", + "$$\n", + "\n", + "The best estimates for the slope $b$ and intercept $a$ of the orthogonal regression line\n", + "\n", + "$$\n", + "Y = bX+a\n", + "$$\n", + "\n", + "are then given by\n", + "\n", + "$$\n", + "\\hat b = {\\sigma^2_X -\\sigma^2_Y + \\sqrt{(\\sigma^2_X -\\sigma^2_Y)^2+4\\sigma^2_{XY}}\\over 2\\sigma_{XY}}, \\quad \\hat a = \\bar Y - b\\bar X.\n", + "$$\n", + "\n", + "As a simple quality criterion we can use the correlation coefficient\n", + "\n", + "$$\\rho = {\\sigma_{XY}\\over \\sigma_X\\sigma_Y}$$\n", + "\n", + "which has value one for a perfect fit.\n", + "\n", + "The following cell codes a function to iterate orthogonal regression over all of the spectral bands of an image pair." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "3bKT0fvXLXYr", + "tags": [] + }, + "outputs": [], + "source": [ + "def orthoregress(current, prev):\n", + " '''\n", + " Iterator function for orthogonal regression\n", + " '''\n", + " k = ee.Number(current)\n", + " prev = ee.Dictionary(prev)\n", + " # Image is concatenation of reference and target.\n", + " image = ee.Image(prev.get('image'))\n", + " coeffs = ee.List(prev.get('coeffs'))\n", + " N = image.bandNames().length().divide(2)\n", + " # Data matrix.\n", + " d = image.select(k, k.add(N))\n", + " means = d.reduceRegion(ee.Reducer.mean(), scale=10, maxPixels=1e10) \\\n", + " .toArray() \\\n", + " .project([0])\n", + " Xm = means.get([0])\n", + " Ym = means.get([1])\n", + " Sigma = ee.Array(d.toArray() \\\n", + " .reduceRegion(ee.Reducer.covariance(), scale=10, maxPixels=1e10) \\\n", + " .get('array'))\n", + " Sxx = Sigma.get([0, 0])\n", + " Syy = Sigma.get([1, 1])\n", + " Sxy = Sigma.get([0, 1])\n", + " # Correlation.\n", + " rho = Sxy.divide(Sxx.multiply(Syy).sqrt())\n", + " # Orthoregress reference onto target.\n", + " b = Syy.subtract(Sxx) \\\n", + " .add(Syy.subtract(Sxx).pow(2).add(Sxy.pow(2).multiply(4)).sqrt()) \\\n", + " .divide(Sxy.multiply(2))\n", + " a = Ym.subtract(b.multiply(Xm))\n", + " coeffs = coeffs.add(ee.List([b, a, rho]))\n", + " return ee.Dictionary({'image': image, 'coeffs': coeffs})\n", + "\n", + "def normalize(ref, target, imadID, pmin=0.9, bandNames=['B2', 'B3', 'B4', 'B8', 'B11', 'B12']):\n", + " '''\n", + " Relative radiometric normalization of target to reference\n", + " using iMAD for 6-band visual/infrared images\n", + "\n", + " Args:\n", + " - ref: Reference image\n", + " - target: Target image\n", + " - imadID: Asset ID of iMAD image\n", + " - pmin: Minimum p-value for identifying no change\n", + " - bandNames: Band names (default: S2 visual/infrared)\n", + "\n", + " Returns:\n", + " - Tuple containing normalized target, coefficients, image stack, and no-change mask\n", + " '''\n", + " try:\n", + " # Get iMAD result and mask out water surfaces.\n", + " waterMask = ee.Image('UMD/hansen/global_forest_change_2015').select('datamask').eq(1)\n", + " res = ee.Image(EXPORT_PATH + imadID).updateMask(waterMask)\n", + " # iMAD image.\n", + " imad = res.select(0, 1, 2, 3, 4, 5)\n", + " # chi-square test statistic.\n", + " z = res.select(6).rename('Z')\n", + " niter = imad.get('niter').getInfo()\n", + " rhos = ee.List(imad.get('rhos')).getInfo()\n", + " print('iMAD result:')\n", + " print('iterations: %i'%niter)\n", + " print('canonical correlations: %s'%rhos)\n", + " print('orthogonal regression:')\n", + " # p-values.\n", + " pval = chi2cdf(z, 6).subtract(1).multiply(-1).rename('pval')\n", + " # no-change mask using p-values greater than pmin.\n", + " noChangeMask = pval.gt(pmin)\n", + " # Stack (concatenate) the reference and target images.\n", + " im_stack = ee.Image.cat(ref, target).clip(aoi).updateMask(noChangeMask)\n", + " # iterate ortho regression over spectral bands.\n", + " first = ee.Dictionary({'image': im_stack, 'coeffs': ee.List([])})\n", + " result = ee.Dictionary(ee.List.sequence(0, 5).iterate(orthoregress, first))\n", + " # Display results.\n", + " coeffs = np.array(result.get('coeffs').getInfo())\n", + " print(' Slope Intercept Rho')\n", + " pprint(coeffs)\n", + " a = coeffs[:, 1].tolist()\n", + " b = coeffs[:, 0].tolist()\n", + " # Normalize target: Y = b X + a =\u003e X = (Y - a) / b\n", + " target_norm = target.subtract(ee.Image.constant(a)).divide(ee.Image.constant(b)).rename(bandNames)\n", + " return (target_norm, coeffs, im_stack, noChangeMask)\n", + " except ee.EEException as e:\n", + " print(\"Error:\", e)\n", + " return (None, None, None, None)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "H0ZfS3Z8LXYs" + }, + "source": [ + "### A Slightly Artificial Example\n", + "\n", + "As a first illustration, we collect two minimally cloudy Sentinel-2 images over the Landkreis Olpe aoi from June, 2022 using the *collect* function from [Part 2](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt2), capturing both SR and TOA reflectances:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "swEp63qrLXYs", + "tags": [] + }, + "outputs": [], + "source": [ + "# Landkreis Olpe.\n", + "aoi = ee.FeatureCollection('projects/google/imad_tutorial_aoi').geometry()\n", + "\n", + "# Spectral band combinations\n", + "#Sentinel-2\n", + "visirbands = ['B2', 'B3', 'B4', 'B8', 'B11', 'B12']\n", + "visnirbands = ['B8','B3','B2']\n", + "visbands = ['B2', 'B3', 'B4']\n", + "irbands = ['B8', 'B11', 'B12']\n", + "rededgebands = ['B5', 'B6', 'B7', 'B8A']\n", + "visnirparams = {'min': [0,0,0], 'max': [6000,1500,1500]}\n", + "visparams = {'min': [0,0,0], 'max': [1500,1500,1500]}\n", + "# Landsat-9 surface reflectance\n", + "visirbands_ls = ['SR_B2','SR_B3','SR_B4','SR_B5','SR_B6','SR_B7']\n", + "visnirbands_ls = ['SR_B5','SR_B3','SR_B2']\n", + "# Rescale LS-9 SR to Sentinel-2 SR\n", + "def rescale_ls(image):\n", + " return image.select(1,2,3,4,5,6).multiply(0.0000275).add(-0.2).multiply(10000)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "HkWNsPJmLXYs", + "tags": [] + }, + "outputs": [], + "source": [ + "im1_sr, im2_toa = collect(aoi, '2022-06-15', '2022-06-17', '2022-06-22', '2022-06-24', 'COPERNICUS/S2_SR_HARMONIZED', coll2 = 'COPERNICUS/S2_HARMONIZED')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "kvEXikFyLXYs" + }, + "source": [ + "Now we will do a relative radiometric normalization of the TOA image to the SR image. This is obviously a constructed example, as the target image, corrected to surface reflectance, is already available.\n", + "\n", + "First we require the iMAD result:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "yg9UDkAyLXYs", + "tags": [] + }, + "outputs": [], + "source": [ + "run_imad(aoi, im1_sr.select(visirbands), im2_toa.select(visirbands), 'MAD_Im1_sr_Im2_toa')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "iBZwDx1WLXYs" + }, + "source": [ + "After the iMAD export is finished (monitor progress at [https://code.earthengine.google.com/tasks](https://code.earthengine.google.com/tasks)), we can apply the normalization procedure, choosing a minimum p-value of 0.9 (the default) for the no change pixels:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "lTIMR6qHLXYs" + }, + "outputs": [], + "source": [ + "target = im2_toa.select(visirbands)\n", + "reference = im1_sr.select(visirbands)\n", + "im2_toa_norm, coeffs, im_stack, noChangeMask = normalize(reference, target, 'MAD_Im1_sr_Im2_toa')\n", + "# No-change image.\n", + "nc = ee.Image.constant(1).clip(aoi).updateMask(noChangeMask)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "2G_Gqhu1LXYs" + }, + "source": [ + "The $\\rho$ values are all $\u003e0.96$ indicating well-defined regressions. Here we have a look at them:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "cTBDTXMlLXYs", + "tags": [] + }, + "outputs": [], + "source": [ + "def plot_orthoregress(coeffs, im_stack, bandNames=visirbands, aoi=None):\n", + " '''Plot the regression fits for 6 spectral bands.'''\n", + " fig, ax = plt.subplots(2, 3, figsize=(12, 8))\n", + " a = coeffs[:, 1].tolist()\n", + " b = coeffs[:, 0].tolist()\n", + "\n", + " for i in range(3):\n", + " # Visual bands.\n", + " try:\n", + " d = np.array(ee.Array(im_stack.select(i, 6 + i)\n", + " .reduceRegion(ee.Reducer.toList(2), geometry=aoi, scale=10, maxPixels=1e10)\n", + " .get('list')).getInfo())\n", + " ax[0, i].plot(d[:, 0], d[:, 1], '.', label=bandNames[i])\n", + " ax[0, i].set_ylim(-50, 2000)\n", + " ax[0, i].set_xlim(-50, 2000)\n", + " ax[0, i].plot([-50, 2000], [a[i] - 50 * b[i], a[i] + 2000 * b[i]], '-r', label='ortho regression')\n", + " ax[0, i].legend()\n", + " ax[0, i].grid()\n", + " except ee.EEException as e:\n", + " print(\"Error fetching visual band data:\", e)\n", + "\n", + " # IR bands.\n", + " try:\n", + " d = np.array(ee.Array(im_stack.select(i + 3, 9 + i)\n", + " .reduceRegion(ee.Reducer.toList(2), geometry=aoi, scale=10, maxPixels=1e10)\n", + " .get('list')).getInfo())\n", + " ax[1, i].plot(d[:, 0], d[:, 1], '.', label=bandNames[i + 3])\n", + " ax[1, i].set_ylim(-50, 6000)\n", + " ax[1, i].set_xlim(-50, 6000)\n", + " ax[1, i].plot([-50, 6000], [a[i + 3] - 50 * b[i + 3], a[i + 3] + 6000 * b[i + 3]], '-r', label='ortho regression')\n", + " ax[1, i].legend()\n", + " ax[1, i].grid()\n", + " except ee.EEException as e:\n", + " print(\"Error fetching IR band data:\", e)\n", + "\n", + " plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "FhEzbe0iLXYs", + "tags": [] + }, + "outputs": [], + "source": [ + "plot_orthoregress( coeffs, im_stack )" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "TfEyqgJmLXYs" + }, + "source": [ + "The next cell compares the normalized and unnormalized TOA images with the reference SR image. The no-change pixels used in the regression are also displayed:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "EcOdPV1VLXYt" + }, + "outputs": [], + "source": [ + "M1 = geemap.Map()\n", + "M1.centerObject(aoi, 13)\n", + "M1.addLayer(im1_sr.select(visnirbands), visnirparams, 'Im1_sr')\n", + "M1.addLayer(im2_toa.select(visnirbands), visnirparams, 'Im2_toa')\n", + "M1.addLayer(im2_toa_norm.select(visnirbands), visnirparams, 'Im2_toa_norm')\n", + "M1.addLayer(nc,{'min': 0, 'max': 1, 'palette': ['yellow']}, 'no change')\n", + "\n", + "M1" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "nCNwpnwXLXYt" + }, + "source": [ + "As we might expect, comparison of the NDVI indices is meaningful after the radiometric normalization, but not before:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "t9eVcACJLXYt" + }, + "outputs": [], + "source": [ + "def ndvi_s2(im):\n", + " '''Sentinel-2 NDVI.'''\n", + " return im.select('B8').subtract(im.select('B4')).divide(\n", + " im.select('B8').add(im.select('B4')))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "GWjtm9pZLXYt" + }, + "outputs": [], + "source": [ + "M2 = geemap.Map()\n", + "M2.centerObject(aoi, 13)\n", + "M2.addLayer(ndvi_s2(im1_sr), {'min': 0, 'max': 1.2, 'palette': ['black', 'yellow']}, 'Im1_ndvi')\n", + "M2.addLayer(ndvi_s2(im2_toa_norm), {'min': 0, 'max': 1.2, 'palette': ['black', 'yellow']}, 'Im2_toa_norm_ndvi')\n", + "M2.addLayer(ndvi_s2(im2_toa), {'min': 0, 'max': 1.2, 'palette': ['black', 'yellow']}, 'Im2_toa_ndvi')\n", + "\n", + "M2" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "x3UezjdpLXYt" + }, + "source": [ + "### A More Realistic Example\n", + "\n", + "Top of atmosphere reflectance images are available for Sentinel-2 from 2015-06-23, surface reflectance images only from 2017-03-28. So we might attempt to normalize historical TOA images to their earliest SR counterparts. Here, for example, we gather an SR image from May, 2017 and a TOA image from May, 2016:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "AUGRc--6LXYt", + "tags": [] + }, + "outputs": [], + "source": [ + "im1_sr, im2_toa = collect(aoi, '2017-05-01', '2017-05-31', '2016-05-01', '2016-05-31', 'COPERNICUS/S2_SR_HARMONIZED', coll2='COPERNICUS/S2_HARMONIZED')" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "TUHFeCKCLXYt" + }, + "source": [ + "Again we first run the iMAD algorithm on the two images:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "RxCFf_wKLXYt", + "tags": [] + }, + "outputs": [], + "source": [ + "run_imad(aoi, im1_sr.select(visirbands), im2_toa.select(visirbands), 'MAD_Im1_sr_Im1_toa', maxiter=200)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "16Pqn574LXYt" + }, + "source": [ + "Once the export completes (check [https://code.earthengine.google.com/tasks](https://code.earthengine.google.com/tasks)), perform the orthogonal regression on the invariant pixels:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "na6m1MV1LXYt" + }, + "outputs": [], + "source": [ + "target = im2_toa.select(visirbands)\n", + "reference = im1_sr.select(visirbands)\n", + "im2_toa_norm, coeffs, im_stack, noChangeMask = normalize(reference, target, 'MAD_Im1_sr_Im1_toa')\n", + "# No-change image.\n", + "nc = ee.Image.constant(1).clip(aoi).updateMask(noChangeMask)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "02buhAX-LXYt", + "tags": [] + }, + "outputs": [], + "source": [ + "plot_orthoregress(coeffs, im_stack)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "0hx3uLejLXYt" + }, + "source": [ + "The above plot for the NIR band B8 is interesting as it seems to indicate that the surface reflectance correction is slightly nonlinear.\n", + "In any case, the regression plots look reasonable, so here again is the comparison of the normalized and unnormalized TOA images with the reference SR image:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "R0ZM9GLALXYt" + }, + "outputs": [], + "source": [ + "M3 = geemap.Map()\n", + "M3.centerObject(aoi, 13)\n", + "M3.addLayer(im1_sr.select(visnirbands), visnirparams, 'Im1_sr')\n", + "M3.addLayer(im2_toa.select(visnirbands), visnirparams, 'Im1_toa')\n", + "M3.addLayer(im2_toa_norm.select(visnirbands), visnirparams, 'Im2_toa_norm')\n", + "M3.addLayer(nc,{'min': 0, 'max': 1, 'palette': ['yellow']}, 'no change')\n", + "\n", + "M3" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "2j0fXpFPLXYu" + }, + "source": [ + "### Harmonization\n", + "\n", + "The term *harmonization*, with reference to Sentinel-2 data, has a special meaning, see [the note](https://developers.google.com/earth-engine/datasets/catalog/sentinel-2) in the GEE data catalogue. However, here we'll use it to characterize relative radiometric normalization across two different remote sensing satellite missions.\n", + "\n", + "Landsat-9 surface reflectance images are available from October, 2021, Landsat-8 SR from March, 2013. Paralleling our investigation of clear cutting in [Part 2](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt2), we gather a cloud-free Landsat-8 image from March 30, 2021 and a Landsat-9 acquisition from March 9, 2022, likewise cloud-free. For convenience, the reflectances are rescaled to the same range $[0,10000]$ as for Sentinel-2." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "saapPrVDLXYu", + "tags": [] + }, + "outputs": [], + "source": [ + "im_ls8, im_ls9 = collect(aoi, '2021-03-29', '2021-03-31', '2022-03-08', '2022-03-10', 'LANDSAT/LC08/C02/T1_L2', coll2='LANDSAT/LC09/C02/T1_L2')\n", + "im_ls8 = rescale_ls(im_ls8)\n", + "im_ls9 = rescale_ls(im_ls9)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "779R1GokLXYu" + }, + "source": [ + "Run the iMAD algorithm:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "uU_VNwdALXYu", + "tags": [] + }, + "outputs": [], + "source": [ + "run_imad(aoi, im_ls8.select(visirbands_ls), im_ls9.select(visirbands_ls), 'MAD_Im_ls8_Im_ls9', scale=30, maxiter=200)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "kdy6BSlBLXYu" + }, + "source": [ + "Once the export completes (check [https://code.earthengine.google.com/tasks](https://code.earthengine.google.com/tasks)), perform the normalization (or rather, the harmonization):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "S4754xXXLXYu" + }, + "outputs": [], + "source": [ + "target = im_ls8.select(visirbands_ls)\n", + "reference = im_ls9.select(visirbands_ls)\n", + "im_ls8_norm, coeffs, im_stack, noChangeMask = normalize(reference, target, 'MAD_Im_ls8_Im_ls9', pmin=0.9, bandNames=visirbands_ls)\n", + "# No-change image.\n", + "nc = ee.Image.constant(1).clip(aoi).updateMask(noChangeMask)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "QKLYqORfLXYu", + "tags": [] + }, + "outputs": [], + "source": [ + "plot_orthoregress(coeffs, im_stack, bandNames=visirbands_ls)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "mj3J_xE5LXYu" + }, + "source": [ + "For most of the bands, the intercepts are small and the slopes fairly close to one. The NIR band 'SR_B5' is the exception. This is apparent in the direct comparison of the original and harmonized Landsat 8 images:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "hwJzmws3LXYv" + }, + "outputs": [], + "source": [ + "M4 = geemap.Map()\n", + "M4.centerObject(aoi, 11)\n", + "M4.addLayer(im_ls9.select(visnirbands_ls), visnirparams, 'Im_ls9')\n", + "M4.addLayer(im_ls8.select(visnirbands_ls), visnirparams, 'Im_ls8')\n", + "M4.addLayer(im_ls8_norm.select(visnirbands_ls), visnirparams, 'Im_ls8_norm')\n", + "M4.addLayer(nc,{'min': 0, 'max': 1, 'palette': ['yellow']}, 'no change')\n", + "\n", + "M4" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "Rhr8B90_LXYv" + }, + "source": [ + "### Conclusion\n", + "\n", + "This ends our three-part tutorial on GEE change detection with the iMAD transformation. We hope it has been useful and informative for anyone facing the task of discriminating changes and/or invariances in multispectral image pairs, particularly in the GEE environment." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "6kaumJrnLXYv" + }, + "source": [ + "### Exercises" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "Z3_PrBpWLXYv" + }, + "source": [ + "1. By clustering the Landsat iMAD image from the last example, try to determine the total forest cover loss due to clear cutting between the two acquisitions, see [Part 2](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt2)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "TSOSFoAjLXYv" + }, + "outputs": [], + "source": [ + "im_imad = ee.Image(EXPORT_PATH + 'MAD_Im_ls8_Im_ls9').select(0, 1, 2, 3, 4, 5)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "UWWKXxHBLXYv" + }, + "source": [ + "2. The following code generates a time series of 8 relatively cloud-free Sentinel-2 TOA images from June to September, 2022. Use the methods discussed above to perform a relative radiometric normalization of the entire series." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "UyH71PxCLXYv" + }, + "outputs": [], + "source": [ + "im1, im2 = collect(aoi, '2022-06-15', '2022-06-17', '2022-06-22', '2022-06-24', 'COPERNICUS/S2_HARMONIZED')\n", + "im3, im4 = collect(aoi, '2022-06-25', '2022-07-09', '2022-07-09', '2022-07-24', 'COPERNICUS/S2_HARMONIZED')\n", + "im5, im6 = collect(aoi, '2022-07-28', '2022-07-31', '2022-08-01', '2022-08-31', 'COPERNICUS/S2_HARMONIZED')\n", + "im7, im8 = collect(aoi, '2022-08-13', '2022-08-31', '2022-09-01', '2022-09-30', 'COPERNICUS/S2_HARMONIZED')" + ] + } + ], + "metadata": { + "colab": { + "provenance": [] + }, + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/tutorials/imad-tutorial-pt3/index.md b/tutorials/imad-tutorial-pt3/index.md deleted file mode 100644 index 9e8fd22f2..000000000 --- a/tutorials/imad-tutorial-pt3/index.md +++ /dev/null @@ -1,660 +0,0 @@ ---- -title: Change Detection in GEE - Part 3 Relative Radiometric Normalization and Harmonization -description: Iteratively re-weighted Multivariate Alteration Detection. -author: mortcanty -tags: change-detection, imad -date_published: 2023-03-24 ---- - - -{% include "earth-engine/_widgets/_notebook_buttons.html" with github_path="mortcanty/eeimad/blob/main/src/gee_change_detection_pt_3.ipynb" %} - -Test of pixel intensities across different image acquisitions requires that the received signals have similar radiometric scales. -As we saw in [Part 2](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt2) of the present tutorial, provided the iMAD algorithm converges satisfactorily, it will isolate the pixels -which have not significantly changed in the two input images presented to it. A regression analysis on these no-change -pixels can then determine how well the radiometric parameters of the two acquisitions compare. -If, for instance, images that were corrected to absolute surface reflectance are examined in this way, we would expect a good match, i.e., a -slope close to one and an intercept close to zero at all spectral wavelengths. In general, -for uncalibrated or poorly corrected images, this will not be the case. Then, the regression -coefficients can be applied to normalize one image (the target, say) to the other -(the reference). This is a prerequisite, for example, for tracing features such as NDVI indices over -a time series when the images have not been reduced to surface reflectances, -see e.g., [Gan et al. (2021)](https://ieeexplore.ieee.org/document/9392311), or indeed for -*harmonizing* the data from two different sensors of the same family, such as Landsat 8 with -Landsat 9. In this final part of the GEE Change Detection Tutorial we'll have a closer look at this approach, which we'll refer to here as *relative radiometric normalization*. For more detailed treatment, see also [Canty and Nielsen (2008)](https://www.sciencedirect.com/science/article/abs/pii/S0034425707003495), [Canty et al. (2004)](https://www.sciencedirect.com/science/article/abs/pii/S0034425704001208), and [Schroeder et al. (2006)](https://www.sciencedirect.com/science/article/abs/pii/S0034425706001179). A different (but similar) method involving identification of so-called *pseudoinvariant features* -is explained in [another GEE community tutorial](https://developers.google.com/earth-engine/tutorials/community/pseudo-invariant-feature-matching) in this series. - -## Preliminaries - -```python -# Enter your own export to assets path name here ----------- -EXPORT_PATH = 'projects/YOUR_GEE_PROJECT_NAME/assets/imad/' -# ------------------------------------------------ -``` - -```python -import ee - -ee.Authenticate(auth_mode='notebook') -ee.Initialize() -``` - -```python -# Import other packages used in the tutorial -%matplotlib inline -import geemap -import numpy as np -import random, time -import matplotlib.pyplot as plt -from scipy.stats import norm, chi2 - -from pprint import pprint # for pretty printing -``` - -```python -#@title Routines from Parts 1 and 2 - -def trunc(values, dec = 3): - '''Truncate a 1-D array to dec decimal places.''' - return np.trunc(values*10**dec)/(10**dec) - -# Display an image in a one percent linear stretch. -def display_ls(image, map, name, centered = False): - bns = image.bandNames().length().getInfo() - if bns == 3: - image = image.rename('B1', 'B2', 'B3') - pb_99 = ['B1_p99', 'B2_p99', 'B3_p99'] - pb_1 = ['B1_p1', 'B2_p1', 'B3_p1'] - img = ee.Image.rgb( - image.select('B1'), image.select('B2'), image.select('B3')) - else: - image = image.rename('B1') - pb_99 = ['B1_p99'] - pb_1 = ['B1_p1'] - img = image.select('B1') - percentiles = image.reduceRegion( - ee.Reducer.percentile([1, 99]), maxPixels=1e11) - mx = percentiles.values(pb_99) - if centered: - mn = ee.Array(mx).multiply(-1).toList() - else: - mn = percentiles.values(pb_1) - map.addLayer(img, {'min': mn, 'max': mx}, name) - -def covarw(image, weights=None, scale=20, maxPixels=1e10): - '''Return the centered image and its weighted covariance matrix.''' - try: - geometry = image.geometry() - bandNames = image.bandNames() - N = bandNames.length() - if weights is None: - weights = image.constant(1) - weightsImage = image.multiply(ee.Image.constant(0)).add(weights) - means = image.addBands(weightsImage) \ - .reduceRegion(ee.Reducer.mean().repeat(N).splitWeights(), - scale=scale, - maxPixels=maxPixels) \ - .toArray() \ - .project([1]) - centered = image.toArray().subtract(means) - B1 = centered.bandNames().get(0) - b1 = weights.bandNames().get(0) - nPixels = ee.Number(centered.reduceRegion(ee.Reducer.count(), - scale=scale, - maxPixels=maxPixels).get(B1)) - sumWeights = ee.Number(weights.reduceRegion(ee.Reducer.sum(), - geometry=geometry, - scale=scale, - maxPixels=maxPixels).get(b1)) - covw = centered.multiply(weights.sqrt()) \ - .toArray() \ - .reduceRegion(ee.Reducer.centeredCovariance(), - geometry=geometry, - scale=scale, - maxPixels=maxPixels) \ - .get('array') - covw = ee.Array(covw).multiply(nPixels).divide(sumWeights) - return (centered.arrayFlatten([bandNames]), covw) - except Exception as e: - print('Error: %s'%e) - -def corr(cov): - '''Transform covariance matrix to correlation matrix.''' - # Diagonal matrix of inverse sigmas. - sInv = cov.matrixDiagonal().sqrt().matrixToDiag().matrixInverse() - # Transform. - corr = sInv.matrixMultiply(cov).matrixMultiply(sInv).getInfo() - # Truncate. - return [list(map(trunc, corr[i])) for i in range(len(corr))] - -def geneiv(C,B): - '''Return the eigenvalues and eigenvectors of the generalized eigenproblem - C*X = lambda*B*X''' - try: - C = ee.Array(C) - B = ee.Array(B) - # Li = choldc(B)^-1 - Li = ee.Array(B.matrixCholeskyDecomposition().get('L')).matrixInverse() - # Solve symmetric, ordinary eigenproblem Li*C*Li^T*x = lambda*x - Xa = Li.matrixMultiply(C) \ - .matrixMultiply(Li.matrixTranspose()) \ - .eigen() - # Eigenvalues as a row vector. - lambdas = Xa.slice(1, 0, 1).matrixTranspose() - # Eigenvectors as columns. - X = Xa.slice(1, 1).matrixTranspose() - # Generalized eigenvectors as columns, Li^T*X - eigenvecs = Li.matrixTranspose().matrixMultiply(X) - return (lambdas, eigenvecs) - except Exception as e: - print('Error: %s'%e) - -# Collect a Sentinel-2 image pair. -def collect(aoi, t1a ,t1b, t2a, t2b, coll1, coll2 = None): - try: - if coll2 == None: - coll2 = coll1 - im1 = ee.Image( ee.ImageCollection(coll1) - .filterBounds(aoi) - .filterDate(ee.Date(t1a), ee.Date(t1b)) - .filter(ee.Filter.contains(rightValue=aoi, leftField='.geo')) - .sort('CLOUDY_PIXEL_PERCENTAGE') - .sort('CLOUD_COVER_LAND') - .first() - .clip(aoi) ) - im2 = ee.Image( ee.ImageCollection(coll2) - .filterBounds(aoi) - .filterDate(ee.Date(t2a), ee.Date(t2b)) - .filter(ee.Filter.contains(rightValue=aoi, leftField='.geo')) - .sort('CLOUDY_PIXEL_PERCENTAGE') - .sort('CLOUD_COVER_LAND') - .first() - .clip(aoi) ) - timestamp = im1.date().format('E MMM dd HH:mm:ss YYYY') - print(timestamp.getInfo()) - timestamp = im2.date().format('E MMM dd HH:mm:ss YYYY') - print(timestamp.getInfo()) - return (im1, im2) - except Exception as e: - print('Error: %s'%e) - -#@title The iMAD code -def chi2cdf(Z,df): - '''Chi-square cumulative distribution function with df degrees of freedom.''' - return ee.Image(Z.divide(2)).gammainc(ee.Number(df).divide(2)) - -def imad(current,prev): - '''Iterator function for iMAD.''' - done = ee.Number(ee.Dictionary(prev).get('done')) - return ee.Algorithms.If(done, prev, imad1(current, prev)) - -def imad1(current,prev): - '''Iteratively re-weighted MAD.''' - image = ee.Image(ee.Dictionary(prev).get('image')) - Z = ee.Image(ee.Dictionary(prev).get('Z')) - allrhos = ee.List(ee.Dictionary(prev).get('allrhos')) - nBands = image.bandNames().length().divide(2) - weights = chi2cdf(Z,nBands).subtract(1).multiply(-1) - scale = ee.Dictionary(prev).getNumber('scale') - niter = ee.Dictionary(prev).getNumber('niter') - # Weighted stacked image and weighted covariance matrix. - centeredImage, covarArray = covarw(image, weights, scale) - bNames = centeredImage.bandNames() - bNames1 = bNames.slice(0, nBands) - bNames2 = bNames.slice(nBands) - centeredImage1 = centeredImage.select(bNames1) - centeredImage2 = centeredImage.select(bNames2) - s11 = covarArray.slice(0, 0, nBands).slice(1, 0, nBands) - s22 = covarArray.slice(0, nBands).slice(1, nBands) - s12 = covarArray.slice(0, 0, nBands).slice(1, nBands) - s21 = covarArray.slice(0, nBands).slice(1, 0, nBands) - c1 = s12.matrixMultiply(s22.matrixInverse()).matrixMultiply(s21) - b1 = s11 - c2 = s21.matrixMultiply(s11.matrixInverse()).matrixMultiply(s12) - b2 = s22 - # Solution of generalized eigenproblems. - lambdas, A = geneiv(c1, b1) - _, B = geneiv(c2, b2) - rhos = lambdas.sqrt().project(ee.List([1])) - # Test for convergence. - lastrhos = ee.Array(allrhos.get(-1)) - done = rhos.subtract(lastrhos) \ - .abs() \ - .reduce(ee.Reducer.max(), ee.List([0])) \ - .lt(ee.Number(0.0001)) \ - .toList() \ - .get(0) - allrhos = allrhos.cat([rhos.toList()]) - # MAD variances. - sigma2s = rhos.subtract(1).multiply(-2).toList() - sigma2s = ee.Image.constant(sigma2s) - # Ensure sum of positive correlations between X and U is positive. - tmp = s11.matrixDiagonal().sqrt() - ones = tmp.multiply(0).add(1) - tmp = ones.divide(tmp).matrixToDiag() - s = tmp.matrixMultiply(s11).matrixMultiply(A).reduce(ee.Reducer.sum(), [0]).transpose() - A = A.matrixMultiply(s.divide(s.abs()).matrixToDiag()) - # Ensure positive correlation. - tmp = A.transpose().matrixMultiply(s12).matrixMultiply(B).matrixDiagonal() - tmp = tmp.divide(tmp.abs()).matrixToDiag() - B = B.matrixMultiply(tmp) - # Canonical and MAD variates. - centeredImage1Array = centeredImage1.toArray().toArray(1) - centeredImage2Array = centeredImage2.toArray().toArray(1) - U = ee.Image(A.transpose()).matrixMultiply(centeredImage1Array) \ - .arrayProject([0]) \ - .arrayFlatten([bNames1]) - V = ee.Image(B.transpose()).matrixMultiply(centeredImage2Array) \ - .arrayProject([0]) \ - .arrayFlatten([bNames2]) - iMAD = U.subtract(V) - # Chi-square image. - Z = iMAD.pow(2) \ - .divide(sigma2s) \ - .reduce(ee.Reducer.sum()) - return ee.Dictionary({'done': done, 'scale': scale, 'niter': niter.add(1), - 'image': image, 'allrhos': allrhos, 'Z': Z, 'iMAD': iMAD}) - -#@title Run iMAD algorithm as export task -def run_imad(aoi, image1, image2, assetFN, scale=10, maxiter=100): - try: - N = image1.bandNames().length().getInfo() - imadnames = ['iMAD'+str(i+1) for i in range(N)] - imadnames.append('Z') - # Maximum iterations. - inputlist = ee.List.sequence(1, maxiter) - first = ee.Dictionary({'done':0, - 'scale': scale, - 'niter': ee.Number(0), - 'image': image1.addBands(image2), - 'allrhos': [ee.List.sequence(1, N)], - 'Z': ee.Image.constant(0), - 'iMAD': ee.Image.constant(0)}) - # Iteration. - result = ee.Dictionary(inputlist.iterate(imad, first)) - # Retrieve results. - iMAD = ee.Image(result.get('iMAD')).clip(aoi) - rhos = ee.String.encodeJSON(ee.List(result.get('allrhos')).get(-1)) - Z = ee.Image(result.get('Z')) - niter = result.getNumber('niter') - # Export iMAD and Z as a single image, including rhos and number of iterations in properties. - iMAD_export = ee.Image.cat(iMAD, Z).rename(imadnames).set('rhos', rhos, 'niter', niter) - assetId = EXPORT_PATH + assetFN - assexport = ee.batch.Export.image.toAsset(iMAD_export, - description='assetExportTask', - assetId=assetId, scale=scale, maxPixels=1e10) - assexport.start() - print('Exporting iMAD to %s\n task id: %s'%(assetId, str(assexport.id))) - except Exception as e: - print('Error: %s'%e) -``` - -### Relative radiometric normalization - -We will illustrate the idea with the German administrative district scene used in the first two parts of this tutorial ([Part 1](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt1), [Part 2](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt2)), the Landkreis Olpe. Once we have isolated what we think to be the no-change pixels for two different Sentinel-2 acquisitions in the aoi, we can perform a regression analysis on them to determine how well their radiometric parameters compare *relative* to the identified no-change pixels in the two scenes. Clearly, both variables involved have measurement uncertainty associated with them. In fact, which acquisition is termed reference and which is termed target data is arbitrary. Therefore it is preferable to use *orthogonal linear regression* rather than ordinary linear regression. The former method treats the data symmetrically whereas the latter assumes uncertainty only in the dependent variable. - -Orthogonal linear regression on two variables is explained in detail in the Appendix of [Canty et al. (2004)](https://www.sciencedirect.com/science/article/abs/pii/S0034425704001208), but can be briefly summarized as follows: let the no-change observations in a given spectral band of the reference and target images be $X_i$ and $Y_i, \ i=1\dots n$, respectively. Then define the $n\times 2$ *data matrix* - -$$ -d = \begin{pmatrix} X_1 & Y_1\cr X_2 & Y_2 \cr \vdots & \vdots \cr X_n & Y_n\end{pmatrix}, -$$ - -and calculate from it the mean values $\bar X,\ \bar Y$ and the covariance matrix - -$$ -\Sigma =\begin{pmatrix} \sigma^2_X & \sigma_{XY} \cr \sigma_{YX}\ & \sigma^2_Y\end{pmatrix}. -$$ - -The best estimates for the slope $b$ and intercept $a$ of the orthogonal regression line - -$$ -Y = bX+a -$$ - -are then given by - -$$ -\hat b = {\sigma^2_X -\sigma^2_Y + \sqrt{(\sigma^2_X -\sigma^2_Y)^2+4\sigma^2_{XY}}\over 2\sigma_{XY}}, \quad \hat a = \bar Y - b\bar X. -$$ - -As a simple quality criterion we can use the correlation coefficient - -$$\rho = {\sigma_{XY}\over \sigma_X\sigma_Y}$$ - -which has value one for a perfect fit. - -The following cell codes a function to iterate orthogonal regression over all of the spectral bands of an image pair. - -```python -def orthoregress(current, prev): - ''' - Iterator function for orthogonal regression - ''' - k = ee.Number(current) - prev = ee.Dictionary(prev) - # Image is concatenation of reference and target. - image = ee.Image(prev.get('image')) - coeffs = ee.List(prev.get('coeffs')) - N = image.bandNames().length().divide(2) - # Data matrix. - d = image.select(k, k.add(N)) - means = d.reduceRegion(ee.Reducer.mean(), scale=10, maxPixels=1e10) \ - .toArray() \ - .project([0]) - Xm = means.get([0]) - Ym = means.get([1]) - Sigma = ee.Array(d.toArray() \ - .reduceRegion(ee.Reducer.covariance(), scale=10, maxPixels=1e10) \ - .get('array')) - Sxx = Sigma.get([0, 0]) - Syy = Sigma.get([1, 1]) - Sxy = Sigma.get([0, 1]) - # Correlation. - rho = Sxy.divide(Sxx.multiply(Syy).sqrt()) - # Orthoregress reference onto target. - b = Syy.subtract(Sxx) \ - .add(Syy.subtract(Sxx).pow(2).add(Sxy.pow(2).multiply(4)).sqrt()) \ - .divide(Sxy.multiply(2)) - a = Ym.subtract(b.multiply(Xm)) - coeffs = coeffs.add(ee.List([b, a, rho])) - return ee.Dictionary({'image': image, 'coeffs': coeffs}) - -def normalize(ref, target, imadID, pmin=0.9, bandNames=['B2', 'B3', 'B4', 'B8', 'B11', 'B12']): - ''' - Relative radiometric normalization of target to reference - using iMAD for 6-band visual/infrared images - - Args: - - ref: Reference image - - target: Target image - - imadID: Asset ID of iMAD image - - pmin: Minimum p-value for identifying no change - - bandNames: Band names (default: S2 visual/infrared) - - Returns: - - Tuple containing normalized target, coefficients, image stack, and no-change mask - ''' - try: - # Get iMAD result and mask out water surfaces. - waterMask = ee.Image('UMD/hansen/global_forest_change_2015').select('datamask').eq(1) - res = ee.Image(EXPORT_PATH + imadID).updateMask(waterMask) - # iMAD image. - imad = res.select(0, 1, 2, 3, 4, 5) - # chi-square test statistic. - z = res.select(6).rename('Z') - niter = imad.get('niter').getInfo() - rhos = ee.List(imad.get('rhos')).getInfo() - print('iMAD result:') - print('iterations: %i'%niter) - print('canonical correlations: %s'%rhos) - print('orthogonal regression:') - # p-values. - pval = chi2cdf(z, 6).subtract(1).multiply(-1).rename('pval') - # no-change mask using p-values greater than pmin. - noChangeMask = pval.gt(pmin) - # Stack (concatenate) the reference and target images. - im_stack = ee.Image.cat(ref, target).clip(aoi).updateMask(noChangeMask) - # iterate ortho regression over spectral bands. - first = ee.Dictionary({'image': im_stack, 'coeffs': ee.List([])}) - result = ee.Dictionary(ee.List.sequence(0, 5).iterate(orthoregress, first)) - # Display results. - coeffs = np.array(result.get('coeffs').getInfo()) - print(' Slope Intercept Rho') - pprint(coeffs) - a = coeffs[:, 1].tolist() - b = coeffs[:, 0].tolist() - # Normalize target: Y = b X + a => X = (Y - a) / b - target_norm = target.subtract(ee.Image.constant(a)).divide(ee.Image.constant(b)).rename(bandNames) - return (target_norm, coeffs, im_stack, noChangeMask) - except ee.EEException as e: - print("Error:", e) - return (None, None, None, None) -``` - -### A Slightly Artificial Example - -As a first illustration, we collect two minimally cloudy Sentinel-2 images over the Landkreis Olpe aoi from June, 2022 using the *collect* function from [Part 2](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt2), capturing both SR and TOA reflectances: - -```python -# Landkreis Olpe area of interest. -aoi = ee.FeatureCollection('projects/google/imad_tutorial_aoi').geometry() - -# Spectral band combinations -#Sentinel-2 -visirbands = ['B2', 'B3', 'B4', 'B8', 'B11', 'B12'] -visnirbands = ['B8','B3','B2'] -visbands = ['B2', 'B3', 'B4'] -irbands = ['B8', 'B11', 'B12'] -rededgebands = ['B5', 'B6', 'B7', 'B8A'] -visnirparams = {'min': [0,0,0], 'max': [6000,1500,1500]} -visparams = {'min': [0,0,0], 'max': [1500,1500,1500]} -# Landsat-9 surface reflectance -visirbands_ls = ['SR_B2','SR_B3','SR_B4','SR_B5','SR_B6','SR_B7'] -visnirbands_ls = ['SR_B5','SR_B3','SR_B2'] -# Rescale LS-9 SR to Sentinel-2 SR -def rescale_ls(image): - return image.select(1,2,3,4,5,6).multiply(0.0000275).add(-0.2).multiply(10000) -``` - -```python -im1_sr, im2_toa = collect(aoi, '2022-06-15', '2022-06-17', '2022-06-22', '2022-06-24', 'COPERNICUS/S2_SR_HARMONIZED', coll2 = 'COPERNICUS/S2_HARMONIZED') -``` - -Now we will do a relative radiometric normalization of the TOA image to the SR image. This is obviously a constructed example, as the target image, corrected to surface reflectance, is already available. - -First we require the iMAD result: - -```python -run_imad(aoi, im1_sr.select(visirbands), im2_toa.select(visirbands), 'MAD_Im1_sr_Im2_toa') -``` - -After the iMAD export is finished (monitor progress at [https://code.earthengine.google.com/tasks](https://code.earthengine.google.com/tasks)), we can apply the normalization procedure, choosing a minimum p-value of 0.9 (the default) for the no change pixels: - -```python -target = im2_toa.select(visirbands) -reference = im1_sr.select(visirbands) -im2_toa_norm, coeffs, im_stack, noChangeMask = normalize(reference, target, 'MAD_Im1_sr_Im2_toa') -# No-change image. -nc = ee.Image.constant(1).clip(aoi).updateMask(noChangeMask) -``` - -The $\rho$ values are all $>0.96$ indicating well-defined regressions. Here we have a look at them: - -```python -def plot_orthoregress(coeffs, im_stack, bandNames=visirbands, aoi=None): - '''Plot the regression fits for 6 spectral bands.''' - fig, ax = plt.subplots(2, 3, figsize=(12, 8)) - a = coeffs[:, 1].tolist() - b = coeffs[:, 0].tolist() - - for i in range(3): - # Visual bands. - try: - d = np.array(ee.Array(im_stack.select(i, 6 + i) - .reduceRegion(ee.Reducer.toList(2), geometry=aoi, scale=10, maxPixels=1e10) - .get('list')).getInfo()) - ax[0, i].plot(d[:, 0], d[:, 1], '.', label=bandNames[i]) - ax[0, i].set_ylim(-50, 2000) - ax[0, i].set_xlim(-50, 2000) - ax[0, i].plot([-50, 2000], [a[i] - 50 * b[i], a[i] + 2000 * b[i]], '-r', label='ortho regression') - ax[0, i].legend() - ax[0, i].grid() - except ee.EEException as e: - print("Error fetching visual band data:", e) - - # IR bands. - try: - d = np.array(ee.Array(im_stack.select(i + 3, 9 + i) - .reduceRegion(ee.Reducer.toList(2), geometry=aoi, scale=10, maxPixels=1e10) - .get('list')).getInfo()) - ax[1, i].plot(d[:, 0], d[:, 1], '.', label=bandNames[i + 3]) - ax[1, i].set_ylim(-50, 6000) - ax[1, i].set_xlim(-50, 6000) - ax[1, i].plot([-50, 6000], [a[i + 3] - 50 * b[i + 3], a[i + 3] + 6000 * b[i + 3]], '-r', label='ortho regression') - ax[1, i].legend() - ax[1, i].grid() - except ee.EEException as e: - print("Error fetching IR band data:", e) - - plt.show() -``` - - -```python -plot_orthoregress( coeffs, im_stack ) -``` - -The next cell compares the normalized and unnormalized TOA images with the reference SR image. The no-change pixels used in the regression are also displayed: - -```python -M1 = geemap.Map() -M1.centerObject(aoi, 13) -M1.addLayer(im1_sr.select(visnirbands), visnirparams, 'Im1_sr') -M1.addLayer(im2_toa.select(visnirbands), visnirparams, 'Im2_toa') -M1.addLayer(im2_toa_norm.select(visnirbands), visnirparams, 'Im2_toa_norm') -M1.addLayer(nc,{'min': 0, 'max': 1, 'palette': ['yellow']}, 'no change') - -M1 -``` - -As we might expect, comparison of the NDVI indices is meaningful after the radiometric normalization, but not before: - -```python -def ndvi_s2(im): - '''Sentinel-2 NDVI.''' - return im.select('B8').subtract(im.select('B4')).divide( - im.select('B8').add(im.select('B4'))) -``` - -```python -M2 = geemap.Map() -M2.centerObject(aoi, 13) -M2.addLayer(ndvi_s2(im1_sr), {'min': 0, 'max': 1.2, 'palette': ['black', 'yellow']}, 'Im1_ndvi') -M2.addLayer(ndvi_s2(im2_toa_norm), {'min': 0, 'max': 1.2, 'palette': ['black', 'yellow']}, 'Im2_toa_norm_ndvi') -M2.addLayer(ndvi_s2(im2_toa), {'min': 0, 'max': 1.2, 'palette': ['black', 'yellow']}, 'Im2_toa_ndvi') - -M2 -``` - -### A More Realistic Example - -Top of atmosphere reflectance images are available for Sentinel-2 from 2015-06-23, surface reflectance images only from 2017-03-28. So we might attempt to normalize historical TOA images to their earliest SR counterparts. Here, for example, we gather an SR image from May, 2017 and a TOA image from May, 2016: - - -```python -im1_sr, im2_toa = collect(aoi, '2017-05-01', '2017-05-31', '2016-05-01', '2016-05-31', 'COPERNICUS/S2_SR_HARMONIZED', coll2='COPERNICUS/S2_HARMONIZED') -``` - -Again we first run the iMAD algorithm on the two images: - -```python -run_imad(aoi, im1_sr.select(visirbands), im2_toa.select(visirbands), 'MAD_Im1_sr_Im1_toa', maxiter=200) -``` - -Once the export completes (check [https://code.earthengine.google.com/tasks](https://code.earthengine.google.com/tasks)), perform the orthogonal regression on the invariant pixels: - -```python -target = im2_toa.select(visirbands) -reference = im1_sr.select(visirbands) -im2_toa_norm, coeffs, im_stack, noChangeMask = normalize(reference, target, 'MAD_Im1_sr_Im1_toa') -# No-change image. -nc = ee.Image.constant(1).clip(aoi).updateMask(noChangeMask) -``` - -```python -plot_orthoregress(coeffs, im_stack) -``` - -The above plot for the NIR band B8 is interesting as it seems to indicate that the surface reflectance correction is slightly nonlinear. -In any case, the regression plots look reasonable, so here again is the comparison of the normalized and unnormalized TOA images with the reference SR image: - -```python -M3 = geemap.Map() -M3.centerObject(aoi, 13) -M3.addLayer(im1_sr.select(visnirbands), visnirparams, 'Im1_sr') -M3.addLayer(im2_toa.select(visnirbands), visnirparams, 'Im1_toa') -M3.addLayer(im2_toa_norm.select(visnirbands), visnirparams, 'Im2_toa_norm') -M3.addLayer(nc,{'min': 0, 'max': 1, 'palette': ['yellow']}, 'no change') - -M3 -``` - -### Harmonization - -The term *harmonization*, with reference to Sentinel-2 data, has a special meaning, see [the note](https://developers.google.com/earth-engine/datasets/catalog/sentinel-2) in the GEE data catalogue. However, here we'll use it to characterize relative radiometric normalization across two different remote sensing satellite missions. - -Landsat-9 surface reflectance images are available from October, 2021, Landsat-8 SR from March, 2013. Paralleling our investigation of clear cutting in [Part 2](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt2), we gather a cloud-free Landsat-8 image from March 30, 2021 and a Landsat-9 acquisition from March 9, 2022, likewise cloud-free. For convenience, the reflectances are rescaled to the same range $[0,10000]$ as for Sentinel-2. - -```python -im_ls8, im_ls9 = collect(aoi, '2021-03-29', '2021-03-31', '2022-03-08', '2022-03-10', 'LANDSAT/LC08/C02/T1_L2', coll2='LANDSAT/LC09/C02/T1_L2') -im_ls8 = rescale_ls(im_ls8) -im_ls9 = rescale_ls(im_ls9) -``` - -Run the iMAD algorithm: - -```python -run_imad(aoi, im_ls8.select(visirbands_ls), im_ls9.select(visirbands_ls), 'MAD_Im_ls8_Im_ls9', scale=30, maxiter=200) -``` - -Once the export completes (check [https://code.earthengine.google.com/tasks](https://code.earthengine.google.com/tasks)), perform the normalization (or rather, the harmonization): - -```python -target = im_ls8.select(visirbands_ls) -reference = im_ls9.select(visirbands_ls) -im_ls8_norm, coeffs, im_stack, noChangeMask = normalize(reference, target, 'MAD_Im_ls8_Im_ls9', pmin=0.9, bandNames=visirbands_ls) -# No-change image. -nc = ee.Image.constant(1).clip(aoi).updateMask(noChangeMask) -``` - -```python -plot_orthoregress(coeffs, im_stack, bandNames=visirbands_ls) -``` - -For most of the bands, the intercepts are small and the slopes fairly close to one. The NIR band 'SR_B5' is the exception. This is apparent in the direct comparison of the original and harmonized Landsat 8 images: - -```python -M4 = geemap.Map() -M4.centerObject(aoi, 11) -M4.addLayer(im_ls9.select(visnirbands_ls), visnirparams, 'Im_ls9') -M4.addLayer(im_ls8.select(visnirbands_ls), visnirparams, 'Im_ls8') -M4.addLayer(im_ls8_norm.select(visnirbands_ls), visnirparams, 'Im_ls8_norm') -M4.addLayer(nc,{'min': 0, 'max': 1, 'palette': ['yellow']}, 'no change') - -M4 -``` - -### Conclusion - -This ends our three-part tutorial on GEE change detection with the iMAD transformation. We hope it has been useful and informative for anyone facing the task of discriminating changes and/or invariances in multispectral image pairs, particularly in the GEE environment. - -### Exercises - -1\. By clustering the Landsat iMAD image from the last example, try to determine the total forest cover loss due to clear cutting between the two acquisitions, see [Part 2](https://developers.google.com/earth-engine/tutorials/community/imad-tutorial-pt2). - -```python -im_imad = ee.Image(EXPORT_PATH + 'MAD_Im_ls8_Im_ls9').select(0, 1, 2, 3, 4, 5) -``` - -2\. The following code generates a time series of 8 relatively cloud-free Sentinel-2 TOA images from June to September, 2022. Use the methods discussed above to perform a relative radiometric normalization of the entire series. - -```python -im1, im2 = collect(aoi, '2022-06-15', '2022-06-17', '2022-06-22', '2022-06-24', 'COPERNICUS/S2_HARMONIZED') -im3, im4 = collect(aoi, '2022-06-25', '2022-07-09', '2022-07-09', '2022-07-24', 'COPERNICUS/S2_HARMONIZED') -im5, im6 = collect(aoi, '2022-07-28', '2022-07-31', '2022-08-01', '2022-08-31', 'COPERNICUS/S2_HARMONIZED') -im7, im8 = collect(aoi, '2022-08-13', '2022-08-31', '2022-09-01', '2022-09-30', 'COPERNICUS/S2_HARMONIZED') -```