tutorials/imad-tutorial-pt2/index.ipynb

{
  "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 in Google Earth Engine - The MAD Transformation (Part 2)\n",
        "Author: 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"
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": [
        "This tutorial will optionally export assets. Edit the following\n",
        "`EXPORT_PATH` variable to specify the location to store the assets.\n",
        "All assets that are needed to complete the tutorial are hosted by Earth Engine,\n",
        "but if you'd like to display assets that you export, replace paths as needed."
      ]
    },
    {
      "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(\n",
        "    'projects/google/imad_tutorial/landkreis_olpe_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. Earth Engine\n",
        "provides the asset for use in this demo, so it is not required that you run the\n",
        "export cell to complete the tutorial."
      ]
    },
    {
      "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, assetId, 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",
        "        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": [
        "asset_path = f'{EXPORT_PATH}LandkreisOlpe'\n",
        "run_imad(aoi, im1.select(visirbands), im2.select(visirbands), asset_path)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "21061d67",
      "metadata": {
        "id": "21061d67"
      },
      "source": [
        "After the export finishes (if you ran it), the number of iterations and the\n",
        "final canonical correlations can be read from properties of the exported image.\n",
        "(Earth Engine-hosted assets are imported here, edit the paths if you'd like to\n",
        "use your copies)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "a67a360d",
      "metadata": {
        "id": "a67a360d"
      },
      "outputs": [],
      "source": [
        "im_imad = ee.Image(\n",
        "    'projects/google/imad_tutorial/LandkreisOlpe').select(0, 1, 2, 3, 4, 5)\n",
        "im_z = ee.Image(\n",
        "    'projects/google/imad_tutorial/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
}