8_Array.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/csaybar/EEwPython/blob/dev/8_Array.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_rUp3DC76g8u"
   },
   "source": [
    "<!--COURSE_INFORMATION-->\n",
    "<img align=\"left\" style=\"padding-right:10px;\" src=\"https://sitejerk.com/images/google-earth-logo-png-5.png\" width=5% >\n",
    "<img align=\"right\" style=\"padding-left:10px;\" src=\"https://colab.research.google.com/img/colab_favicon_256px.png\" width=6% >\n",
    "\n",
    "\n",
    ">> *This notebook is part of the free course [EEwPython](https://colab.research.google.com/github/csaybar/EEwPython/blob/master/index.ipynb); the content is available [on GitHub](https://github.com/csaybar/EEwPython)* and released under the [Apache 2.0 License](https://www.gnu.org/licenses/gpl-3.0.en.html). 99% of this material has been adapted from [Google Earth Engine Guides](https://developers.google.com/earth-engine/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "g-iMfR17vqXj"
   },
   "source": [
    "<!--NAVIGATION-->\n",
    " < [Chart](7_Chart.ipynb) | [Contents](index.ipynb) |  [Specialized Algorithms](9_SpecializedAlgorithms.ipynb)>\n",
    "\n",
    "<a href=\"https://colab.research.google.com/github/csaybar/EEwPython/blob/master/8_Array.ipynb\"><img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open and Execute in Google Colaboratory\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "h-SbkLxfVRa7"
   },
   "source": [
    "<center>\n",
    "<h1>Google Earth Engine with Python </h1>\n",
    "<h2> Array </h2>\n",
    "</center>\n",
    "<h2> Topics:</h2>\n",
    "\n",
    "1. Array Overview\n",
    "2. Arrays and Array Images\n",
    "3. Array Transformations\n",
    "4. Eigen Analysis\n",
    "5. Array Sorting and Reducing\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "fpViOi9lnJMk"
   },
   "source": [
    "## Connecting GEE with Google Services"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "eq43ZIzVnLU6"
   },
   "source": [
    "- **Authenticate to Earth Engine**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "qvqKvZEnnM3B"
   },
   "outputs": [],
   "source": [
    "!pip install earthengine-api #earth-engine Python API"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "MkELcluCnNw6"
   },
   "outputs": [],
   "source": [
    "!earthengine authenticate "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "KdebAMvEnQL8"
   },
   "source": [
    "- **Authenticate to Google Drive (OPTIONAL)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "y1l8QKUinPAC"
   },
   "outputs": [],
   "source": [
    "from google.colab import drive\n",
    "drive.mount('/content/drive')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "A_9MwmnunShU"
   },
   "source": [
    "- **Authenticate to Google Cloud (OPTIONAL)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "MA146y_5nSKz"
   },
   "outputs": [],
   "source": [
    "from google.colab import auth\n",
    "auth.authenticate_user()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "hn_XrTB2nW-9"
   },
   "source": [
    "## Testing the software setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "YxzruSDonXud"
   },
   "outputs": [],
   "source": [
    "# Earth Engine Python API\n",
    "import ee\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "dlT4E9qPnY0t"
   },
   "outputs": [],
   "source": [
    "import folium\n",
    "\n",
    "# Define the URL format used for Earth Engine generated map tiles.\n",
    "EE_TILES = 'https://earthengine.googleapis.com/map/{mapid}/{{z}}/{{x}}/{{y}}?token={token}'\n",
    "\n",
    "print('Folium version: ' + folium.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "id": "hzaT4wLNne48"
   },
   "outputs": [],
   "source": [
    "#@title Mapdisplay: Display GEE Features or Images using folium.\n",
    "def Mapdisplay(center, dicc, Tiles=\"OpensTreetMap\",zoom_start=10):\n",
    "    '''\n",
    "    :param center: Center of the map (Latitude and Longitude).\n",
    "    :param dicc: Earth Engine Geometries or Tiles dictionary\n",
    "    :param Tiles: Mapbox Bright,Mapbox Control Room,Stamen Terrain,Stamen Toner,stamenwatercolor,cartodbpositron.\n",
    "    :zoom_start: Initial zoom level for the map.\n",
    "    :return: A folium.Map object.\n",
    "    '''\n",
    "    mapViz = folium.Map(location=center,tiles=Tiles, zoom_start=zoom_start)\n",
    "    for k,v in dicc.items():\n",
    "      if ee.image.Image in [type(x) for x in v.values()]:\n",
    "        folium.TileLayer(\n",
    "            tiles = v[\"tile_fetcher\"].url_format,\n",
    "            attr  = 'Google Earth Engine',\n",
    "            overlay =True,\n",
    "            name  = k\n",
    "          ).add_to(mapViz)\n",
    "      else:\n",
    "        folium.GeoJson(\n",
    "        data = v,\n",
    "        name = k\n",
    "          ).add_to(mapViz)\n",
    "    mapViz.add_child(folium.LayerControl())\n",
    "    return mapViz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 336
    },
    "colab_type": "code",
    "id": "B2nDCvpal1zn",
    "outputId": "83632f11-428e-4117-834d-be1351a09ab0"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<center><iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/-qo8L5GmKO0\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe></center>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 8,
     "metadata": {
      "tags": []
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#@title # Array Overview\n",
    "from IPython.display import HTML\n",
    "HTML('<center><iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/-qo8L5GmKO0\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe></center>')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VxW_4vQjg6tr"
   },
   "source": [
    "Earth Engine represents 1-D vectors, 2-D matrices, 3-D cubes, and higher dimensional hypercubes with the **ee.Array** type. Arrays are a flexible data structure, but in exchange for the power they offer, they do not scale as well as other data structures in Earth Engine. If the problem can be solved without using arrays, the result will be computed faster and more efficiently. But if the problem requires a higher dimension model, flexible linear algebra, or anything else arrays are uniquely suited to, you can use the **Array** class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "wmgIgnMYnEVE"
   },
   "source": [
    "## Array dimension, shape and size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DWiNfoESnH5L"
   },
   "source": [
    "The dimension of an array refers to the number of axes along which the underlying data varies. For example, 0-D arrays are scalar numbers, 1-D arrays are vectors, 2-D arrays are matrices, 3-D arrays are cubes, and >3-D arrays are hyper-cubes. For an N-dimensional array, there are N axes from 0 to N-1. The shape of the array is determined by the lengths of the axes. The length of an axis is the number of positions along it. The array size, or number of total elements in the array, equals the product of the axis lengths. Each value at every position on every axis must have a valid number, since sparse or ragged arrays are not currently supported. The array’s element type indicates what kind of number each element is; all elements of the array will have the same type."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HWulkd7Cn43S"
   },
   "source": [
    "# Arrays and Array Images\n",
    "\n",
    "Arrays in Earth Engine are constructed from lists of numbers and lists of lists. The degree of nesting determines the number of dimensions. To get started with a simple, motivated example, consider the following example of an Array created from Landsat 5 tasseled cap (TC) coefficients ([Crist and Cicone 1984](http://dx.doi.org/10.1109/TGRS.1984.350619))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "_nYPIgINoFZh"
   },
   "outputs": [],
   "source": [
    "# Create an Array of Tasseled Cap coefficients.\n",
    "coefficients = ee.Array([\n",
    "  [0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863],\n",
    "  [-0.2848, -0.2435, -0.5436, 0.7243, 0.0840, -0.1800],\n",
    "  [0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572],\n",
    "  [-0.8242, 0.0849, 0.4392, -0.0580, 0.2012, -0.2768],\n",
    "  [-0.3280, 0.0549, 0.1075, 0.1855, -0.4357, 0.8085],\n",
    "  [0.1084, -0.9022, 0.4120, 0.0573, -0.0251, 0.0238]\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "gdX5q-qZoIyl"
   },
   "source": [
    "Confirm that this is a 6x6, 2-D Array using length(), which will return the lengths of each axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "a-D1ZAv0oKIG"
   },
   "outputs": [],
   "source": [
    "# Print the dimensions.\n",
    "print(coefficients.length())\n",
    "# [6,6]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "2G4iWF95oPZ_"
   },
   "source": [
    "The following table illustrates the arrangement of the matrix entries along the 0-axis and the 1-axis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DcVfuZ4coht5"
   },
   "source": [
    "<tbody><tr>\n",
    "        <th></th><th></th><th colspan=\"6\">1-axis -&gt;</th>\n",
    "      </tr>\n",
    "      <tr class=\"alt\">\n",
    "        <td></td><td></td><td><b>0</b></td><td><b>1</b></td><td><b>2</b></td><td><b>3</b></td><td><b>4</b></td><td><b>5</b></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "        <td></td><td><b>0</b></td><td>0.3037</td><td>0.2793</td><td>0.4743</td><td>0.5585</td><td>0.5082</td><td>0.1863</td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "        <td></td><td><b>1</b></td><td>-0.2848</td><td>-0.2435</td><td>-0.5436</td><td>0.7243</td><td>0.0840</td><td>-0.1800</td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "        <td><b>0-axis</b></td><td><b>2</b></td><td>0.1509</td><td>0.1973</td><td>0.3279</td><td>0.3406</td><td>-0.7112</td><td>-0.4572</td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "        <td></td><td><b>3</b></td><td>-0.8242</td><td>0.0849</td><td>0.4392</td><td>-0.0580</td><td>0.2012</td><td>-0.2768</td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "        <td></td><td><b>4</b></td><td>-0.3280</td><td>0.0549</td><td>0.1075</td><td>0.1855</td><td>-0.4357</td><td>0.8085</td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "        <td></td><td><b>5</b></td><td>0.1084</td><td>-0.9022</td><td>0.4120</td><td>0.0573</td><td>-0.0251</td><td>0.0238</td>\n",
    "      </tr>\n",
    "    </tbody>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "0pjTiyZho-tR"
   },
   "source": [
    "The indices on the left of the table indicate positions along the 0-axis. The n-th element within each list on the 0-axis is in the n-th position along the 1-axis. For example, the entry at coordinate [3,1] of the array is 0.0849. Suppose ‘greenness’ is the TC component of interest. You can get the greenness sub-matrix using slice()."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "qiMOrK1johUS"
   },
   "outputs": [],
   "source": [
    "# Get the 1x6 greenness slice, display it.\n",
    "greenness = coefficients.slice(axis=0, start=1, end=2, step=1)\n",
    "print(greenness)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "usbZfV6gpFUo"
   },
   "source": [
    "The 2-D greenness matrix should look something like."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "I2sKe6Rpqxrl"
   },
   "source": [
    "**[[-0.2848, -0.2435, -0.5436, 0.7243, 0.084, -0.18 ]]**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "S3sLQL9-pHdI"
   },
   "source": [
    "Observe that the **start** and end parameters of **slice()** correspond to the 0-axis indices displayed in the table (**start** is inclusive and end is exclusive).\n",
    "\n",
    "To get a greenness image, matrix multiply the bands of a Landsat 5 image by the greenness matrix. To do that, first convert the multi-band Landsat image into an “Array Image”, where each pixel is an **Array** of band values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "l2HSPwfYpQsj"
   },
   "outputs": [],
   "source": [
    "# Load a Landsat 5 image, select the bands of interest.\n",
    "image = ee.Image('LANDSAT/LT05/C01/T1_TOA/LT05_044034_20081011').select(['B1', 'B2', 'B3', 'B4', 'B5', 'B7'])\n",
    "\n",
    "# Make an Array Image, with a 1-D Array per pixel.\n",
    "arrayImage1D = image.toArray()\n",
    "\n",
    "# Make an Array Image with a 2-D Array per pixel, 6x1.\n",
    "arrayImage2D = arrayImage1D.toArray(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "o8MGhBBlpYTT"
   },
   "source": [
    "In this example, note that **toArray()** converts **image** to an array image in which each pixel is a 1-D vector, the entries of which correspond to the 6 values at the corresponding positions in the bands of **image**. An array image of 1-D vectors created in this manner has no concept of 2-D shape. To perform 2-D only operations such as matrix multiplication, convert it into a 2-D array per-pixel image with **toArrray(1)**. In each pixel of the 2-D array image, there is a 6x1 matrix of band values. To see this, consider the following toy example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "7_V9PDlbpg0M"
   },
   "outputs": [],
   "source": [
    "array1D = ee.Array([1, 2, 3])            # [1,2,3]\n",
    "array2D = ee.Array.cat([array1D], 1)     # [[1],[2],[3]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZEE3ShmZplqN"
   },
   "source": [
    "Observe that the **array1D** vector varies along the 0-axis. The **array2D** matrix does as well, but it’s got an extra dimension. Calling **toArray(1)** on the array image is like calling cat(**bandVector**, 1) on every pixel. Using the 2-D array image, left multiply by an image where each pixel contains a 2-D matrix of greenness coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "SHYSDygepk8o"
   },
   "outputs": [],
   "source": [
    "# Do a matrix multiplication: 1x6 times 6x1.\n",
    "# Cast the greenness Array to an Image prior to multiplication.\n",
    "greennessArrayImage = ee.Image(greenness).matrixMultiply(arrayImage2D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3UJrecaQpwTf"
   },
   "source": [
    "The result is a new array image where every pixel is the 1x1 matrix that results from matrix multiplying the 1x6 greenness matrix (left) and the 6x1 band matrix (right). For display purposes, convert to a regular, one-band image with **arrayGet()**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "V1mHDR4hpy6o"
   },
   "outputs": [],
   "source": [
    "# Get the result from the 1x1 array in each pixel of the 2-D array image.\n",
    "greennessImage = greennessArrayImage.arrayGet([0, 0])\n",
    "\n",
    "# Display the input imagery with the greenness result.\n",
    "dicc = {\n",
    "    'image': image.getMapId({'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.5}),\n",
    "    'greenness': greennessImage.getMapId({'min': -0.1, 'max': 0.1})\n",
    "}\n",
    "\n",
    "coords = [37.562, -122.3]\n",
    "Mapdisplay(coords, dicc, \"Stamen Terrain\", 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Y7eakxGfp7Nw"
   },
   "source": [
    "Here is a complete example, which uses the entire coefficients array to compute multiple tasseled cap components at once and display the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ENl3em5fp80p"
   },
   "outputs": [],
   "source": [
    "# Define an Array of Tasseled Cap coefficients.\n",
    "coefficients = ee.Array([\n",
    "  [0.3037, 0.2793, 0.4743, 0.5585, 0.5082, 0.1863],\n",
    "  [-0.2848, -0.2435, -0.5436, 0.7243, 0.0840, -0.1800],\n",
    "  [0.1509, 0.1973, 0.3279, 0.3406, -0.7112, -0.4572],\n",
    "  [-0.8242, 0.0849, 0.4392, -0.0580, 0.2012, -0.2768],\n",
    "  [-0.3280, 0.0549, 0.1075, 0.1855, -0.4357, 0.8085],\n",
    "  [0.1084, -0.9022, 0.4120, 0.0573, -0.0251, 0.0238]\n",
    "])\n",
    "\n",
    "# Load a Landsat 5 image, select the bands of interest.\n",
    "image = ee.Image('LANDSAT/LT05/C01/T1_TOA/LT05_044034_20081011').select(['B1', 'B2', 'B3', 'B4', 'B5', 'B7'])\n",
    "\n",
    "# Make an Array Image, with a 1-D Array per pixel.\n",
    "arrayImage1D = image.toArray()\n",
    "\n",
    "# Make an Array Image with a 2-D Array per pixel, 6x1.\n",
    "arrayImage2D = arrayImage1D.toArray(1)\n",
    "\n",
    "# Do a matrix multiplication: 6x6 times 6x1.\n",
    "# Get rid of the extra dimensions.\n",
    "componentsImage = ee.Image(coefficients).matrixMultiply(arrayImage2D).arrayProject([0]).arrayFlatten([['brightness', 'greenness', 'wetness', 'fourth', 'fifth', 'sixth']])\n",
    "\n",
    "# Display the first three bands of the result and the input imagery.\n",
    "vizParams = {\n",
    "  'bands': ['brightness', 'greenness', 'wetness'],\n",
    "  'min': -0.1, 'max': [0.5, 0.1, 0.1]\n",
    "}\n",
    "\n",
    "# Display the input imagery with the greenness result.\n",
    "dicc = {\n",
    "    'image': image.getMapId({'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 0.5}),\n",
    "    'components': componentsImage.getMapId(vizParams)\n",
    "}\n",
    "\n",
    "coords = [37.562, -122.3]\n",
    "Mapdisplay(coords, dicc, \"Stamen Terrain\", 10)\n",
    "\n",
    "# Tasseled cap components “brightness” (red), “greenness” (green), and “wetness” (blue)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rbSAdlV5qDLa"
   },
   "source": [
    "Note that when getting bands from an array image, first get rid of the extra dimensions with **project()**, then convert it back to a regular image with **arrayFlatten()**. The output should look something like."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kFXnRngiy5OO"
   },
   "source": [
    "# Array Transformations\n",
    "\n",
    "Earth Engine supports array transformations such as transpose, inverse and pseudo-inverse. As an example, consider an ordinary least squares (OLS) regression of a time series of images. In the following example, an image with bands for predictors and a response is converted to an array image, then “solved” to obtain least squares coefficients estimates three ways. First, assemble the image data and convert to arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "4z1MsXbzzr_Y"
   },
   "outputs": [],
   "source": [
    "# This function masks the input with a threshold on the simple cloud score.\n",
    "def cloudMask (img):\n",
    "  cloudscore = ee.Algorithms.Landsat.simpleCloudScore(img).select('cloud')\n",
    "  return img.updateMask(cloudscore.lt(50))\n",
    "\n",
    "# Load a Landsat 5 image collection.\n",
    "  # Filter to get only two years of data.\n",
    "  # Filter to get only imagery at a point of interest.\n",
    "  # Mask clouds by mapping the cloudMask function over the collection.\n",
    "  # Select NIR and red bands only.\n",
    "  # Sort the collection in chronological order.\n",
    "collection = ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA').filterDate('2008-04-01', '2010-04-01').filterBounds(ee.Geometry.Point(-122.2627, 37.8735)).map(cloudMask).select(['B4', 'B3']).sort('system:time_start', True)\n",
    "\n",
    "# This function computes the predictors and the response from the input.\n",
    "  # Compute time of the image in fractional years relative to the Epoch.\n",
    "  # Compute the season in radians, one cycle per year.\n",
    "  # Return an image of the predictors followed by the response.\n",
    "  # 0. constant\n",
    "  # 1. linear trend\n",
    "  # 2. seasonal\n",
    "  # 3. seasonal\n",
    "  # 4. response\n",
    "def makeVariables(image):\n",
    "  year = ee.Image(image.date().difference(ee.Date('1970-01-01'), 'year'))\n",
    "  season = year.multiply(2 * math.pi)\n",
    "  return image.select().addBands(ee.Image(1)).addBands(year.rename('t')).addBands(season.sin().rename('sin')).addBands(season.cos().rename('cos')).addBands(image.normalizedDifference().rename('NDVI')).toFloat()\n",
    "\n",
    "# Define the axes of variation in the collection array.\n",
    "imageAxis = 0\n",
    "bandAxis = 1\n",
    "\n",
    "# Convert the collection to an array.\n",
    "array = collection.map(makeVariables).toArray()\n",
    "\n",
    "# Check the length of the image axis (number of images).\n",
    "arrayLength = array.arrayLength(imageAxis)\n",
    "# Update the mask to ensure that the number of images is greater than or\n",
    "# equal to the number of predictors (the linear model is solveable).\n",
    "array = array.updateMask(arrayLength.gt(4))\n",
    "\n",
    "# Get slices of the array according to positions along the band axis.\n",
    "predictors = array.arraySlice(bandAxis, 0, 4)\n",
    "response = array.arraySlice(bandAxis, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "2NfXcUL9z2Ys"
   },
   "source": [
    "Note that **arraySlice()** returns all the images in the time series for the range of indices specified along the **bandAxis** (the 1-axis). At this point, matrix algebra can be used to solve for the OLS coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ZhqiOeFWz-yk"
   },
   "outputs": [],
   "source": [
    "# Compute coefficients the hard way.\n",
    "coefficients1 = predictors.arrayTranspose().matrixMultiply(predictors).matrixInverse().matrixMultiply(predictors.arrayTranspose()).matrixMultiply(response)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "vKA2En3v0EyV"
   },
   "source": [
    "Although this method works, it is inefficient and makes for difficult to read code. A better way is to use the **pseudoInverse()** method (**matrixPseudoInverse()** for an array image)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mxgLLAxT0Jhl"
   },
   "outputs": [],
   "source": [
    "# Compute coefficients the easy way.\n",
    "coefficients2 = predictors.matrixPseudoInverse().matrixMultiply(response)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3FnUu1-C0O6G"
   },
   "source": [
    "From a readability and computational efficiency perspective, the best way to get the OLS coefficients is **solve()** (**matrixSolve()** for an array image). The **solve()** function determines how to best solve the system from characteristics of the inputs, using the pseudo-inverse for overdetermined systems, the inverse for square matrices and special techniques for nearly singular matrices:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "NVoCBETe0cVo"
   },
   "outputs": [],
   "source": [
    "# Compute coefficients the easiest way.\n",
    "coefficients3 = predictors.matrixSolve(response)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "yysZ-BCN0e75"
   },
   "source": [
    "To get a multi-band image, project the array image into a lower dimensional space, then flatten it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "pa9DH3Js0gXg"
   },
   "outputs": [],
   "source": [
    "# Turn the results into a multi-band image.\n",
    "  # Get rid of the extra dimensions.\n",
    "coefficientsImage = coefficients3.arrayProject([0]).arrayFlatten([['constant', 'trend', 'sin', 'cos']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "CURYQgmh0nj0"
   },
   "source": [
    "Examine the outputs of the three methods and observe that the resultant matrix of coefficients is the same regardless of the solver. That **solve()** is flexible and efficient makes it a good choice for general purpose linear modeling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VE5g2rV_0re0"
   },
   "source": [
    "# Eigen Analysis\n",
    "\n",
    "The [principal components (PC) transform](https://en.wikipedia.org/wiki/Principal_component_analysis)  (also known as the Karhunen-Loeve transform) is a spectral rotation that takes spectrally correlated image data and outputs uncorrelated data. The PC transform accomplishes this by diagonalizing the input band correlation matrix through Eigen-analysis. To do this in Earth Engine, use a covariance reducer on an array image and the **eigen()** command on the resultant covariance array. Consider the following function for that purpose (which is part of a complete example):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "EjzS3J8h1Ago"
   },
   "outputs": [],
   "source": [
    "def getPrincipalComponents(centered, scale, region):\n",
    "  # Collapse the bands of the image into a 1D array per pixel.\n",
    "  arrays = centered.toArray()\n",
    "  # Compute the covariance of the bands within the region.\n",
    "  co= 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()).arrayProject([0]).arrayFlatten([getNewBandNames('sd')])\n",
    "\n",
    "  # Turn the PCs into a P-band image, normalized by SD.\n",
    "      # Throw out an an unneeded dimension, [[]] -> [].\n",
    "      # Make the one band array image a multi-band image, [] -> image.\n",
    "      # Normalize the PCs by their SDs.\n",
    "  return principalComponents.arrayProject([0]).arrayFlatten([getNewBandNames('pc')]).divide(sdImage)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rftv9bxO1YUq"
   },
   "source": [
    "The input to the function is a mean zero image, a scale and a region over which to perform the analysis. Note that the input imagery first needs to be converted to a 1-D array image and then reduced using **ee.Reducer.centeredCovariance()**. The array returned by this reduction is the symmetric variance-covariance matrix of the input. Use the **eigen()** command to get the eigenvalues and eigenvectors of the covariance matrix. The matrix returned by **eigen()** contains the eigenvalues in the 0-th position of the 1-axis. As shown in the above function, use **slice()** to separate the eigenvalues and the eigenvectors. Each element along the 0-axis of the eigenVectors matrix is an eigenvector. As in the tasseled cap (TC) example, perform the transformation by matrix multiplying the **arrayImage** by the eigenvectors. In this example, each eigenvector multiplication results in a PC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JvG0-28j1q7-"
   },
   "source": [
    "# Array Sorting and Reducing\n",
    "\n",
    "Array sorting is useful for obtaining custom quality mosaics which involve reducing a subset of image bands according to the values in a different band. The following example sorts by a cloud index, then gets the mean of the least cloudy subset of images in the collection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "k0syIZd1231R"
   },
   "outputs": [],
   "source": [
    "# Define an arbitrary region of interest as a point.\n",
    "roi = ee.Geometry.Point(-122.26032, 37.87187)\n",
    "\n",
    "# Use these bands.\n",
    "bandNames = ee.List(['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11'])\n",
    "\n",
    "# Load a Landsat 8 collection.\n",
    "  # Select the bands of interest to avoid taking up memory.\n",
    "  # Filter to get only imagery at a point of interest.\n",
    "  # Filter to get only six months of data.\n",
    "  # Mask clouds by mapping the cloudMask function over the collection.\n",
    "  # This will add a cloud score band called 'cloud' to every image.\n",
    "\n",
    "def simpleCloudImage(image):\n",
    "  return ee.Algorithms.Landsat.simpleCloudScore(image)\n",
    "\n",
    "collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA').select(bandNames).filterBounds(roi).filterDate('2014-06-01', '2014-12-31').map(ee.Algorithms.Landsat.simpleCloudScore)\n",
    "\n",
    "# Convert the collection to an array.\n",
    "array = collection.toArray()\n",
    "\n",
    "# Label of the axes.\n",
    "imageAxis = 0\n",
    "bandAxis = 1\n",
    "\n",
    "# Get the cloud slice and the bands of interest.\n",
    "bands = array.arraySlice(bandAxis, 0, bandNames.length())\n",
    "clouds = array.arraySlice(bandAxis, bandNames.length())\n",
    "\n",
    "# Sort by cloudiness.\n",
    "sorted = bands.arraySort(clouds)\n",
    "\n",
    "# Get the least cloudy images, 20% of the total.\n",
    "numImages = sorted.arrayLength(imageAxis).multiply(0.2).int()\n",
    "leastCloudy = sorted.arraySlice(imageAxis, 0, numImages)\n",
    "\n",
    "# Get the mean of the least cloudy images by reducing along the image axis.\n",
    "mean = leastCloudy.arrayReduce(\n",
    "  reducer= ee.Reducer.mean(),\n",
    "  axes= [imageAxis]\n",
    ")\n",
    "\n",
    "# Turn the reduced array image into a multi-band image for display.\n",
    "meanImage = mean.arrayProject([bandAxis]).arrayFlatten([bandNames])\n",
    "\n",
    "# Display the input imagery with the greenness result.\n",
    "dicc = {\n",
    "    'roi': roi.getInfo(),\n",
    "    'meanImage': meanImage.getMapId({'bands': ['B5', 'B4', 'B2'], 'min': 0, 'max': 0.5})\n",
    "}\n",
    "\n",
    "coords = [37.562, -122.3]\n",
    "Mapdisplay(coords, dicc, \"Stamen Terrain\", 9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "0d9hxtgj7DDu"
   },
   "source": [
    "As in the linear modeling example, separate the bands of interest from the sort index using **arraySlice()** along the band axis. Then sort the bands of interest by the cloud index using **arraySort()**. After the pixels have been sorted by ascending cloudiness, use **arraySlice()** along the **imageAxis** to get 20% of the least cloudy pixels. Lastly, apply **arrayReduce()** along the **imageAxis** with a mean reducer to get the mean of the least cloudy pixels. The final step converts the array image back to a multi-band image for display.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zJIDLJzUhK_0"
   },
   "source": [
    "<!--NAVIGATION-->\n",
    " < [Chart](7_Chart.ipynb) | [Contents](index.ipynb) |  [Specialized Algorithms](9_SpecializedAlgorithms.ipynb)>\n",
    "\n",
    "<a href=\"https://colab.research.google.com/github/csaybar/EEwPython/blob/master/8_Array.ipynb\"><img align=\"left\" src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open in Colab\" title=\"Open and Execute in Google Colaboratory\"></a>"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [
    "h-SbkLxfVRa7",
    "wmgIgnMYnEVE",
    "u789tt7jtFcI",
    "HWulkd7Cn43S",
    "kFXnRngiy5OO",
    "VE5g2rV_0re0",
    "JvG0-28j1q7-"
   ],
   "include_colab_link": true,
   "name": "ee_Array.ipynb",
   "provenance": [],
   "toc_visible": true,
   "version": "0.3.2"
  },
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}