Array/eigen_analysis.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<table class=\"ee-notebook-buttons\" align=\"left\">\n",
        "    <td><a target=\"_blank\"  href=\"https://github.com/giswqs/earthengine-py-notebooks/tree/master/Array/eigen_analysis.ipynb\"><img width=32px src=\"https://www.tensorflow.org/images/GitHub-Mark-32px.png\" /> View source on GitHub</a></td>\n",
        "    <td><a target=\"_blank\"  href=\"https://nbviewer.jupyter.org/github/giswqs/earthengine-py-notebooks/blob/master/Array/eigen_analysis.ipynb\"><img width=26px src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/3/38/Jupyter_logo.svg/883px-Jupyter_logo.svg.png\" />Notebook Viewer</a></td>\n",
        "    <td><a target=\"_blank\"  href=\"https://colab.research.google.com/github/giswqs/earthengine-py-notebooks/blob/master/Array/eigen_analysis.ipynb\"><img src=\"https://www.tensorflow.org/images/colab_logo_32px.png\" /> Run in Google Colab</a></td>\n",
        "</table>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Install Earth Engine API and geemap\n",
        "Install the [Earth Engine Python API](https://developers.google.com/earth-engine/python_install) and [geemap](https://geemap.org). The **geemap** Python package is built upon the [ipyleaflet](https://github.com/jupyter-widgets/ipyleaflet) and [folium](https://github.com/python-visualization/folium) packages and implements several methods for interacting with Earth Engine data layers, such as `Map.addLayer()`, `Map.setCenter()`, and `Map.centerObject()`.\n",
        "The following script checks if the geemap package has been installed. If not, it will install geemap, which automatically installs its [dependencies](https://github.com/giswqs/geemap#dependencies), including earthengine-api, folium, and ipyleaflet."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Installs geemap package\n",
        "import subprocess\n",
        "\n",
        "try:\n",
        "    import geemap\n",
        "except ImportError:\n",
        "    print('Installing geemap ...')\n",
        "    subprocess.check_call([\"python\", '-m', 'pip', 'install', 'geemap'])"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "import ee\n",
        "import geemap"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Create an interactive map \n",
        "The default basemap is `Google Maps`. [Additional basemaps](https://github.com/giswqs/geemap/blob/master/geemap/basemaps.py) can be added using the `Map.add_basemap()` function. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "Map = geemap.Map(center=[40,-100], zoom=4)\n",
        "Map"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Add Earth Engine Python script "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Add Earth Engine dataset\n",
        "# Compute the Principal Components of a Landsat 8 image.\n",
        "\n",
        "\n",
        "# Load a landsat 8 image, select the bands of interest.\n",
        "image = ee.Image('LANDSAT/LC8_L1T/LC80440342014077LGN00') \\\n",
        "  .select(['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11'])\n",
        "\n",
        "# Display the input imagery and the region in which to do the PCA.\n",
        "region = image.geometry()\n",
        "Map.centerObject(ee.FeatureCollection(region), 10)\n",
        "Map.addLayer(ee.Image().paint(region, 0, 2), {}, 'Region')\n",
        "Map.addLayer(image, {'bands': ['B5', 'B4', 'B2'], 'min': 0, 'max': 20000}, 'Original Image')\n",
        "\n",
        "# Set some information about the input to be used later.\n",
        "scale = 30\n",
        "bandNames = image.bandNames()\n",
        "\n",
        "# Mean center the data to enable a faster covariance reducer\n",
        "# and an SD stretch of the principal components.\n",
        "meanDict = image.reduceRegion(**{\n",
        "    'reducer': ee.Reducer.mean(),\n",
        "    'geometry': region,\n",
        "    'scale': scale,\n",
        "    'maxPixels': 1e9\n",
        "})\n",
        "means = ee.Image.constant(meanDict.values(bandNames))\n",
        "centered = image.subtract(means)\n",
        "\n",
        "# This helper function returns a list of new band names.\n",
        "def getNewBandNames(prefix):\n",
        "  seq = ee.List.sequence(1, bandNames.length())\n",
        "  return seq.map(lambda b: ee.String(prefix).cat(ee.Number(b).int().format()))\n",
        "\n",
        "\n",
        "# This function accepts mean centered imagery, a scale and\n",
        "# a region in which to perform the analysis.  It returns the\n",
        "# Principal Components (PC) in the region as a new image.\n",
        "def getPrincipalComponents(centered, scale, region):\n",
        "  # Collapse the bands of the image into a 1D array per pixel.\n",
        "  arrays = centered.toArray()\n",
        "\n",
        "  # Compute the covariance of the bands within the region.\n",
        "  covar= arrays.reduceRegion(**{\n",
        "    'reducer': ee.Reducer.centeredCovariance(),\n",
        "    'geometry': region,\n",
        "    'scale': scale,\n",
        "    'maxPixels': 1e9\n",
        "  })\n",
        "\n",
        "  # Get the 'array' covariance result and cast to an array.\n",
        "  # This represents the band-to-band covariance within the region.\n",
        "  covarArray = ee.Array(covar.get('array'))\n",
        "\n",
        "  # Perform an eigen analysis and slice apart the values and vectors.\n",
        "  eigens = covarArray.eigen()\n",
        "\n",
        "  # This is a P-length vector of Eigenvalues.\n",
        "  eigenValues = eigens.slice(1, 0, 1)\n",
        "  # This is a PxP matrix with eigenvectors in rows.\n",
        "  eigenVectors = eigens.slice(1, 1)\n",
        "\n",
        "  # Convert the array image to 2D arrays for matrix computations.\n",
        "  arrayImage = arrays.toArray(1)\n",
        "\n",
        "  # Left multiply the image array by the matrix of eigenvectors.\n",
        "  principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage)\n",
        "\n",
        "  # Turn the square roots of the Eigenvalues into a P-band image.\n",
        "  sdImage = ee.Image(eigenValues.sqrt()) \\\n",
        "    .arrayProject([0]).arrayFlatten([getNewBandNames('sd')])\n",
        "\n",
        "  # Turn the PCs into a P-band image, normalized by SD.\n",
        "  return principalComponents \\\n",
        "    .arrayProject([0]) \\\n",
        "    .arrayFlatten([getNewBandNames('pc')]) \\\n",
        "    .divide(sdImage) \\\n",
        "\n",
        "\n",
        "# Get the PCs at the specified scale and in the specified region\n",
        "pcImage = getPrincipalComponents(centered, scale, region)\n",
        "\n",
        "Map.addLayer(pcImage.select(0), {}, 'Image')\n",
        "\n",
        "for i in range(0, bandNames.length().getInfo()):\n",
        "    band = pcImage.bandNames().get(i).getInfo()\n",
        "    Map.addLayer(pcImage.select([band]), {'min': -2, 'max': 2}, band)\n",
        "\n",
        "\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Display Earth Engine data layers "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "Map.addLayerControl() # This line is not needed for ipyleaflet-based Map.\n",
        "Map"
      ],
      "outputs": [],
      "execution_count": null
    }
  ],
  "metadata": {
    "anaconda-cloud": {},
    "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.6.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}