S2_download_and_preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ee\n",
    "from IPython.display import clear_output\n",
    "import numpy as np\n",
    "import rasterio\n",
    "from glob import glob\n",
    "import time\n",
    "import subprocess\n",
    "import json\n",
    "import pandas as pd\n",
    "import sys\n",
    "sys.path.append('../func')\n",
    "from wrapper_rasterio import export_tiff\n",
    "# from pathlib import Path\n",
    "# from urllib.request import urlretrieve\n",
    "# import zipfile\n",
    "\n",
    "# ee.Authenticate()\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# S2 Download"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define several functions to be applied on EE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def maskS2clouds(image):\n",
    "    \"\"\"\n",
    "    masking cloud for sentinel2 with the band \"QA60\".\n",
    "    \"\"\"\n",
    "    qa = image.select('QA60')\n",
    "    cloudBitMask = 1 << 10\n",
    "    cirrusBitMask = 1 << 11\n",
    "    mask = qa.bitwiseAnd(cloudBitMask).eq(0).And(qa.bitwiseAnd(cirrusBitMask).eq(0))\n",
    "    return image.updateMask(mask).select(\"B.*\").copyProperties(image, [\"system:time_start\"])\n",
    "\n",
    "def calc_NDVI(image):\n",
    "    NDVI = image.normalizedDifference(['B8', 'B4']).rename('NDVI')\n",
    "\n",
    "    return NDVI.copyProperties(image, [\"system:time_start\"])\n",
    "\n",
    "def calc_ECARR(image):    \n",
    "\n",
    "    ECARR = image.expression(\n",
    "        '0.0008 * (10000 * B8A / (B3 * B5)) ** 1.2403', {\n",
    "            'B3': image.select('B3'),\n",
    "            'B5': image.select('B5'),\n",
    "            'B8A': image.select('B8A')\n",
    "        }).rename('ECARR')\n",
    "\n",
    "\n",
    "#     # Divide does not work\n",
    "#     ECARR = image.divide(10000).expression(\n",
    "#         '0.0008 * (B8A / (B3 * B5)) ** 1.2403', {\n",
    "#             'B3': image.select('B3'),\n",
    "#             'B5': image.select('B5'),\n",
    "#             'B8A': image.select('B8A')\n",
    "#         }).rename('ECARR')\n",
    "    \n",
    "    return ECARR.copyProperties(image, [\"system:time_start\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load an AOI for downloading  \n",
    "Make sure the CRS is **EPSG:4326**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# with open('./data/aoi/Albany_sub.geojson','r') as f:\n",
    "#     aoi_json = json.load(f)\n",
    "\n",
    "with open('./data/aoi/Albany_150km_rec.geojson','r') as f:\n",
    "    aoi_json = json.load(f)\n",
    "\n",
    "# with open('./data/aoi/Esperance_150km.geojson','r') as f:\n",
    "#     aoi_json = json.load(f)\n",
    "\n",
    "coords = aoi_json['features'][0]['geometry']['coordinates'][0][0]\n",
    "aoi = ee.Geometry.Polygon(coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download an ECARR image and export to a bucket"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "epsg = 28350\n",
    "resolution = 20\n",
    "gs_bucket = 'takahata-dev'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'state': 'COMPLETED', 'description': 'Download S2 data', 'creation_timestamp_ms': 1639018590397, 'update_timestamp_ms': 1639019119329, 'start_timestamp_ms': 1639018616079, 'task_type': 'EXPORT_IMAGE', 'destination_uris': ['https://console.developers.google.com/storage/browser/takahata-dev/sentinel2/'], 'attempt': 1, 'id': 'PKVI3QIZ6ECV2DGL7TJ7KHPQ', 'name': 'projects/earthengine-legacy/operations/PKVI3QIZ6ECV2DGL7TJ7KHPQ'}\n"
     ]
    }
   ],
   "source": [
    "img_collection = (ee.ImageCollection(f'COPERNICUS/S2_SR')\n",
    "                  .select(['B3','B5','B8A','QA60'])\n",
    "                  .filterBounds(aoi)\n",
    "                  .filterDate(ee.Date('2019-12-01'), ee.Date('2019-12-31'))\n",
    "                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))\n",
    "#                   .filter(ee.Filter.lt('NODATA_PIXEL_PERCENTAGE',5))\n",
    "                  .map(maskS2clouds)\n",
    "                  .map(calc_ECARR))\n",
    "\n",
    "info = img_collection.getInfo()\n",
    "info_date = info['features'][0]['properties']['system:index']\n",
    "\n",
    "img = ee.Image(img_collection.median()) # Take the median for a month\n",
    "\n",
    "task = ee.batch.Export.image.toCloudStorage(image=img,  # an ee.Image object.\n",
    "                                            region=aoi,  # an ee.Geometry object.\n",
    "                                            description='Download S2 data',\n",
    "                                            bucket=gs_bucket,   # Must be a root directory. E.g. \"d1/d2\" does not work\n",
    "                                            fileNamePrefix=f'sentinel2/ECARR_Albany_{info_date}', # Target directory can be specified here\n",
    "                                            crs=f'EPSG:{epsg}',\n",
    "                                            scale=resolution, # meter per pixel\n",
    "                                            maxPixels=1e+10\n",
    "                                           )\n",
    "\n",
    "task.start()\n",
    "while True:\n",
    "    clear_output(wait=True)\n",
    "    status = task.status()\n",
    "    print(status)\n",
    "    if status['state'] == 'COMPLETED' or status['state'] == 'FAILED':\n",
    "        break\n",
    "    time.sleep(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'state': 'COMPLETED', 'description': 'Download S2 data', 'creation_timestamp_ms': 1639019120860, 'update_timestamp_ms': 1639021029567, 'start_timestamp_ms': 1639019133015, 'task_type': 'EXPORT_IMAGE', 'destination_uris': ['https://console.developers.google.com/storage/browser/takahata-dev/sentinel2/'], 'attempt': 1, 'id': 'SGELGXZ6BWL5FFZP4JUULAUG', 'name': 'projects/earthengine-legacy/operations/SGELGXZ6BWL5FFZP4JUULAUG'}\n"
     ]
    }
   ],
   "source": [
    "img_collection = (ee.ImageCollection(f'COPERNICUS/S2_SR')\n",
    "                  .select(['B3','B5','B8A','QA60'])\n",
    "                  .filterBounds(aoi)\n",
    "                  .filterDate(ee.Date('2020-12-01'), ee.Date('2020-12-31'))\n",
    "                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))\n",
    "#                   .filter(ee.Filter.lt('NODATA_PIXEL_PERCENTAGE',5))\n",
    "                  .map(maskS2clouds)\n",
    "                  .map(calc_ECARR))\n",
    "\n",
    "info = img_collection.getInfo()\n",
    "info_date = info['features'][0]['properties']['system:index']\n",
    "\n",
    "img = ee.Image(img_collection.median()) # Take the median for a month\n",
    "\n",
    "task = ee.batch.Export.image.toCloudStorage(image=img,  # an ee.Image object.\n",
    "                                            region=aoi,  # an ee.Geometry object.\n",
    "                                            description='Download S2 data',\n",
    "                                            bucket=gs_bucket,   # Must be a root directory. E.g. \"d1/d2\" does not work\n",
    "                                            fileNamePrefix=f'sentinel2/ECARR_Albany_{info_date}', # Target directory can be specified here\n",
    "                                            crs=f'EPSG:{epsg}',\n",
    "                                            scale=resolution, # meter per pixel\n",
    "                                            maxPixels=1e+10\n",
    "                                           )\n",
    "\n",
    "task.start()\n",
    "while True:\n",
    "    clear_output(wait=True)\n",
    "    status = task.status()\n",
    "    print(status)\n",
    "    if status['state'] == 'COMPLETED' or status['state'] == 'FAILED':\n",
    "        break\n",
    "    time.sleep(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Quick visualization (available only when the image is small)    \n",
    "# visualization = {\n",
    "#     \"region\": aoi,\n",
    "#     \"scale\": 20,\n",
    "#     \"crs\": 'EPSG:28350'\n",
    "# }\n",
    "# url = img.getDownloadURL(visualization)\n",
    "# urlretrieve(url,'data/image/img_tmp.zip')\n",
    "\n",
    "# with zipfile.ZipFile(\"data/image/img_tmp.zip\",\"r\") as zip_ref:\n",
    "#     zip_ref.extractall(f\"data/image/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download an NDVI image and export to a bucket"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'state': 'COMPLETED', 'description': 'Download S2 data', 'creation_timestamp_ms': 1639021032452, 'update_timestamp_ms': 1639021631926, 'start_timestamp_ms': 1639021050597, 'task_type': 'EXPORT_IMAGE', 'destination_uris': ['https://console.developers.google.com/storage/browser/takahata-dev/sentinel2/'], 'attempt': 1, 'id': '42SIQVXP646WCCZT6RA2256G', 'name': 'projects/earthengine-legacy/operations/42SIQVXP646WCCZT6RA2256G'}\n"
     ]
    }
   ],
   "source": [
    "img_collection = (ee.ImageCollection(f'COPERNICUS/S2_SR')\n",
    "                  .select(['B4','B8','QA60'])\n",
    "                  .filterBounds(aoi)\n",
    "                  .filterDate(ee.Date('2019-12-01'), ee.Date('2019-12-31'))\n",
    "                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))\n",
    "                  .map(maskS2clouds)\n",
    "                  .map(calc_NDVI))\n",
    "\n",
    "info = img_collection.getInfo()\n",
    "info_date = info['features'][0]['properties']['system:index']\n",
    "\n",
    "\n",
    "img = ee.Image(img_collection.median())\n",
    "task = ee.batch.Export.image.toCloudStorage(image=img,\n",
    "                                            region=aoi,\n",
    "                                            description='Download S2 data',\n",
    "                                            bucket=gs_bucket,\n",
    "                                            fileNamePrefix=f'sentinel2/NDVI_Albany_{info_date}',\n",
    "                                            crs=f'EPSG:{epsg}',\n",
    "                                            scale=resolution, # meter per pixel\n",
    "                                            maxPixels=1e+10)\n",
    "\n",
    "task.start()\n",
    "while True:\n",
    "    clear_output(wait=True)\n",
    "    status = task.status()\n",
    "    print(status)\n",
    "    if status['state'] == 'COMPLETED' or status['state'] == 'FAILED':\n",
    "        break\n",
    "    time.sleep(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'state': 'COMPLETED', 'description': 'Download S2 data', 'creation_timestamp_ms': 1639022809141, 'update_timestamp_ms': 1639023193969, 'start_timestamp_ms': 1639022829717, 'task_type': 'EXPORT_IMAGE', 'destination_uris': ['https://console.developers.google.com/storage/browser/takahata-dev/sentinel2/'], 'attempt': 1, 'id': '2A3WS6SUZ7W6EIBYYHROWUXU', 'name': 'projects/earthengine-legacy/operations/2A3WS6SUZ7W6EIBYYHROWUXU'}\n"
     ]
    }
   ],
   "source": [
    "img_collection = (ee.ImageCollection(f'COPERNICUS/S2_SR')\n",
    "                  .select(['B4','B8','QA60'])\n",
    "                  .filterBounds(aoi)\n",
    "                  .filterDate(ee.Date('2020-12-01'), ee.Date('2020-12-31'))\n",
    "                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))\n",
    "                  .map(maskS2clouds)\n",
    "                  .map(calc_NDVI))\n",
    "\n",
    "info = img_collection.getInfo()\n",
    "info_date = info['features'][0]['properties']['system:index']\n",
    "\n",
    "\n",
    "img = ee.Image(img_collection.median())\n",
    "task = ee.batch.Export.image.toCloudStorage(image=img,\n",
    "                                            region=aoi,\n",
    "                                            description='Download S2 data',\n",
    "                                            bucket=gs_bucket,\n",
    "                                            fileNamePrefix=f'sentinel2/NDVI_Albany_{info_date}',\n",
    "                                            crs=f'EPSG:{epsg}',\n",
    "                                            scale=resolution, # meter per pixel\n",
    "                                            maxPixels=1e+10)\n",
    "\n",
    "task.start()\n",
    "while True:\n",
    "    clear_output(wait=True)\n",
    "    status = task.status()\n",
    "    print(status)\n",
    "    if status['state'] == 'COMPLETED' or status['state'] == 'FAILED':\n",
    "        break\n",
    "    time.sleep(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After the export is completed, copy the files from the bucket to your local/VM machine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Copying gs://takahata-dev/sentinel2/ECARR_20191202T015619_20191202T020454_T50HNG.tif...\n",
      "Copying gs://takahata-dev/sentinel2/ECARR_20201202T021349_20201202T022052_T50HMG.tif...\n",
      "Copying gs://takahata-dev/sentinel2/ECARR_Esperance_20191202T015619_20191202T020454_T50HQJ.tif...\n",
      "Copying gs://takahata-dev/sentinel2/ECARR_Esperance_20201201T015621_20201201T015638_T50HQJ.tif...\n",
      "\\ [4 files][  3.4 GiB/  3.4 GiB]   33.6 MiB/s                                   \n",
      "==> NOTE: You are performing a sequence of gsutil operations that may\n",
      "run significantly faster if you instead use gsutil -m cp ... Please\n",
      "see the -m section under \"gsutil help options\" for further information\n",
      "about when gsutil -m can be advantageous.\n",
      "\n",
      "Copying gs://takahata-dev/sentinel2/NDVI_20191202T015619_20191202T020454_T50HNG.tif...\n",
      "Copying gs://takahata-dev/sentinel2/NDVI_20201202T021349_20201202T022052_T50HMG.tif...\n",
      "Copying gs://takahata-dev/sentinel2/NDVI_Esperance_20191202T015619_20191202T020454_T50HQJ.tif...\n",
      "/ [7 files][  4.7 GiB/  4.7 GiB]   49.6 MiB/s                                   \n",
      "Operation completed over 7 objects/4.7 GiB.                                      \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Copy the downloaded data to the local/VM machine\n",
    "loc = 'data/image/'\n",
    "cmd = f'gsutil -m cp gs://takahata-dev/sentinel2/*.tif {loc}'\n",
    "subprocess.call(cmd.split())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-processing for labels from the government data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the government data, each pixel only has a unique ID where we need to refer to the definition table to get forest type, owner, etc.  \n",
    "For our purpose, we convert it to a map where each pixel has a forest code so that we are able to know the forest type immediately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "label_gov = rasterio.open('data/aus_for20_indigenous/ind_for20_Albany_28350.tif').read(1)\n",
    "value_types = np.sort(np.unique(label_gov))\n",
    "value_types = value_types[(value_types > 0) & (value_types < 1642)]\n",
    "lookup = pd.read_csv('data/aus_for20_indigenous/def_table.csv')\n",
    "\n",
    "label_gov_forcode = np.zeros(label_gov.shape, dtype=np.uint16)\n",
    "for v in value_types:\n",
    "    code = lookup.loc[lookup.Value == v, 'FOR_CODE'].values[0]\n",
    "    label_gov_forcode[label_gov == v] = code\n",
    "\n",
    "export_tiff(label_gov_forcode,\n",
    "            target='data/label/ind_for20_Albany_28350_forcode.tif',\n",
    "            reference='data/aus_for20_indigenous/ind_for20_Albany_28350.tif',\n",
    "            nodata=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-processing for S2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image files and the labels need to have exactly the same width and height.  \n",
    "If not, they need to be adjusted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "src = rasterio.open('data/label/ind_for20_Albany_28350_forcode.tif')\n",
    "ulx, lry, lrx, uly = src.bounds\n",
    "\n",
    "loc_source = 'data/image'\n",
    "images = glob(f'{loc_source}/*Albany*.tif')\n",
    "for image in images:\n",
    "    product_name = image[len(loc_source)+1:-4]\n",
    "    loc_input = f'{loc_source}/{product_name}.tif'\n",
    "    loc_output = f'{loc_source}/{product_name}_tmp.tif'\n",
    "    cmd = f'gdal_translate -projwin {ulx} {uly} {lrx} {lry} {loc_input} {loc_output}'\n",
    "    subprocess.call(cmd.split())\n",
    "    cmd = f'rm {loc_input}'\n",
    "    subprocess.call(cmd.split())\n",
    "    cmd = f'mv {loc_output} {loc_input}'\n",
    "    subprocess.call(cmd.split())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}