GeoMining_EE_Hops.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "view-in-github"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/neriiacopo/GeoMining-EE-Hops/blob/main_v2/GeoMining_EE_Hops.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NEtimk1sx3X9"
      },
      "source": [
        "# Install and Import Libraries\n",
        "Some libraries are not pre-installed in Colab, export them to Drive to make the process faster for future uses.\n",
        "\n",
        "**On the first run**, run the following two cells"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "vpTTc4k9MXCh",
        "outputId": "c1ca430e-5133-4ac6-a886-f9d9d425d7af"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Mounted at /content/drive\n"
          ]
        }
      ],
      "source": [
        "import os, sys\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive')\n",
        "nb_path = '/content/notebooks'\n",
        "os.symlink('/content/drive/My Drive/Colab Notebooks', nb_path)\n",
        "sys.path.insert(0,nb_path)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0fkd9JRrg0Z_"
      },
      "outputs": [],
      "source": [
        "!pip install ghhops_server pyproj flask flask-ngrok"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9td4Xle6xymD"
      },
      "source": [
        "# Authenticate to Earth Engine"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "7NJ_b591gwNG"
      },
      "outputs": [],
      "source": [
        "from flask import Flask\n",
        "from flask_ngrok import run_with_ngrok\n",
        "import ghhops_server as hs\n",
        "import numpy as np\n",
        "import ee\n",
        "from pyproj import CRS\n",
        "from pyproj.aoi import AreaOfInterest\n",
        "from pyproj.database import query_utm_crs_info\n",
        "from pyproj import Transformer"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MfaEOk1Ly-eM"
      },
      "source": [
        "follow the link to access with your account to Earth Engine and copy and paste the authentication key"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "NDKTF8NGjjt7"
      },
      "outputs": [],
      "source": [
        "ee.Authenticate()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "W7KmRy4FyCTm"
      },
      "source": [
        "# Run the Flask app to connect HOPS to EE"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "IMhasSNDjwrG",
        "outputId": "54464296-2e35-4467-d450-cb9161d09821"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            " * Serving Flask app \"__main__\" (lazy loading)\n",
            " * Environment: production\n",
            "\u001b[31m   WARNING: This is a development server. Do not use it in a production deployment.\u001b[0m\n",
            "\u001b[2m   Use a production WSGI server instead.\u001b[0m\n",
            " * Debug mode: off\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "INFO:werkzeug: * Running on http://127.0.0.1:5000/ (Press CTRL+C to quit)\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            " * Running on http://adc0-34-125-38-152.ngrok.io\n",
            " * Traffic stats available on http://127.0.0.1:4040\n"
          ]
        }
      ],
      "source": [
        "# App Setup -----------------------------------------------------------------------------------------------------\n",
        "\n",
        "# register hops app as middleware\n",
        "app = Flask(__name__)\n",
        "run_with_ngrok(app)\n",
        "\n",
        "# register hops app as middleware\n",
        "hops: hs.HopsFlask = hs.Hops(app)\n",
        "\n",
        "# Initialize earth engine\n",
        "try:\n",
        "    ee.Initialize()\n",
        "except Exception as e:\n",
        "    ee.Authenticate(auth_mode='notebook')\n",
        "    ee.Initialize()\n",
        "\n",
        "\n",
        "# Global Functions ----------------------------------------------------------------------------------------------\n",
        "\n",
        "def castToFloat(o):\n",
        "    list = []\n",
        "    for i in o:\n",
        "        list.append(float(i))\n",
        "\n",
        "    return list\n",
        "\n",
        "def pts_bbox(pts):\n",
        "    pts_py = []\n",
        "    for p in pts:\n",
        "        pts_py.append([p.X, p.Y])\n",
        "\n",
        "    aoi = ee.Geometry.Polygon(\n",
        "        [[[pts_py[0][0], pts_py[0][1]],\n",
        "        [pts_py[1][0], pts_py[1][1]],\n",
        "        [pts_py[2][0], pts_py[2][1]],\n",
        "        [pts_py[3][0], pts_py[3][1]]]], None, False)\n",
        "\n",
        "    return aoi;\n",
        "\n",
        "def pts_polygon(pts):\n",
        "\n",
        "    coords = []\n",
        "    for p in pts:\n",
        "\n",
        "        coord = [p.X,p.Y]\n",
        "        coords.append(coord)\n",
        "\n",
        "    geom_polygon = ee.Geometry.Polygon([ee.Geometry.LinearRing(coords)])\n",
        "\n",
        "    return geom_polygon\n",
        "\n",
        "def pts_multipts(pts):\n",
        "\n",
        "    coords = []\n",
        "    for p in pts:\n",
        "\n",
        "        coord = [p.X,p.Y]\n",
        "        coords.append(coord)\n",
        "\n",
        "    geom_pts = ee.Geometry.MultiPoint(coords)\n",
        "\n",
        "    return geom_pts\n",
        "\n",
        "def img_scaleTrim(image, band, mode, proj, scale, pts):\n",
        "\n",
        "    # Create reducer for sampling mode\n",
        "    reducer = \"ee.Reducer.\"+ mode +\"()\"\n",
        "\n",
        "    # Resample the image to assign custom scale\n",
        "    imageResampled = image \\\n",
        "        .reduceResolution(reducer= eval(reducer), maxPixels = 5000) \\\n",
        "        .reproject(crs='EPSG:4326', scale=scale)\n",
        "\n",
        "    # Extract the bounding region from the locations\n",
        "    aoi = pts_bbox(pts)\n",
        "\n",
        "    # Sample the image within the region\n",
        "    rec_sample = imageResampled.sampleRectangle(region=aoi, defaultValue=0)\n",
        "\n",
        "    # Prepare the outputs\n",
        "    np_values = np.array(rec_sample.get(band).getInfo())\n",
        "    values = np_values.flatten().tolist()\n",
        "\n",
        "    # Change from pixel to node (vertices)\n",
        "    H = np_values.shape[0] -1\n",
        "    W = np_values.shape[1] -1\n",
        "\n",
        "    return values,W,H;\n",
        "\n",
        "# App Components -----------------------------------------------------------------------------------------------\n",
        "\n",
        "@hops.component(\n",
        "    \"/ee_image\",\n",
        "    inputs=[\n",
        "        hs.HopsString(\"layer\", \"layer\"),\n",
        "        hs.HopsString(\"bands\", \"bands\"),\n",
        "        hs.HopsString(\"mode\",\"mode\", \"resampling mode to be applied. Default is mean\", default=\"mean\"),\n",
        "        hs.HopsNumber(\"scale\", \"scale\"),\n",
        "        hs.HopsPoint(\"pts\",\"pts\",\"the bounding box representing the area of analaysis. Note, provide it in the following order: min.Lon(X), max.Lon(X), min.Lat(Y), max.Lat(Y) aka LeftBottom, RightBottom, RightTop, LeftTop\", hs.HopsParamAccess.LIST)],\n",
        "    outputs=[\n",
        "        hs.HopsNumber(\"values\"),\n",
        "        hs.HopsNumber(\"W\"),\n",
        "        hs.HopsNumber(\"H\")\n",
        "    ],\n",
        ")\n",
        "\n",
        "def ee_image(layer,bands,mode,scale,pts):\n",
        "\n",
        "    # Select image layer\n",
        "    image = ee.Image(layer)\\\n",
        "        .select(bands)\n",
        "\n",
        "    # Extract original projection\n",
        "    proj = image.projection()\n",
        "\n",
        "    return img_scaleTrim(image, bands, mode, proj, scale, pts);\n",
        "\n",
        "\n",
        "@hops.component(\n",
        "    \"/ee_filterDate\",\n",
        "    inputs=[\n",
        "        hs.HopsString(\"layer\", \"layer\"),\n",
        "        hs.HopsString(\"dates\", \"dates\", \"first and last day to download\", hs.HopsParamAccess.LIST)],\n",
        "    outputs=[\n",
        "        hs.HopsString(\"layers\"),\n",
        "        hs.HopsString(\"days\")\n",
        "    ],\n",
        ")\n",
        "\n",
        "def ee_filterDate(layer,dates):\n",
        "\n",
        "    # Filter image collection\n",
        "    images = ee.ImageCollection(layer) \\\n",
        "             .filter(ee.Filter.date(dates[0],dates[1]))\n",
        "\n",
        "    # Get metadata\n",
        "    meta = images.getInfo()\n",
        "    length = len(meta[\"features\"])\n",
        "    layers = []\n",
        "    days = []\n",
        "\n",
        "    for i in range(length):\n",
        "        layer = meta[\"features\"][i][\"id\"]\n",
        "        day = meta[\"features\"][i][\"properties\"][\"system:index\"]\n",
        "        layers.append(layer)\n",
        "        days.append(day)\n",
        "\n",
        "    return (layers, days);\n",
        "\n",
        "\n",
        "@hops.component(\n",
        "    \"/ee_nd\",\n",
        "    inputs=[       \n",
        "        hs.HopsString(\"layer\", \"layer\"),\n",
        "        hs.HopsString(\"band1\", \"band1\"),        \n",
        "        hs.HopsString(\"band2\", \"band2\"),\n",
        "        hs.HopsString(\"mode\",\"mode\", \"resampling mode to be applied. Default is mean\", default=\"mean\"),\n",
        "        hs.HopsNumber(\"scale\", \"scale\"),\n",
        "        hs.HopsPoint(\"pts\",\"pts\",\"the region to sample. can be bounding box or polygon depending on bbox boolean\", hs.HopsParamAccess.LIST)],\n",
        "    outputs=[\n",
        "        hs.HopsNumber(\"values\"),\n",
        "        hs.HopsNumber(\"W\"),\n",
        "        hs.HopsNumber(\"H\")\n",
        "    ],\n",
        ")\n",
        "\n",
        "def ee_ND(layer,band1,band2,mode,scale,pts):\n",
        "\n",
        "    # Select the two bands to subtract\n",
        "    image1 = ee.Image(layer).select(band1)\n",
        "    image2 = ee.Image(layer).select(band2)\n",
        "\n",
        "    # Compute Normalized Difference\n",
        "    imageND = image1.subtract(image2) \\\n",
        "              .divide(image1.add(image2)) \\\n",
        "              .rename(\"ND\")\n",
        "\n",
        "    # Extract original projection\n",
        "    proj = image1.projection()\n",
        "\n",
        "    return img_scaleTrim(imageND, \"ND\", mode, proj, scale, pts)\n",
        "\n",
        "@hops.component(\n",
        "    \"/ee_cumCost\",\n",
        "    inputs=[\n",
        "        hs.HopsString(\"layer\", \"layer\", \"the layer from which to calculate the cost to traverse\"),\n",
        "        hs.HopsString(\"cost\", \"cost\", \"the cost band\"),\n",
        "        hs.HopsPoint(\"sources\", \"sources\", \"the location to where to calculate the proximity\", hs.HopsParamAccess.LIST),\n",
        "        hs.HopsNumber(\"maxdistance\", \"maxd\"),\n",
        "        hs.HopsNumber(\"scale\", \"scale\"),\n",
        "        hs.HopsPoint(\"pts\",\"bbox\",\"the bounding box representing the area of analaysis. Note, provide it in the following order: min.Lon(X), max.Lon(X), min.Lat(Y), max.Lat(Y) aka LeftBottom, RightBottom, RightTop, LeftTop\", hs.HopsParamAccess.LIST)],\n",
        "    outputs=[\n",
        "        hs.HopsNumber(\"values\"),\n",
        "        hs.HopsNumber(\"W\"),\n",
        "        hs.HopsNumber(\"H\")\n",
        "    ],\n",
        ")\n",
        "\n",
        "def ee_cumCost(layer,cost,sources,maxd,scale,pts):\n",
        "\n",
        "    # Compute the centroid location\n",
        "    geom_source = pts_multipts(sources)\n",
        "\n",
        "    # Rasterize location into an image where the geometry is 1, everything else is 0\n",
        "    imageSources = ee.Image().toByte().paint(geom_source, 1)\n",
        "\n",
        "    # Mask the sources image with itself.\n",
        "    imageSources = imageSources.updateMask(imageSources)\n",
        "\n",
        "    # Select cost data\n",
        "    cost = ee.Image(layer)\\\n",
        "        .select(cost)\n",
        "\n",
        "    # Extract original projection\n",
        "    proj = cost.projection()\n",
        "\n",
        "    # Resample resolution\n",
        "    costResample = cost \\\n",
        "        .reduceResolution(**{\n",
        "        'reducer': ee.Reducer.mean(),\n",
        "        'maxPixels': 5000\n",
        "        }) \\\n",
        "        .reproject(**{\n",
        "        'crs': proj,\n",
        "        'scale': scale\n",
        "        })\n",
        "\n",
        "    # Compute the cumulative cost\n",
        "    cumulativeCost = costResample.cumulativeCost(source=imageSources, maxDistance=ee.Number.float(maxd)) \\\n",
        "        .rename(\"cumcost\")\n",
        "\n",
        "    return img_scaleTrim(ee.Image(cumulativeCost), \"cumcost\", \"mean\", proj, scale, pts)\n",
        "\n",
        "@hops.component(\n",
        "    \"/ee_cumCostExtra\",\n",
        "    inputs=[\n",
        "        hs.HopsString(\"layer\", \"layer\", \"the layer from which to calculate the cost to traverse\"),\n",
        "        hs.HopsString(\"cost\", \"cost\", \"the cost band\"),\n",
        "        hs.HopsString(\"remap\", \"remap\", \"pairs of x=y values to remap the image into minute/meter cost values\", hs.HopsParamAccess.LIST),\n",
        "        hs.HopsNumber(\"default\",\"default\",\"default\"),\n",
        "        hs.HopsPoint(\"pts\", \"paint\", \"points to paint on top of cost map\", hs.HopsParamAccess.LIST),\n",
        "        hs.HopsNumber(\"val\",\"val\",\"cost applied to the paint in minute/meter\"),\n",
        "        hs.HopsPoint(\"sources\", \"sources\", \"the location to where to calculate the proximity\", hs.HopsParamAccess.LIST),\n",
        "        hs.HopsNumber(\"maxdistance\", \"maxd\"),\n",
        "        hs.HopsNumber(\"scale\", \"scale\"),\n",
        "        hs.HopsPoint(\"bounding box\",\"bbox\",\"the bounding box representing the area of analaysis. Note, provide it in the following order: min.Lon(X), max.Lon(X), min.Lat(Y), max.Lat(Y) aka LeftBottom, RightBottom, RightTop, LeftTop\", hs.HopsParamAccess.LIST)],\n",
        "    outputs=[\n",
        "        hs.HopsNumber(\"values\"),\n",
        "        hs.HopsNumber(\"W\"),\n",
        "        hs.HopsNumber(\"H\")\n",
        "    ],\n",
        ")\n",
        "\n",
        "def ee_cumCostExtra(layer,cost,remap,default,paint,val,sources,maxd,scale,pts):\n",
        "\n",
        "    # Compute the centroid location\n",
        "    geom_source = pts_multipts(sources)\n",
        "\n",
        "    # Rasterize location into an image where the geometry is 1, everything else is 0\n",
        "    imageSources = ee.Image().toByte().paint(geom_source, 1)\n",
        "\n",
        "    # Mask the sources image with itself.\n",
        "    imageSources = imageSources.updateMask(imageSources)\n",
        "    \n",
        "    # Prepare values for remapping\n",
        "    fromVals = []\n",
        "    toVals = []\n",
        "\n",
        "    for m in remap:\n",
        "        fromVals.append(m.split(\"=\")[0])\n",
        "        toVals.append(m.split(\"=\")[1])\n",
        "\n",
        "    # Convert int values\n",
        "    fromVals = castToFloat(fromVals)\n",
        "    toVals = castToFloat(toVals)\n",
        "\n",
        "     # Create the cost layer\n",
        "    cost = ee.Image(layer).select(cost) \\\n",
        "        .remap(fromVals, toVals, defaultValue=default)\n",
        "\n",
        "    # Extract original projection\n",
        "    proj = cost.projection()\n",
        "\n",
        "    # Resample resolution\n",
        "    costResample = cost \\\n",
        "        .reduceResolution(**{\n",
        "        'reducer': ee.Reducer.mean(),\n",
        "        'maxPixels': 5000\n",
        "        }) \\\n",
        "        .reproject(**{\n",
        "        'crs': proj,\n",
        "        'scale': scale\n",
        "        })\n",
        "\n",
        "    # Paint cost\n",
        "    geom_pts = pts_multipts(paint)\n",
        "\n",
        "    costResampledPaint = costResample.paint(geom_pts, val).rename(\"cumcost\")\n",
        "\n",
        "    # Compute the cumulative cost\n",
        "    cumulativeCost = costResampledPaint.cumulativeCost(source=imageSources, maxDistance=ee.Number.float(maxd)) \\\n",
        "        .rename(\"cumcost\")\n",
        "\n",
        "    return img_scaleTrim(ee.Image(cumulativeCost), \"cumcost\", \"mean\", proj, scale, pts)\n",
        "\n",
        "\n",
        "\n",
        "@hops.component(\n",
        "    \"/reproject_UTM\",\n",
        "    name=\"reproject based on bounding box\",\n",
        "    description=\"reproject locations from 4326 to local UTM\",\n",
        "\n",
        "    inputs=[\n",
        "        hs.HopsPoint(\"points\",\"pts\",\"the projected points\", hs.HopsParamAccess.LIST),\n",
        "        hs.HopsBoolean(\"bool\",\"xyz?\",\"should the boundingbox be moved to the origin?\", default=False)\n",
        "    ],\n",
        "    outputs=[\n",
        "        hs.HopsString(\"points\",\"p\",\"the projected points\")\n",
        "    ]\n",
        ")\n",
        "\n",
        "\n",
        "def reproject_UTM(pts, bool):\n",
        "\n",
        "    # Compute the centroid location\n",
        "    xs = []\n",
        "    ys = []\n",
        "\n",
        "    for p in pts:\n",
        "        xs.append(p.X)\n",
        "        ys.append(p.Y)\n",
        "\n",
        "    p_mean = [np.asarray(xs).mean(),np.asarray(ys).mean()]\n",
        "\n",
        "    # Extract UTM of the centroid location\n",
        "    utm_crs_list = query_utm_crs_info(\n",
        "        datum_name=\"WGS 84\",\n",
        "        area_of_interest=AreaOfInterest(\n",
        "            south_lat_degree=p_mean[1],\n",
        "            west_lon_degree=p_mean[0],\n",
        "            north_lat_degree=p_mean[1],\n",
        "            east_lon_degree=p_mean[0]\n",
        "        ),\n",
        "    )\n",
        "\n",
        "    # Project from 4326 to local UTM\n",
        "    WSG84_crs = CRS.from_epsg(4326)\n",
        "    utm_crs = CRS.from_epsg(utm_crs_list[0].code)\n",
        "    transformer = Transformer.from_crs(WSG84_crs, utm_crs, always_xy=True)\n",
        "\n",
        "    pts_UTM = []\n",
        "\n",
        "    p_leftbottom = transformer.transform(pts[0].X,pts[0].Y)\n",
        "\n",
        "    for p in pts:\n",
        "\n",
        "        # provide first X then Y -- Lon then Lat\n",
        "        p_UTM = transformer.transform(p.X,p.Y)\n",
        "\n",
        "        # Move to origin using left bottom corner\n",
        "        if bool == True:\n",
        "            p_UTM = [p_UTM[0] - p_leftbottom[0], p_UTM[1] - p_leftbottom[1]]\n",
        "\n",
        "        pts_UTM.append(str(\"{\" + str(p_UTM[0]) +\",\" + str(p_UTM[1]) + \",0}\"))\n",
        "\n",
        "    return pts_UTM;\n",
        "\n",
        "\n",
        "# Run App ------------------------------------------------------------------------------------------------------\n",
        "\n",
        "if __name__ == \"__main__\":\n",
        "    app.run()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RwnQvFQVyLZN"
      },
      "source": [
        "# Copy the URL to Grasshopper to live connect\n",
        "that is reported after:\n",
        "* Running on ..........................ngrok.io"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "authorship_tag": "ABX9TyO6lqUVNzpXHLY49IcVEEQK",
      "collapsed_sections": [],
      "include_colab_link": true,
      "name": "GeoMining_EE_Hops.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "name": "python",
      "version": "3.8.0 (tags/v3.8.0:fa919fd, Oct 14 2019, 19:37:50) [MSC v.1916 64 bit (AMD64)]"
    },
    "vscode": {
      "interpreter": {
        "hash": "d5f8ca89cbc44b49c451fd70bb201e08944c93046f7633edd91372dbac4a6a85"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}