AU Prediction (ML)/GEE_export.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "identified-advocate",
   "metadata": {},
   "outputs": [],
   "source": [
    "# install dependencies\n",
    "import sys\n",
    "!{sys.executable} -m pip install gdal shapely geopandas rasterio xarray \\\n",
    "    earthengine-api matplotlib contextily --upgrade >> /dev/null"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "excited-anchor",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import shapely\n",
    "from osgeo import gdal, ogr\n",
    "import urllib\n",
    "import shutil\n",
    "import glob\n",
    "import zipfile\n",
    "import ee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ambient-quarterly",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "from IPython.core.display import HTML\n",
    "import contextily as ctx\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "split-graham",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define Pandas display settings\n",
    "pd.set_option('display.max_rows', 5)\n",
    "pd.set_option('display.max_columns', None)\n",
    "pd.set_option('display.width', None)\n",
    "pd.set_option('display.max_colwidth', None)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tropical-miracle",
   "metadata": {},
   "source": [
    "## Define helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "former-dairy",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create worldfile to define image coordinates\n",
    "def worldfile_tofile(fname, area, dimensions):\n",
    "    name, ext = os.path.splitext(fname)\n",
    "    # use QGIS worldfile names convention \n",
    "    jext = ext[1] + ext[-1] + 'w'\n",
    "    fname = os.path.join(str(os.extsep).join([name,jext]))\n",
    "    with open(fname, 'w') as outfile:\n",
    "        xres = (area[2]-area[0])/dimensions[0]\n",
    "        yres = (area[1]-area[3])/dimensions[1]\n",
    "        coefficients = [xres, 0, 0, yres, area[0], area[3]]\n",
    "        print('\\n'.join(map(str, coefficients)), file=outfile)\n",
    "\n",
    "# download GEE URL and save to file\n",
    "def geeurl_tofile(GEEurl, fname):\n",
    "    with urllib.request.urlopen(GEEurl) as response, open(fname, 'wb') as outfile:\n",
    "        shutil.copyfileobj(response, outfile)\n",
    "\n",
    "def gee_preview_tofile(GEEimage, vis, dimensions, fname=None):\n",
    "    GEEurl = GEEimage\\\n",
    "        .visualize(**vis)\\\n",
    "        .getThumbURL({'dimensions':dimensions, 'format': 'jpg'})\n",
    "    #print (GEEurl)\n",
    "    if fname is not None:\n",
    "        geeurl_tofile(GEEurl, fname)\n",
    "        worldfile_tofile(fname, area, dimensions)\n",
    "    return {'url': GEEurl, 'width': dimensions[0], 'height': dimensions[1]}\n",
    "\n",
    "def split_rect(rect, n):\n",
    "    lats = np.linspace(rect[0], rect[2], n+1)\n",
    "    lons = np.linspace(rect[1], rect[3], n+1)\n",
    "    #print (lats, lons)\n",
    "    cells = []\n",
    "    for lt1, lt2 in zip(lats.ravel()[:-1], lats.ravel()[1:]):\n",
    "        for ll1, ll2 in zip(lons.ravel()[:-1], lons.ravel()[1:]):\n",
    "            cell = [lt1, ll1, lt2, ll2]\n",
    "            cells.append(cell)\n",
    "    return cells\n",
    "\n",
    "def zipsbands2image(files):\n",
    "    dss = []\n",
    "    # merge separate file areas\n",
    "    for fname in sorted(files):\n",
    "        #print ('fname', fname)\n",
    "        zip = zipfile.ZipFile(fname)\n",
    "        # merge separate file to dataset\n",
    "        ds = xr.Dataset()\n",
    "        for bandname in zip.namelist():\n",
    "            varname = bandname.split('.')[1]\n",
    "            da = xr.open_rasterio(f'/vsizip/{fname}/{bandname}').squeeze(drop=True)\n",
    "            ds[varname] = da\n",
    "            da.close()\n",
    "        dss.append(ds)\n",
    "    return xr.merge(dss)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "imposed-lincoln",
   "metadata": {},
   "source": [
    "## Define GEE functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "accepting-magnitude",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to mask clouds using the Sentinel-2 QA band.\n",
    "def GEEmaskS2clouds(image):\n",
    "    # Get the pixel QA band.\n",
    "    qa = image.select('QA60')\n",
    "\n",
    "    # Bits 10 and 11 are clouds and cirrus, respectively.\n",
    "    cloudBitMask = 1 << 10\n",
    "    cirrusBitMask = 1 << 11\n",
    "\n",
    "    # Both flags should be set to zero, indicating clear conditions.\n",
    "    cloudMask = qa.bitwiseAnd(cloudBitMask).eq(0)\n",
    "    cirrusMask = qa.bitwiseAnd(cirrusBitMask).eq(0)\n",
    "\n",
    "    # Return the masked and scaled data, without the QA bands.\n",
    "    return image\\\n",
    "        .updateMask(cloudMask)\\\n",
    "        .updateMask(cirrusMask)\\\n",
    "        .divide(10000)\\\n",
    "        .select(\"B.*\")\\\n",
    "        .copyProperties(image, [\"system:time_start\"])\n",
    "\n",
    "# Function to cloud mask from the pixel_qa band of Landsat 8 SR data.\n",
    "def GEEmaskL8sr(image):\n",
    "    # Get the pixel QA band.\n",
    "    qa = image.select('pixel_qa')\n",
    "    \n",
    "    # Bits 3 and 5 are cloud shadow and cloud, respectively.\n",
    "    cloudShadowBitMask = 1 << 3\n",
    "    cloudsBitMask = 1 << 5\n",
    "\n",
    "    # Both flags should be set to zero, indicating clear conditions.\n",
    "    cloudShadowMask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)\n",
    "    cloudsMask = qa.bitwiseAnd(cloudsBitMask).eq(0)\n",
    "\n",
    "    # Return the masked image, scaled to reflectance, without the QA bands.\n",
    "    return image\\\n",
    "        .updateMask(cloudShadowMask)\\\n",
    "        .updateMask(cloudsMask)\\\n",
    "        .divide(10000)\\\n",
    "        .select(\"B[0-9]*\")\\\n",
    "        .copyProperties(image, [\"system:time_start\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "opposite-herald",
   "metadata": {},
   "source": [
    "## Connect to GEE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "otherwise-performer",
   "metadata": {},
   "outputs": [],
   "source": [
    "service_account = 'console@gee-export-308512.iam.gserviceaccount.com'\n",
    "service_key = '/Users/mbg/gee-export.json'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "tight-armstrong",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for GDAL\n",
    "os.environ[\"GOOGLE_APPLICATION_CREDENTIALS\"] = service_key"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "simple-lucas",
   "metadata": {},
   "outputs": [],
   "source": [
    "credentials = ee.ServiceAccountCredentials(service_account, service_key)\n",
    "ee.Initialize(credentials)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sensitive-biology",
   "metadata": {},
   "source": [
    "## Define work area and scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "fleet-scottish",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'type': 'Polygon',\n",
       " 'coordinates': [[[116.789084, -9.031312],\n",
       "   [117.291524, -9.031312],\n",
       "   [117.291524, -8.632405],\n",
       "   [116.789084, -8.632405],\n",
       "   [116.789084, -9.031312]]]}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# work area, decimal degrees\n",
    "area = [116.789084,  -9.031312, 117.291524,  -8.632405]\n",
    "# processing scale, m\n",
    "scale = 30\n",
    "\n",
    "GEEarea = ee.Geometry.Rectangle(*area)\n",
    "GEEarea.getInfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "spanish-paris",
   "metadata": {},
   "source": [
    "## Prepare image bands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "induced-hazard",
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://developers.google.com/earth-engine/datasets/catalog/CSP_ERGo_1_0_Global_ALOS_topoDiversity\n",
    "GEEdiversity = ee.Image(\"CSP/ERGo/1_0/Global/ALOS_topoDiversity\").rename('topoDiversity')\n",
    "# https://developers.google.com/earth-engine/datasets/catalog/CSP_ERGo_1_0_Global_ALOS_mTPI\n",
    "GEEposition = ee.Image(\"CSP/ERGo/1_0/Global/ALOS_mTPI\").rename('mTPI')\n",
    "# https://developers.google.com/earth-engine/datasets/catalog/JAXA_ALOS_PALSAR_YEARLY_FNF\n",
    "# 1 - Forest, 2 - Non-Forest, 3 - Water\n",
    "GEEfnf = ee.ImageCollection(\"JAXA/ALOS/PALSAR/YEARLY/FNF\").mosaic().select('fnf').rename('FNF')\n",
    "# https://developers.google.com/earth-engine/datasets/catalog/JAXA_ALOS_AW3D30_V3_2\n",
    "GEEdsm = ee.ImageCollection(\"JAXA/ALOS/AW3D30/V3_2\").mosaic().select('DSM')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "located-disease",
   "metadata": {},
   "source": [
    "### Sentinel-2 composite image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "broken-potter",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collected Sentinel-2 SR images 177\n",
      "Collected Sentinel-2 SR bands ['S2SRB1', 'S2SRB2', 'S2SRB3', 'S2SRB4', 'S2SRB5', 'S2SRB6', 'S2SRB7', 'S2SRB8', 'S2SRB8A', 'S2SRB9', 'S2SRB11', 'S2SRB12']\n"
     ]
    }
   ],
   "source": [
    "# https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR\n",
    "GEEs2srcol = ee.ImageCollection('COPERNICUS/S2_SR')\\\n",
    "    .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))\\\n",
    "    .map(GEEmaskS2clouds)\\\n",
    "    .filterBounds(GEEarea)\n",
    "GEEs2sr = GEEs2srcol.median().rename(GEEs2srcol.first().bandNames().map(lambda name: ee.String('S2SR').cat(name)))\n",
    "print ('Collected Sentinel-2 SR images', GEEs2srcol.size().getInfo())\n",
    "print ('Collected Sentinel-2 SR bands', GEEs2sr.bandNames().getInfo())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "convertible-postage",
   "metadata": {},
   "source": [
    "### Landsat 8 composite image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "tested-christian",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collected Landsat-8 TOA images 239\n",
      "Collected Landsat-8 TOA bands ['LC8TOAB1', 'LC8TOAB2', 'LC8TOAB3', 'LC8TOAB4', 'LC8TOAB5', 'LC8TOAB6', 'LC8TOAB7', 'LC8TOAB8', 'LC8TOAB9', 'LC8TOAB10', 'LC8TOAB11']\n"
     ]
    }
   ],
   "source": [
    "# https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_TOA\n",
    "GEElc8toacol = ee.ImageCollection('LANDSAT/LC08/C01/T1')\\\n",
    "    .filterBounds(GEEarea)\n",
    "# https://developers.google.com/earth-engine/apidocs/ee-algorithms-landsat-simplecomposite\n",
    "GEElc8toacomposite = ee.Algorithms.Landsat.simpleComposite(\n",
    "        collection=GEElc8toacol,\n",
    "        cloudScoreRange=20,\n",
    "        maxDepth=40,\n",
    "        asFloat=True\n",
    "    )\n",
    "GEElc8toa = GEElc8toacomposite.rename(GEElc8toacomposite.bandNames().map(lambda name: ee.String('LC8TOA').cat(name)))\n",
    "print ('Collected Landsat-8 TOA images', GEElc8toacol.size().getInfo())\n",
    "print ('Collected Landsat-8 TOA bands', GEElc8toa.bandNames().getInfo())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "literary-extra",
   "metadata": {},
   "source": [
    "#### Additional code to fetch Landsat-8 SR product\n",
    "#### https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_SR\n",
    "GEElc8srcol = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')\\\n",
    "    .filter(ee.Filter.lt('CLOUD_COVER', 20))\\\n",
    "    .map(GEEmaskL8sr)\\\n",
    "    .filterBounds(GEEarea)\n",
    "\n",
    "GEElc8sr = GEElc8srcol.median().rename(GEElc8srcol.first().bandNames().map(lambda name: ee.String('LC8SR').cat(name)))\n",
    "\n",
    "print ('Collected Landsat-8 SR images', GEElc8srcol.size().getInfo())\n",
    "\n",
    "print ('Collected Landsat-8 SR bands', GEElc8sr.bandNames().getInfo())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "color-andrews",
   "metadata": {},
   "source": [
    "## Build image mosaic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "instrumental-property",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "collected bands ['S2SRB1', 'S2SRB2', 'S2SRB3', 'S2SRB4', 'S2SRB5', 'S2SRB6', 'S2SRB7', 'S2SRB8', 'S2SRB8A', 'S2SRB9', 'S2SRB11', 'S2SRB12', 'LC8TOAB1', 'LC8TOAB2', 'LC8TOAB3', 'LC8TOAB4', 'LC8TOAB5', 'LC8TOAB6', 'LC8TOAB7', 'LC8TOAB8', 'LC8TOAB9', 'LC8TOAB10', 'LC8TOAB11', 'topoDiversity', 'mTPI', 'FNF', 'DSM']\n"
     ]
    }
   ],
   "source": [
    "GEEimage = \\\n",
    "    GEEs2sr\\\n",
    "    .addBands(GEElc8toa)\\\n",
    "    .addBands(GEEdiversity)\\\n",
    "    .addBands(GEEposition)\\\n",
    "    .addBands(GEEfnf)\\\n",
    "    .addBands(GEEdsm)\\\n",
    "    .clip(GEEarea)\n",
    "print ('collected bands', GEEimage.bandNames().getInfo())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worse-membrane",
   "metadata": {},
   "source": [
    "### Make image preview and save on local filesystem\n",
    "Note: we need to check the image preview quality before the image usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "official-opposition",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define GEE Landsat 8 visualization parameters\n",
    "#LC8vis = {'bands':['LC8SRB4', 'LC8SRB3', 'LC8SRB2'], 'min':0.03, 'max':0.18, 'gamma':1.4}\n",
    "LC8vis = {'bands':['LC8TOAB4', 'LC8TOAB3', 'LC8TOAB2'], 'min':0.03, 'max':0.18, 'gamma':1.4}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "romance-launch",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/2c41e47556a09a79c62803a5f2d83408-3ca612fa6c2b9f1aec93ed08fe140f10:getPixels\" width=\"160\" height=\"160\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# show small and fast preview low scale\n",
    "Image(**gee_preview_tofile(GEEimage, LC8vis, (160, 160)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "saved-memorial",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/1fa72922ac66f72165e368485bd77412-9ab0f75ddcaf046943d8a5b57dd4c65d:getPixels\n",
      "CPU times: user 51.7 ms, sys: 14.3 ms, total: 66 ms\n",
      "Wall time: 2min 16s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "print (gee_preview_tofile(GEEimage, LC8vis, (800, 800), 'data/preview.800x800.jpg')['url'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "greenhouse-bailey",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/1dadef2a9cb0cf6ae6d3171cd36d39bb-6d74250258c0499d5ec4f1b095bb13d7:getPixels\n",
      "CPU times: user 77.9 ms, sys: 37.1 ms, total: 115 ms\n",
      "Wall time: 4min 14s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "print (gee_preview_tofile(GEEimage, LC8vis, (1600, 1600), 'data/preview.1600x1600.jpg')['url'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "confidential-pharmacy",
   "metadata": {},
   "source": [
    "### Extract raw image bands and save on local filesystem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "relative-seller",
   "metadata": {},
   "outputs": [],
   "source": [
    "# split the area for n*n cells to follow download limits\n",
    "n = 3\n",
    "for idx, cell in enumerate(split_rect(area, n)):\n",
    "    GEEurl = GEEimage.getDownloadURL({'filePerBand':True, 'scale': scale, 'region': ee.Geometry.Rectangle(*cell)})\n",
    "    print (f'downloading cell {idx+1} from {n*n}', GEEurl)\n",
    "    fname = f'data/image_{idx}.{scale}m.zip'\n",
    "    geeurl_tofile(GEEurl, fname)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "gentle-patent",
   "metadata": {},
   "source": [
    "### Convert downloaded zipped images chunks to compressed NetCDF file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "removed-library",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "ds = zipsbands2image(glob.glob(f'data/image_*.{scale}m.zip'))\n",
    "# fix some variables data\n",
    "ds.mTPI.values[ds.mTPI==ds.mTPI.nodatavals] = np.nan\n",
    "ds.topoDiversity.values[ds.topoDiversity == -np.inf] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "visible-reynolds",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check data variables\n",
    "for varname in ds.data_vars:\n",
    "    print (varname, float(ds[varname].min()), float(ds[varname].max()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dependent-variance",
   "metadata": {},
   "outputs": [],
   "source": [
    "encoding = {var: dict(zlib=True, complevel=6) for var in ds.data_vars}\n",
    "ds.to_netcdf(f'data/image.{scale}m.nc', encoding=encoding)\n",
    "print (ds.dims)"
   ]
  }
 ],
 "metadata": {
  "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}