Earth_Engine_PCA.ipynb

earth_engine_pca.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyNt6TiVgQLxbeDHzuMrRwwi",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/jdbcode/6304b8cf121148fe5f2b101e24e0cce3/earth_engine_pca.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fSIfBsgi8dNK"
      },
      "source": [
        "#@title Copyright 2024 Google LLC. { display-mode: \"form\" }\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."
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import ee\n",
        "import geemap\n",
        "\n",
        "ee.Authenticate()\n",
        "ee.Initialize(project='MY-PROJECT')"
      ],
      "metadata": {
        "id": "B66Qy8ems9nT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Use these bands.\n",
        "band_names = ee.List(['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11'])\n",
        "\n",
        "# Load a landsat 8 image and select the bands of interest.\n",
        "image = ee.Image('LANDSAT/LC08/C02/T1/LC08_044034_20140318').select(band_names)\n",
        "\n",
        "# Display the input imagery and the region in which to do the PCA.\n",
        "region = image.geometry()\n",
        "m = geemap.Map()\n",
        "m.center_object(region, 10)\n",
        "m.add_layer(ee.Image().paint(region, 0, 2), {}, 'Region')\n",
        "m.add_layer(\n",
        "    image,\n",
        "    {'bands': ['B5', 'B4', 'B2'], 'min': 0, 'max': 20000},\n",
        "    'Original Image',\n",
        ")\n",
        "display(m)\n",
        "\n",
        "# Set an appropriate scale for Landsat data.\n",
        "scale = 30\n",
        "\n",
        "# Mean center the data to enable a faster covariance reducer\n",
        "# and an SD stretch of the principal components.\n",
        "mean_dict = image.reduceRegion(\n",
        "    reducer=ee.Reducer.mean(), geometry=region, scale=scale, maxPixels=1e9\n",
        ")\n",
        "means = mean_dict.toImage(band_names)\n",
        "centered = image.subtract(means)\n",
        "\n",
        "\n",
        "# This helper function returns a list of new band names.\n",
        "def get_new_band_names(prefix):\n",
        "  seq = ee.List.sequence(1, band_names.length())\n",
        "\n",
        "  def add_prefix_and_number(b):\n",
        "    return ee.String(prefix).cat(ee.Number(b).int())\n",
        "\n",
        "  return seq.map(add_prefix_and_number)\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 get_principal_components(centered, scale, region):\n",
        "  # Collapse bands into 1D array\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",
        "  covar_array = ee.Array(covar.get('array'))\n",
        "\n",
        "  # Perform an eigen analysis and slice apart the values and vectors.\n",
        "  eigens = covar_array.eigen()\n",
        "\n",
        "  # This is a P-length vector of Eigenvalues.\n",
        "  eigen_values = eigens.slice(1, 0, 1)\n",
        "  # This is a PxP matrix with eigenvectors in rows.\n",
        "  eigen_vectors = eigens.slice(1, 1)\n",
        "\n",
        "  # Convert the array image to 2D arrays for matrix computations.\n",
        "  array_image = arrays.toArray(1)\n",
        "\n",
        "  # Left multiply the image array by the matrix of eigenvectors.\n",
        "  principal_components = ee.Image(eigen_vectors).matrixMultiply(array_image)\n",
        "\n",
        "  # Turn the square roots of the Eigenvalues into a P-band image.\n",
        "  sd_image = (\n",
        "      ee.Image(eigen_values.sqrt())\n",
        "      .arrayProject([0])\n",
        "      .arrayFlatten([get_new_band_names('sd')])\n",
        "  )\n",
        "\n",
        "  # Turn the PCs into a P-band image, normalized by SD.\n",
        "  return (\n",
        "      # Throw out an an unneeded dimension, [[]] -> [].\n",
        "      principal_components.arrayProject([0])\n",
        "      # Make the one band array image a multi-band image, [] -> image.\n",
        "      .arrayFlatten([get_new_band_names('pc')])\n",
        "      # Normalize the PCs by their SDs.\n",
        "      .divide(sd_image)\n",
        "  )\n",
        "\n",
        "\n",
        "\n",
        "# Get the PCs at the specified scale and in the specified region\n",
        "pc_image = get_principal_components(centered, scale, region)\n",
        "\n",
        "# Plot each PC as a new layer\n",
        "for i, band in enumerate(pc_image.bandNames().getInfo()):\n",
        "  m.add_layer(pc_image.select([band]), {'min': -2, 'max': 2}, band)"
      ],
      "metadata": {
        "id": "MOVk7mp1s1ma"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}