Skip to content

Instantly share code, notes, and snippets.

@mortcanty

mortcanty/imadtutorial.ipynb

Created August 1, 2021 18:07
Show Gist options
  • Save mortcanty/1f5b1924dd655944c9866ec74c657820 to your computer and use it in GitHub Desktop.
Save mortcanty/1f5b1924dd655944c9866ec74c657820 to your computer and use it in GitHub Desktop.
Download ZIP
iMADTutorial.ipynb
Display the source blob
Display the rendered blob
Raw
imadtutorial.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "iMADTutorial.ipynb",
"provenance": [],
"collapsed_sections": [],
"authorship_tag": "ABX9TyMvcsBgX14yeA1LWBsG7TKt",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/mortcanty/1f5b1924dd655944c9866ec74c657820/imadtutorial.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "pJrMq0aNrF-u"
},
"source": [
"#Bitemporal Change Detection on GEE with the iMAD Algorithm\n",
"## Mort Canty\n",
"##April, 2021\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "4xFjivyNFe7U"
},
"source": [
"## Context\n",
"\n",
"The Multivariate Alteration Detection (MAD) algorithm 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 iterative version](https://www2.imm.dtu.dk/pubdb/pubs/4695-full.html) (iMAD). It has since found widespread application for the generation of change information as well as for performing [relative radiometric normalization](http://www2.imm.dtu.dk/pubdb/pubs/5362-full.html) of optical/infrared imagery.\n",
"\n",
"Particularly with the development of the fantastic tools available in Google Earth Engine (GEE) it is easy to generate reflectance time series animations from satellite imagery which show changes quite dramatically. Such visualizations are very good as eye-openers to give a conceptual, qualitative impression of change. This may, however, be perceived differently from person to person. With rigorous statistical methods on the other hand one can derive formal, quantitative and reproducible change information on both pixel and patch/field levels. The iMAD algorithm is one such method. \n",
"\n",
"\n",
"\n",
"This tutorial is based on the description of iMAD given in Chapter 9 of my textbook [Image Anaylsis, Classification and Change Detection in Remote Sensing](https://www.taylorfrancis.com/books/image-analysis-classification-change-detection-remote-sensing-morton-john-canty/10.1201/9780429464348). In preparing it I have tried to take as much advantage as possible of the GEE platform to illustrate and demonstrate the theory interactively .\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "cClM-tAWSnuS"
},
"source": [
"## Preliminaries\n",
"\n",
"The cells below execute the necessary formalities for accessing the Earth Engine and the Google Drive from this Colab notebook. They also import the various Python add-ons needed, including the Folium interactive map package."
]
},
{
"cell_type": "code",
"metadata": {
"id": "ZdyBUBN4F1IA"
},
"source": [
"#@title\n",
"import ee\n",
" \n",
"# Trigger the authentication flow.\n",
"ee.Authenticate()\n",
" \n",
"# Initialize the library.\n",
"ee.Initialize()"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "l1iCQUIDH6wR"
},
"source": [
"from google.colab import drive, files\n",
"drive.mount('/content/drive')"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "X7Z-Enl6Id4n"
},
"source": [
"!ls -l /content/drive/MyDrive/gee"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "NqcghIcEIWnZ"
},
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import random, time\n",
"from osgeo import gdal\n",
"from osgeo.gdalconst import GA_ReadOnly\n",
"import matplotlib.pyplot as plt\n",
"import mpl_toolkits.mplot3d.axes3d as p3\n",
"# Import the Folium library.\n",
"import folium\n",
"\n",
"# Define a method for displaying Earth Engine image tiles to folium map.\n",
"def add_ee_layer(self, ee_image_object, vis_params, name):\n",
" map_id_dict = ee.Image(ee_image_object).getMapId(vis_params)\n",
" folium.raster_layers.TileLayer(\n",
" tiles = map_id_dict['tile_fetcher'].url_format,\n",
" attr = 'Map Data &copy; <a href=\"https://earthengine.google.com/\">Google Earth Engine</a>',\n",
" name = name,\n",
" overlay = True,\n",
" control = True\n",
" ).add_to(self)\n",
"\n",
"# Add EE drawing method to folium.\n",
"folium.Map.add_ee_layer = add_ee_layer"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "8V4fQLWuqTs0"
},
"source": [
"[This folder](https://drive.google.com/drive/folders/1D2RAEu14Zcl9NzqWkSg4J7MilHGHi0BH) folder contains a few of the scripts that will be used in the course of the tutorial."
]
},
{
"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 at 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 examine the difference. Symbolically\n",
" \n",
"$$\n",
"Z = X-Y\n",
"$$\n",
" \n",
"where \n",
"$$\n",
"X= \\pmatrix{X_1\\cr X_2\\cr\\vdots\\cr X_N},\\quad Y= \\pmatrix{Y_1\\cr Y_2\\cr\\vdots\\cr Y_N}\n",
"$$ \n",
" \n",
"are ``typical´´ pixel vectors (more formally: _random vectors_) representing the first and second image, respectively. For example, here is an AOI northwest of Sydney, Australia, where in late March, 2021 extensive flooding occured:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "Xs39zZmNSdru",
"cellView": "form"
},
"source": [
"#@title Area of interest\n",
"geoJSON = {\n",
" \"type\": \"FeatureCollection\",\n",
" \"features\": [\n",
" {\n",
" \"type\": \"Feature\",\n",
" \"properties\": {},\n",
" \"geometry\": {\n",
" \"type\": \"Polygon\",\n",
" \"coordinates\": [\n",
" [\n",
" [\n",
" 150.6562042236328,\n",
" -33.661496156438254\n",
" ],\n",
" [\n",
" 150.9748077392578,\n",
" -33.661496156438254\n",
" ],\n",
" [\n",
" 150.9748077392578,\n",
" -33.40622404437543\n",
" ],\n",
" [\n",
" 150.6562042236328,\n",
" -33.40622404437543\n",
" ],\n",
" [\n",
" 150.6562042236328,\n",
" -33.661496156438254\n",
" ]\n",
" ]\n",
" ]\n",
" }\n",
" }\n",
" ]\n",
"}\n",
"coords = geoJSON['features'][0]['geometry']['coordinates']\n",
"aoi = ee.Geometry.Polygon(coords)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "rwCBEe49d-xx"
},
"source": [
"We first collect two Landsat 8 TOA reflectance scenes"
]
},
{
"cell_type": "code",
"metadata": {
"id": "JzDg-puTSybq",
"cellView": "form"
},
"source": [
"#@title Image collection filters\n",
"im1toa = ee.Image( ee.ImageCollection(\"LANDSAT/LC08/C01/T1_RT_TOA\")\n",
" .filterBounds(aoi) \n",
" .filterDate(ee.Date('2020-10-01'), ee.Date('2021-03-01')) \n",
" .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo')) \n",
" .sort('CLOUD_COVER')\n",
" .first() \n",
" .clip(aoi) )\n",
"im2toa = ee.Image( ee.ImageCollection(\"LANDSAT/LC08/C01/T1_RT_TOA\")\n",
" .filterBounds(aoi) \n",
" .filterDate(ee.Date('2021-03-01'), ee.Date('2021-04-30')) \n",
" .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo')) \n",
" .sort('CLOUD_COVER')\n",
" .first() \n",
" .clip(aoi) )\n",
"timestamp = im1toa.get('system:time_start').getInfo()\n",
"timestamp = time.gmtime(int(timestamp)/1000)\n",
"print(time.strftime('%c', timestamp))\n",
"timestamp = im2toa.get('system:time_start').getInfo()\n",
"timestamp = time.gmtime(int(timestamp)/1000)\n",
"print(time.strftime('%c', timestamp))"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "zekLrceuaNtK"
},
"source": [
" and display them, together with the first three optical bands of the difference image $Z=X-Y$:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "kc_pwScrphK1",
"cellView": "form"
},
"source": [
"#@title Image difference\n",
"location = aoi.centroid().coordinates().getInfo()[::-1]\n",
"m = folium.Map(location=location, zoom_start=11, height=800, width=1000)\n",
"\n",
"percentiles = im1toa.reduceRegion(ee.Reducer.percentile([1,99]),scale=30,maxPixels=1e10) \n",
"mn = percentiles.values(['B2_p1','B3_p1','B4_p1'])\n",
"mx = percentiles.values(['B2_p99','B3_p99','B4_p99'])\n",
"rgb = ee.Image.rgb(im1toa.select('B2'),im1toa.select('B3'),im1toa.select('B4'))\n",
"m.add_ee_layer(rgb, {'min': mn, 'max': mx}, '24-09-2018')\n",
"\n",
"percentiles = im2toa.reduceRegion(ee.Reducer.percentile([1,99]),scale=30,maxPixels=1e10) \n",
"mn = percentiles.values(['B2_p1','B3_p1','B4_p1'])\n",
"mx = percentiles.values(['B2_p99','B3_p99','B4_p99'])\n",
"rgb = ee.Image.rgb(im2toa.select('B2'),im2toa.select('B3'),im2toa.select('B4'))\n",
"m.add_ee_layer(rgb, {'min': mn, 'max': mx}, '20-04-2019')\n",
"\n",
"imd = im1toa.subtract(im2toa)\n",
"percentiles = imd.reduceRegion(ee.Reducer.percentile([1,99]),scale=30,maxPixels=1e10) \n",
"mn = percentiles.values(['B2_p1','B3_p1','B4_p1'])\n",
"mx = percentiles.values(['B2_p99','B3_p99','B4_p99'])\n",
"rgb = ee.Image.rgb(imd.select('B2'),imd.select('B3'),imd.select('B4'))\n",
"m.add_ee_layer(rgb, {'min': mn, 'max': mx}, 'Difference bands 2,3,4')\n",
"\n",
"# Add a layer control panel to the map.\n",
"m.add_child(folium.LayerControl())\n",
"\n",
"# Display the map.\n",
"display(m)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "h2d-vg-ezwpe"
},
"source": [
"The most obvious changes are due to the flooding as well as to cloud (and cloud shadow) movement. Small intensity differences indicate no change, large\n",
"positive or negative values indicate change, and decision thresholds\n",
"can be set to define significance. The thresholds are usually expressed\n",
"in terms of standard deviations from the mean difference value, which is taken to correspond to no change (the dark gray background in the difference image above). If the detected signals are uncorrelated, the variances of the difference\n",
"images are simply\n",
"$$\n",
"{\\rm var}(Z_i) = {\\rm var}(X_i) + {\\rm var}(Y_i),\\quad i=1\\dots N,\n",
"$$\n",
"or about twice as noisy as the individual image bands. 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_."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "rhMR6KLhkar5"
},
"source": [
"## The Multivariate Alteration Detection (MAD) Transformation\n",
"\n",
"While the inundated areas are easily identified in the simple difference images, it is possible to derive a better and more quantitative change map. This involves taking greater advantage of the statistical properties of the images.\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 = a^\\top X.\n",
"$$\n",
"\n",
"The symbol $a^\\top$ denotes the transpose of the vector $a$, i.e., the row vector $(a_1,a_2 \\dots a_N)$ and the expression $a^\\top X$ is an _inner vector product_. The vector of coefficients $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 = b^\\top Y,\n",
"$$\n",
"\n",
"and then look at the scalar difference $U-V$. Change information\n",
"is now contained in a single image rather than distributed among all $N$ bands.\n",
"\n",
"One has, of course, to choose the coefficient vectors $a$ and\n",
"$b$ in some suitable way. In [Nielsen et al. (1998)](https://www2.imm.dtu.dk/pubdb/pubs/1220-full.html) it is suggested\n",
"that they be determined by making the transformed images $U$ and $V$ _as similar as possible_ before taking their difference. This can be accomplished with standard _canonical correlation analysis_ (CCA), a technique first described by Harold Hotelling in 1936. When performing CCA on a bi-temporal image, one maximizes the correlation $\\rho$ between the random variables $U$ and $V$, given by\n",
"\n",
"$$\n",
"\\rho = {{\\rm cov}(U,V)\\over \\sqrt{{\\rm var}(U)}\\sqrt{{\\rm var}(V)}}. \\tag{1}\n",
"$$\n",
"\n",
"Arbitrary multiples of $U$ and $V$ would clearly have the same correlation (they cancel off in numerator and denominator),\n",
"so a constraint must be chosen. A convenient one is\n",
"\n",
"$$\n",
"{\\rm var}(U)={\\rm var}(V)=1. \\tag{2}\n",
"$$\n",
"\n",
"The details of CCA are given [my textbook](https://www.taylorfrancis.com/books/image-analysis-classification-change-detection-remote-sensing-morton-john-canty/10.1201/9780429464348), but it all boils down to an eigenvalue problem for determining the transformation vectors $a$ and $b$. The solution consists of $N$ sets of eigenvectors $(a^i, b^i),\\ i=1\\dots N$ and, correspondingly, $N$ sets of so-called _canonical variates_ $(U_i, V_i),\\ i= 1\\dots N$. Unlike the bands of the original images $X$ and $Y$, they are ordered by similarity rather than wavelength. The difference images \n",
"\n",
"$$\n",
"M_i = U_i - V_i = (a^i)^\\top X - (b^i)^\\top Y,\\quad i=1\\dots N, \\tag{3}\n",
"$$\n",
"\n",
"are called the MAD variates. \n",
"\n",
"The cononical and MAD variates have very nice properties indeed. The former are _all mutually uncorrelated_ except for the\n",
"pairs $(U_i,V_i)$, and these are ordered by decreasing correlation:\n",
"\n",
"$$\n",
"{\\rm cov}(U_i, V_i) = \\rho_i,\\quad i=1\\dots N, \\tag{4}\n",
"$$\n",
"\n",
"with $\\rho_1 > \\rho_2 >\\dots >\\rho_N$.\n",
"\n",
"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{5}\n",
"$$\n",
"\n",
"and their variances by\n",
"\n",
"$$\n",
"\\sigma_{M_i}^2={\\rm var}(U_i-V_i)=2(1-\\rho_i),\\quad i=1\\dots N. \\tag{6}\n",
"$$\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 but 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\n",
"bands."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "Td12EFgEyB7C"
},
"source": [
"## Iterative re-weighting: The iMAD algorithm\n",
" \n",
"Let's imagine two images of the same scene, acquired at different times\n",
"under similar conditions, but for which no ground reflectance\n",
"changes have occurred whatsoever. Then the only differences\n",
"between them will be due to random effects like instrument noise\n",
"and atmospheric fluctuation. In such a case we would expect that\n",
"the histogram of any difference component that we generate will\n",
"be very nearly Gaussian. In particular, the MAD variates, being\n",
"uncorrelated, should follow a multivariate, zero mean normal\n",
"distribution with diagonal covariance matrix. Change observations\n",
"would deviate more or less strongly from such a distribution.\n",
"We might therefore expect an improvement in the sensitivity of the MAD\n",
"transformation if we can establish a better background\n",
"of no change against which to detect change.\n",
"This can be done in an iteration scheme [(Nielsen 2007)](https://www2.imm.dtu.dk/pubdb/pubs/4695-full.html) in which, when\n",
"calculating the statistics for each successive iteration of the MAD transformation, observations are weighted in some appropriate fashion.\n",
" \n",
"Let the variable $Z$ represent 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, \\tag{7}\n",
"$$\n",
" \n",
"where $\\sigma_{M_i}$ is given by Equation (6).\n",
"Then, since the no-change observations are expected to be\n",
"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. In fact it is easy to show that $Z$ is a _likelihood ratio test_ for change, see again [my textbook](https://www.taylorfrancis.com/books/image-analysis-classification-change-detection-remote-sensing-morton-john-canty/10.1201/9780429464348), Chapter 9 and the discussion of statistical hypothesis testing in my [Google Earth Engine Tutorial.](https://developers.google.com/earth-engine/tutorials/community/detecting-changes-in-sentinel-1-imagery-pt-2)\n",
" \n",
"The so-called $p$-_value_ for an observation of the test statistic $Z$ is defined to be the probability that a sample $z$ drawn from the chi-square distribution could be that large or larger. This is given by\n",
" \n",
"$$\n",
"p(z) = 1-P_{\\chi^2;N}(z),\n",
"$$\n",
" \n",
"where $P_{\\chi^2;N}$ is the cumulative chi square probability distribution.\n",
"Small $p$-values are unlikely if no change has occured at a given pixel location and so small $p$-values indicate change. Therefore, after each iteration, 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 gradually reduces the influence of the change observations on the MAD transformation. Iteration continues until some stopping criterion is met, such as lack of significant change in the canonical correlations. The whole procedure constitutes the iMAD algorithm.\n",
" \n",
"The following input cell is an implementation of the iMAD algorithm written in the GEE Python API:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "n6e3lvyZKsNG",
"cellView": "form"
},
"source": [
"#@title The iMAD code\n",
"import ee\n",
"\n",
"def covarw(image, weights, scale):\n",
" '''Return the weighted centered image and its weighted covariance matrix''' \n",
" geometry = image.geometry()\n",
" bandNames = image.bandNames()\n",
" N = bandNames.length()\n",
" weightsImage = image.multiply(ee.Image.constant(0)).add(weights)\n",
" means = image.addBands(weightsImage) \\\n",
" .reduceRegion(ee.Reducer.mean().repeat(N).splitWeights(), scale=scale, maxPixels=1e10) \\\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(), scale=scale).get(B1)) \n",
" sumWeights = ee.Number(weights.reduceRegion(ee.Reducer.sum(),geometry=geometry, scale=scale).get(b1))\n",
" covw = centered.multiply(weights.sqrt()) \\\n",
" .toArray() \\\n",
" .reduceRegion(ee.Reducer.centeredCovariance(), geometry=geometry, scale=scale) \\\n",
" .get('array')\n",
" covw = ee.Array(covw).multiply(nPixels).divide(sumWeights)\n",
" return (centered.arrayFlatten([bandNames]), covw)\n",
"\n",
"def chi2cdf(chi2,df):\n",
" '''Chi square cumulative distribution function with df degrees of freedom'''\n",
" return ee.Image(chi2.divide(2)).gammainc(ee.Number(df).divide(2))\n",
"\n",
"def geneiv(C,B):\n",
" '''Return the eignvalues and eigenvectors of the Generalized eigenproblem\n",
" C*X = lambda*B*X '''\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 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",
"\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",
" chi2 = ee.Image(ee.Dictionary(prev).get('chi2')) \n",
" allrhos = ee.List(ee.Dictionary(prev).get('allrhos'))\n",
" region = image.geometry() \n",
" nBands = image.bandNames().length().divide(2) \n",
" weights = chi2cdf(chi2,nBands).subtract(1).multiply(-1)\n",
" scale = ee.Number(ee.Dictionary(prev).get('scale'))\n",
" niter = ee.Number(ee.Dictionary(prev).get('niter'))\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.001)) \\\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",
" MAD = U.subtract(V)\n",
"# chi square image\n",
" chi2 = MAD.pow(2) \\\n",
" .divide(sigma2s) \\\n",
" .reduce(ee.Reducer.sum()) \\\n",
" .clip(region) \n",
" return ee.Dictionary({'done':done,'scale':scale,'niter':niter.add(1),'image':image,'allrhos':allrhos,'chi2':chi2,'MAD':MAD})"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "d6FI_gi5UiKC"
},
"source": [
"and here we run it on the optical/infrared bands of our image pair:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "hyxPR_GRZ04i",
"cellView": "form"
},
"source": [
"#@title Run iMAD algorithm \n",
"%%time\n",
"alpha = 0.0001\n",
"im1 = im1toa.select('B2','B3','B4','B5','B6','B7')\n",
"im2 = im2toa.select('B2','B3','B4','B5','B6','B7')\n",
"madnames = ['MAD'+str(i+1) for i in range(6)]\n",
"inputlist = ee.List.sequence(1,100)\n",
"first = ee.Dictionary({'done':ee.Number(0),\n",
" 'scale':ee.Number(60),\n",
" 'niter':ee.Number(0),\n",
" 'image':im1.addBands(im2).clip(aoi),\n",
" 'allrhos': [ee.List.sequence(1,6)],\n",
" 'chi2':ee.Image.constant(0),\n",
" 'MAD':ee.Image.constant(0)}) \n",
"result = ee.Dictionary(inputlist.iterate(imad,first)) \n",
"MAD = ee.Image(result.get('MAD')).rename(madnames)\n",
"chi2 = ee.Image(result.get('chi2')).rename(['chi2']) \n",
"pval = chi2cdf(chi2,6).subtract(1).multiply(-1).rename('pval')\n",
"nochange = pval.gt(alpha)\n",
" \n",
"print('Iterations: %i'%ee.Number(result.get('niter')).getInfo())\n",
"print('Canonical correlations:\\n%s'%str(ee.List(result.get('allrhos')).get(-1).getInfo()))"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "9iUiafXLU-7H"
},
"source": [
"The next cell displays the first three iMAD bands in RGB and also shows the chi square image, Equation (7) and the no change pixels at _significance level_ $\\alpha=0.0001$, i.e., those pixels for which \n",
"\n",
"$$\n",
"p(z) > 0.0001.\n",
"$$"
]
},
{
"cell_type": "code",
"metadata": {
"id": "GGM_81zQTu0Q",
"cellView": "form"
},
"source": [
"#@title Display the results of the iMAD transformation\n",
"location = aoi.centroid().coordinates().getInfo()[::-1]\n",
"m = folium.Map(location=location, zoom_start=11, height=800, width=1000)\n",
"\n",
"# overlay MAD, chi square and no change\n",
"#percentiles = MAD.reduceRegion(ee.Reducer.percentile([1,99]),scale=30) \n",
"MAD.updateMask(nochange)\n",
"rgbd = ee.Image.rgb(MAD.select(1),MAD.select(0),MAD.select(2))\n",
"m.add_ee_layer(rgbd, {'min': -10, 'max': 10}, 'iMAD 2,1,3')\n",
"m.add_ee_layer(chi2, {'min': 0, 'max': 20000}, 'chi2 image')\n",
"m.add_ee_layer(nochange, {'min': 0, 'max': 1}, 'nochange')\n",
"\n",
"\n",
"# Add a layer control panel to the map.\n",
"m.add_child(folium.LayerControl())\n",
"\n",
"# Display the map.\n",
"display(m)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "aU6UonxyaeMt"
},
"source": [
"Middle gray pixels in the iMAD RGB representation above again point to no change, the bright colored pixels mark the inundated areas. Notice that the river bed also exhibits a strong (black) change signal, due to the silting or muddying of the water. This was also apparent in the simple difference image."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "cJlCMSywTHcE"
},
"source": [
"\n",
"## Unsupervised clustering\n",
"In order to generate a more easily interpetable change map we will next apply an _unsupervised clustering algorithm_ to the iMAD image. It will generally be necessary to experiment with the number of clusters, because clustering forces a partitioning of the change/no-change pixels which may or may not be natural. We prefer to use the _Expectation Maximization_ algorithm described in Chapter 8 of [Image Anaylsis, Classification and Change Detection im Remote Sensing](https://www.taylorfrancis.com/books/image-analysis-classification-change-detection-remote-sensing-morton-john-canty/10.1201/9780429464348) because it includes multi-resolution, simuated annealing and spatial class membership. It is not available on GEE but is coded in the program _em.py_ in the _scripts_ folder. In order to use it, we have to download the iMAD image to our Google Drive, calling it _MADWindsor_ (the town of Windsor got the brunt of the flood):"
]
},
{
"cell_type": "code",
"metadata": {
"id": "TiqC9Oe51IBx"
},
"source": [
"drexport = ee.batch.Export.image.toDrive(MAD,\n",
" description='driveExportTask', \n",
" folder = 'gee',\n",
" fileNamePrefix='MADWindsor',scale=30,maxPixels=1e10)\n",
"drexport.start()"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "VC7sU4XNW5j9"
},
"source": [
"Here is the result for 6 clusters and an $800\\times 800$ pixel spatial subset of the iMAD image:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "J5BuH5dG4n0v"
},
"source": [
"%run /content/drive/MyDrive/scripts/em -K 6 -d [100,100,800,800] /content/drive/MyDrive/gee/MADWindsor.tif"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "yoBfFQs9XnZp"
},
"source": [
"Also included in the _scripts_ folder is the program _dispms.py_ for displaying multispectral images and classification maps so we'll use it to view the change map:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "FZ1QRfm88lnH"
},
"source": [
"%run /content/drive/MyDrive/scripts/dispms -f /content/drive/MyDrive/gee/MADWindsor_em.tif -c \\\n",
"-r \"['no change','2','3','4','5','6']\""
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "SY9jtyB3YPK9"
},
"source": [
"Thus we have managed to quantize the change image into six categories. Evidently the yellow and orange clusters (classes 4 and 5) correspond to the flooding whereas most cloud changes, for example, are in the the sixth cluster (red). Here is some quick and dirty code to display the various clusters in a 3-dimensional plot, in this instance classes 1 (black = no change) and 5 (orange = change):"
]
},
{
"cell_type": "code",
"metadata": {
"id": "hJ2Pu-YW2vsP",
"cellView": "form"
},
"source": [
"#@title 3D Cluster Plot\n",
"\n",
"''' Make a 3d plot of em clusters '''\n",
"\n",
"gdal.AllRegister()\n",
"num = 2000\n",
"rows = 800\n",
"cols = 800\n",
"x0 = 100\n",
"y0 = 100\n",
"clusters = [1,5]\n",
"\n",
"# read the first 3 iMAD bands into array imad\n",
"inDataset = gdal.Open('/content/drive/MyDrive/gee/MADWindsor.tif',GA_ReadOnly)\n",
"imad = np.zeros((rows*cols,3))\n",
"for b in range(3):\n",
" band = inDataset.GetRasterBand(b+1)\n",
" imad[:,b] = band.ReadAsArray(x0,y0,rows,cols).astype(float).ravel()\n",
"inDataset = None \n",
"# and the classified image into array em \n",
"inDataset = gdal.Open('/content/drive/MyDrive/gee/MADWindsor_em.tif',GA_ReadOnly)\n",
"band = inDataset.GetRasterBand(1)\n",
"em = band.ReadAsArray(0,0,cols,rows).ravel()\n",
"inDataset = None\n",
"k = np.max(em)\n",
"if clusters==None:\n",
" clusters = range(0,min(5,k))\n",
"else:\n",
" for i in range(len(clusters)):\n",
" clusters[i] -= 1 \n",
"fig = plt.figure(figsize=(10, 8))\n",
"ax = p3.Axes3D(fig) \n",
"ax.set_xlabel('MAD1')\n",
"ax.set_ylabel('MAD2')\n",
"ax.set_zlabel('MAD3')\n",
"colors = ['black','blue','cyan','yellow','orange','brown']\n",
"labels = ['no change','change2','change3','change4','change5','change6']\n",
"for i in clusters: \n",
" idx = np.where(em==(i+1))[0] \n",
" random.shuffle(idx) \n",
" idx = idx[:num]\n",
" ax.scatter(imad[idx,0], imad[idx,1], imad[idx,2], s=10, c=colors[i], label=labels[i])\n",
"plt.title('iMAD clusters')\n",
"plt.legend(loc=2)\n",
"plt.show()"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "gMWvdYYhZeYx"
},
"source": [
"Now the advantage of the iMAD algorithm can be clearly seen: the no-change cluster is small and very compact and the changes are well discriminated. \n",
"\n",
"Finally, after downloading the cluster image to my own PC and ingesting it into my GEE assets, here is an overlay of the flood classes onto a folium map. You can do the same with your own assets directory:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "y2Z8uIWzP5Qo"
},
"source": [
"location = aoi.centroid().coordinates().getInfo()[::-1]\n",
"m = folium.Map(location=location, zoom_start=11, height=800, width=1000)\n",
"\n",
"# overlay cluster image, clusters 4 and 5 only\n",
"im = ee.Image('users/mortcanty/tutorial/MADWindsor_em')\n",
"im = im.updateMask(im.eq(5).Or(im.eq(4)))\n",
"\n",
"m.add_ee_layer(im,{'min':1,'max':6, 'palette':['black','blue','cyan','yellow','orange','brown']},'MAD_EM')\n",
"m.add_child(folium.LayerControl())\n",
"\n",
"display(m)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "kv9WzmXdoXgW"
},
"source": [
"## Kangaroo Island\n",
"Here is another area of interest, this time Kangaroo Island, also Australia, where in January, 2020, bush fires destroyed almost 50% of the island :"
]
},
{
"cell_type": "code",
"metadata": {
"id": "yxQjqBUMofmA",
"cellView": "form"
},
"source": [
"#@title Area of interest\n",
"geoJSON = {\n",
" \"type\": \"FeatureCollection\",\n",
" \"features\": [\n",
" {\n",
" \"type\": \"Feature\",\n",
" \"properties\": {},\n",
" \"geometry\": {\n",
" \"type\": \"Polygon\",\n",
" \"coordinates\": [\n",
" [\n",
" [\n",
" 136.53533935546875,\n",
" -36.07352228885536\n",
" ],\n",
" [\n",
" 138.1365966796875,\n",
" -36.07352228885536\n",
" ],\n",
" [\n",
" 138.1365966796875,\n",
" -35.55680897384459\n",
" ],\n",
" [\n",
" 136.53533935546875,\n",
" -35.55680897384459\n",
" ],\n",
" [\n",
" 136.53533935546875,\n",
" -36.07352228885536\n",
" ]\n",
" ]\n",
" ]\n",
" }\n",
" }\n",
" ]\n",
"}\n",
"coords = geoJSON['features'][0]['geometry']['coordinates']\n",
"aoi = ee.Geometry.Polygon(coords)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "RogcKTmn_j5s"
},
"source": [
"We collect two images, one from January, 2019 and one from February, 2020."
]
},
{
"cell_type": "code",
"metadata": {
"id": "ocjSpGVbpCUr",
"cellView": "form"
},
"source": [
"#@title Image collection filters\n",
"im1toa = ee.Image( ee.ImageCollection(\"LANDSAT/LC08/C01/T1_RT_TOA\")\n",
" .filterBounds(aoi) \n",
" .filterDate(ee.Date('2019-01-15'), ee.Date('2019-03-01')) \n",
" .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo')) \n",
" .sort('CLOUD_COVER')\n",
" .first() \n",
" .clip(aoi) )\n",
"im2toa = ee.Image( ee.ImageCollection(\"LANDSAT/LC08/C01/T1_RT_TOA\")\n",
" .filterBounds(aoi) \n",
" .filterDate(ee.Date('2020-01-15'), ee.Date('2020-03-01')) \n",
" .filter(ee.Filter.contains(rightValue=aoi,leftField='.geo')) \n",
" .sort('CLOUD_COVER')\n",
" .first() \n",
" .clip(aoi) )\n",
"timestamp = im1toa.get('system:time_start').getInfo()\n",
"timestamp = time.gmtime(int(timestamp)/1000)\n",
"print(time.strftime('%c', timestamp))\n",
"timestamp = im2toa.get('system:time_start').getInfo()\n",
"timestamp = time.gmtime(int(timestamp)/1000)\n",
"print(time.strftime('%c', timestamp))"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "1oQKOiX1_9VB"
},
"source": [
"and, as before, display the difference image:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "YKkOq66wqpSf",
"cellView": "form"
},
"source": [
"#@title Image difference\n",
"location = aoi.centroid().coordinates().getInfo()[::-1]\n",
"m = folium.Map(location=location, zoom_start=10, height=800, width=1000)\n",
"\n",
"percentiles = im1toa.reduceRegion(ee.Reducer.percentile([1,99]),scale=30,maxPixels=1e10) \n",
"mn = percentiles.values(['B2_p1','B3_p1','B4_p1'])\n",
"mx = percentiles.values(['B2_p99','B3_p99','B4_p99'])\n",
"rgb = ee.Image.rgb(im1toa.select('B2'),im1toa.select('B3'),im1toa.select('B4'))\n",
"m.add_ee_layer(rgb, {'min': mn, 'max': mx}, '22-01-2019')\n",
"\n",
"percentiles = im2toa.reduceRegion(ee.Reducer.percentile([1,99]),scale=30,maxPixels=1e10) \n",
"mn = percentiles.values(['B2_p1','B3_p1','B4_p1'])\n",
"mx = percentiles.values(['B2_p99','B3_p99','B4_p99'])\n",
"rgb = ee.Image.rgb(im2toa.select('B2'),im2toa.select('B3'),im2toa.select('B4'))\n",
"m.add_ee_layer(rgb, {'min': mn, 'max': mx}, '10-02-2020')\n",
"\n",
"imd = im1toa.subtract(im2toa)\n",
"percentiles = imd.reduceRegion(ee.Reducer.percentile([1,99]),scale=30,maxPixels=1e10) \n",
"mn = percentiles.values(['B2_p1','B3_p1','B4_p1'])\n",
"mx = percentiles.values(['B2_p99','B3_p99','B4_p99'])\n",
"rgb = ee.Image.rgb(imd.select('B2'),imd.select('B3'),imd.select('B4'))\n",
"m.add_ee_layer(rgb, {'min': mn, 'max': mx}, 'Difference bands 2,3,4')\n",
"\n",
"# Add a layer control panel to the map.\n",
"m.add_child(folium.LayerControl())\n",
"\n",
"# Display the map.\n",
"display(m)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "vHTJZ6E2vwSE"
},
"source": [
"Because of the size of the AOI we will need to run iMAD as an export task and, moreover, we will mask out the sea:"
]
},
{
"cell_type": "code",
"metadata": {
"id": "JKciF0AVqG1b",
"cellView": "form"
},
"source": [
"#@title Run iMAD algorithm \n",
"\n",
"water_mask = ee.Image('UMD/hansen/global_forest_change_2015') \\\n",
" .select('datamask').eq(1)\n",
"\n",
"alpha = 0.01\n",
"im1 = im1toa.select('B2','B3','B4','B5','B6','B7')\n",
"im2 = im2toa.select('B2','B3','B4','B5','B6','B7')\n",
"madnames = ['MAD'+str(i+1) for i in range(6)]\n",
"inputlist = ee.List.sequence(1,100)\n",
"first = ee.Dictionary({'done':ee.Number(0),\n",
" 'scale':ee.Number(60),\n",
" 'niter':ee.Number(0),\n",
" 'image':im1.addBands(im2).clip(aoi).updateMask(water_mask),\n",
" 'allrhos': [ee.List.sequence(1,6)],\n",
" 'chi2':ee.Image.constant(0),\n",
" 'MAD':ee.Image.constant(0)}) \n",
"result = ee.Dictionary(inputlist.iterate(imad,first)) \n",
"MAD = ee.Image(result.get('MAD')).rename(madnames)\n",
" \n",
"assetId='users/mortcanty/tutorial/kgroo',\n",
"assexport = ee.batch.Export.image.toAsset(MAD,\n",
" description='assetExportTask', \n",
" assetId='users/mortcanty/tutorial/kgroo',\n",
" scale=30, \n",
" maxPixels=1e10)\n",
"assexport.start()\n",
"\n",
"print('Exporting change map to %s'%assetId)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "pRjEoaKKARn4"
},
"source": [
"Displaying an RGB composite of three MAD components, the burned over areas are clearly visible in yellow/pink colors."
]
},
{
"cell_type": "code",
"metadata": {
"id": "e95GJPWq2bAE",
"cellView": "form"
},
"source": [
"#@title Display the results of the iMAD transformation\n",
"location = aoi.centroid().coordinates().getInfo()[::-1]\n",
"m = folium.Map(location=location, zoom_start=9, height=800, width=1000)\n",
"MAD = ee.Image('users/mortcanty/tutorial/kgroo')\n",
"# overlay MAD\n",
"rgbd = ee.Image.rgb(MAD.select('MAD4'),MAD.select('MAD2'),MAD.select('MAD5'))\n",
"m.add_ee_layer(rgbd, {'min': -40, 'max': 40}, 'iMAD 4,2,5')\n",
"\n",
"# Add a layer control panel to the map.\n",
"m.add_child(folium.LayerControl())\n",
"\n",
"# Display the map.\n",
"display(m)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "IwzNfbpmcLI1"
},
"source": [
"## Conclusion and Outlook\n",
"While simple image comparison or differencing can be useful, the statistical transformations and iterations implicit in the iMAD algorithm offer a more powerful and quantifiable means of analyzing and categorizing changes in bitemporal image data. The main caveat is that, in some instances, the iMAD iteration will not converge simply because there are too few no-change pixels in the scene. A remedy is to focus on a spatial subset with relatively more no-change observations and then, if the algorithm converges, use the resulting eigenvectors $a^i, b^i$ to transform the full images.\n",
"\n",
"In a future tutorial we will have a look at _relative radiometric image normalization_ using the no-change pixels identified by the iMAD algorithm."
]
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment