notebooks/MODIS_L2P_SST_DataCube.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "# The following notebook is no longer supported, see [this related notebook](https://github.com/podaac/tutorials/blob/master/notebooks/meetings_workshops/workshop_osm_2022/ECCO_ssh_sst_corr.ipynb) highlighting Zarr Data Cubes\n",
    "\n",
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Before-you-start\" data-toc-modified-id=\"Before-you-start-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Before you start</a></span></li><li><span><a <li><span><a href=\"#Hands-off-workflow\" data-toc-modified-id=\"Hands-off-workflow-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Hands-off workflow</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis Ready Data (ARD) workflow for MODIS Aqua L2P SST collection\n",
    "\n",
    "This notebook demonstrates how to created a gridded \"Data Cube\", essentialy an ARD, from native Level 2P sea surface temperature (SST) data from the MODIS Aqua (https://doi.org/10.5067/GHMDA-2PJ19) collection or dataset.  It can also be applied to Terra L2P SST and other similar L2 satellite collections.\n",
    "\n",
    "\n",
    "## Before you start\n",
    "\n",
    "Before you beginning this tutorial, make sure you have an account in the Earthdata Login, which is required to access data from the NASA Earthdata system. Please visit https://urs.earthdata.nasa.gov to register for an Earthdata Login account. It is free to create and only takes a moment to set up.\n",
    "\n",
    "You will also need a netrc file containing your NASA Earthdata Login credentials in order to execute this notebook. A netrc file can be created manually within text editor and saved to your home directory. For additional information see: [Authentication for NASA Earthdata](https://nasa-openscapes.github.io/2021-Cloud-Hackathon/tutorials/04_NASA_Earthdata_Authentication.html#authentication-via-netrc-file)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from urllib import request, parse\n",
    "from http.cookiejar import CookieJar\n",
    "import getpass\n",
    "import netrc\n",
    "import requests\n",
    "import urllib\n",
    "import json\n",
    "import pprint\n",
    "import time\n",
    "import os\n",
    "from os import makedirs, path\n",
    "from os.path import isdir, basename\n",
    "from urllib.parse import urlencode\n",
    "from urllib.request import urlopen, urlretrieve\n",
    "from datetime import datetime, timedelta\n",
    "from json import dumps, loads\n",
    "import shutil\n",
    "from osgeo import gdal, gdalconst\n",
    "#from nco import Nco\n",
    "import glob\n",
    "#import shutil\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "from pathlib import Path\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Constants\n",
    "\n",
    "# the local download directory\n",
    "download_dir = \"./modis-datacube-output\"\n",
    "\n",
    "# URS, CMR, Harmony roots\n",
    "edl = \"urs.earthdata.nasa.gov\"\n",
    "cmr = \"cmr.earthdata.nasa.gov\"\n",
    "harmony_root = 'https://harmony.earthdata.nasa.gov'\n",
    "\n",
    "download_dir = os.path.abspath(download_dir) + os.path.sep\n",
    "Path(download_dir).mkdir(parents=True, exist_ok=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ARD workflow\n",
    "\n",
    "The idea of this workflow is to:\n",
    "1) Discover MODIS granules within a region of interest (ROI)\n",
    "\n",
    "2) Make spatial subsets of the dicovered granules\n",
    "\n",
    "3) Remap those subsets to a common grid\n",
    "\n",
    "4) Aggregate the maps to a \"Data Cube\": a co-registed set of images aggregated into a single file the cover the ROI\n",
    "\n",
    "This is done using a combination NASA Earthdata Search and Harmony services, and the python modules for the GDAL image transformation software, and the NCO toolkit. \n",
    "\n",
    "The first code block sets the input datasets, and time and space bounds "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-08-11T17:34:08.712987Z",
     "start_time": "2020-08-11T17:34:08.710091Z"
    },
    "scrolled": true,
    "tags": [
     "parameters"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2020-07-01T00:01:15Z,2020-07-13T00:01:15Z\n",
      "-163,15,-153,25\n"
     ]
    }
   ],
   "source": [
    "# GHRSST MODIS Aqua L2P SST v2019.0 collection concept-id (CMR)\n",
    "ccid = \"C1940473819-POCLOUD\"\n",
    "\n",
    "# set the time  bounds  for search. \n",
    "start_time = '2020-07-01T00:01:15Z'\n",
    "stop_time = '2020-07-13T00:01:15Z'\n",
    "temporal_coverage = start_time + ',' + stop_time\n",
    "print(temporal_coverage)\n",
    "\n",
    "# set the spatial bounds\n",
    "# mid Atlantic region\n",
    "#west = -35. ; south = -5. ; east = -25. ; north = 5.\n",
    "\n",
    "# Hawaiian Is.\n",
    "west = -163 ; south = 15 ; east = -153 ; north = 25\n",
    "\n",
    "# Southern California Bight\n",
    "#west = -118 ; south = 32.5 ; east = -117 ; north = 33.5\n",
    "\n",
    "spatial_coverage =  str(west) + ',' + str(south) + ',' + str(east) + ',' + str(north)\n",
    "print(spatial_coverage)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-25T12:03:56.780074Z",
     "start_time": "2020-07-25T12:03:56.777273Z"
    }
   },
   "source": [
    "## Perform the subsetting operation using the Earthdata Harmony service\n",
    "\n",
    "Harmony allows us to use a CMR style query to find and subset all matching granules for a **collection, datetime, and spatial bounding box**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-08-11T17:34:21.088917Z",
     "start_time": "2020-08-11T17:34:21.085442Z"
    }
   },
   "outputs": [],
   "source": [
    "try:  \n",
    "    harmonyConfig = {\n",
    "        'collection_id': ccid,     \n",
    "        'ogc-api-coverages_version': '1.0.0',\n",
    "        'variable': 'all',\n",
    "        'lat': '(' + str(south) + \":\" + str(north) + ')',\n",
    "        'lon': '(' + str(west) + \":\" + str(east) + ')',        \n",
    "        'start': start_time,\n",
    "        'stop': stop_time\n",
    "    }\n",
    "\n",
    "    # subset granule\n",
    "    async_url = harmony_root + '/{collection_id}/ogc-api-coverages/{ogc-api-coverages_version}/collections/{variable}/coverage/rangeset?subset=lat{lat}&subset=lon{lon}&subset=time(\"{start}\":\"{stop}\")'.format(**harmonyConfig)\n",
    "\n",
    "    print('Request URL', async_url)\n",
    "    async_response = request.urlopen(async_url)\n",
    "    async_results = async_response.read()\n",
    "    async_json = json.loads(async_results)\n",
    "    pprint.pprint(async_json)\n",
    "\n",
    "\n",
    "except urllib.error.HTTPError as e:\n",
    "    print(f\"    [{datetime.now()}] FAILURE: {f}\\n\\n{e}\\n{e.read()}\\n\")\n",
    "    raise e                     \n",
    "except Exception as e:\n",
    "    print(f\"    [{datetime.now()}] FAILURE: {f}\\n\\n{e}\\n\")\n",
    "    raise e\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Monitor the Harmony job"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jobConfig = {\n",
    "    'jobID': async_json['jobID']\n",
    "}\n",
    "\n",
    "job_url = harmony_root+'/jobs/{jobID}'.format(**jobConfig)\n",
    "print('Job URL', job_url)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "job_response = request.urlopen(job_url)\n",
    "job_results = job_response.read()\n",
    "job_json = json.loads(job_results)\n",
    "\n",
    "print('Job response:')\n",
    "print()\n",
    "pprint.pprint(job_json)\n",
    "while job_json['status'] == 'running' and job_json['progress'] < 100:\n",
    "    print('Job status is running. Progress is ', job_json['progress'], '%. Trying again.')\n",
    "    time.sleep(10)\n",
    "    loop_response = request.urlopen(job_url)\n",
    "    loop_results = loop_response.read()\n",
    "    job_json = json.loads(loop_results)\n",
    "    if job_json['status'] == 'running':\n",
    "        continue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    Download subsets to local computer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Download the data\n",
    "file_urls = []\n",
    "\n",
    "for job_result in job_json['links']:\n",
    "\n",
    "    download_url = job_result['href']\n",
    "    file_name = job_result['title']\n",
    "    if file_name == 'Job Status':\n",
    "        continue\n",
    "        \n",
    "    if file_name == 'STAC catalog':\n",
    "        continue\n",
    "        \n",
    "    print(\"downloading \" + file_name)\n",
    "    \n",
    "    #store output in a defined directory with a meaninful filename based on orginal name\n",
    "    out_file = download_dir + file_name\n",
    "    print(out_file)\n",
    "    file_urls.append(out_file)\n",
    "    with request.urlopen(download_url) as response, open(out_file, 'wb') as output:\n",
    "      shutil.copyfileobj(response, output) \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Perform resampling/reprojection on subsets using the GDAL module and gdal.Warp().** GDAL will only work on one variable (\"layer\") at a time and also strip out important CF metadata and coordinate variables. Therefore we will use **NCO tricks** to correct these artifacts here and in the next steps. Use the 'pynco' module for NCO python bindings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nco = Nco()\n",
    "\n",
    "# Keyword args for gdal.Warp():\n",
    "# Set the output size in decimal degrees close to 1 km native resolution.\n",
    "# Use 'bilinear' interpolation (another option is 'cubicspline')\n",
    "kwargs = {'format': 'netCDF', 'copyMetadata': True, \n",
    "           'outputBounds': [west, south, east, north], \n",
    "           'xRes': 0.01,\n",
    "           'yRes' : 0.01,\n",
    "           'dstSRS':'+proj=longlat +datum=WGS84 +no_defs',\n",
    "           'resampleAlg': 'bilinear',\n",
    "         }\n",
    "print(kwargs)\n",
    "\n",
    "nc_vars = ['sea_surface_temperature', 'quality_level']\n",
    "\n",
    "# Loop through subsetted files (use file_urls as the loop list) and warp into defined region from kwargs{}\n",
    "\n",
    "for i in range(len(file_urls)):\n",
    "    for j in range(len(nc_vars)): \n",
    "      variable = nc_vars[j]\n",
    "     \n",
    "      # input filename\n",
    "      src_filename = file_urls[i] \n",
    "      print(\"source filename: \", src_filename)\n",
    "        \n",
    "      # load the netCDF 'layer' like sea_surface_temperature (variable)\n",
    "      nc_file = 'NETCDF:' + src_filename + ':' + variable\n",
    "      print(nc_file)\n",
    "    \n",
    "      # try/catch for GDAL steps. Dont work on empty netCDF file from subsetting operation  \n",
    "      try:\n",
    "        src = gdal.Open(nc_file, gdalconst.GA_ReadOnly)\n",
    "\n",
    "        # set output filename\n",
    "        out_filename = download_dir + 'subset_reproject-' + variable +  '-' + basename(file_urls[i]) \n",
    "        ds = gdal.Warp(out_filename, src, **kwargs)\n",
    "        print(\"\")\n",
    "\n",
    "        del ds\n",
    "        del src     \n",
    "        \n",
    "        # add time dimenson to 'Band1' using NCO\n",
    "        nco.ncecat(input=out_filename, output='tmp.nc', options=['-v Band1 -u time'])\n",
    "        nco.ncks(input='tmp.nc', output=out_filename, options=['-v Band1'])\n",
    "\n",
    "        # use NCO to copy (append with -A ) the time variable to the output and make it a record dimension\n",
    "        nco.ncks(input=src_filename, output=out_filename, options=['-v time -A --mk_rec_dmn time'])\n",
    "\n",
    "      except Exception as e:\n",
    "        #Errors can happen when downloading subsetted files because some might \n",
    "        #not actually fall in the bounding box, and are returned 'empty'\n",
    "        print(f\"    [{datetime.now()}] FAILURE: \\n\")\n",
    "     \n",
    "      else:\n",
    "        print(f\"    [{datetime.now()}]  SUCCESS for: {out_filename}\")\n",
    "        print()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Add variable level CF metadata information that GDAL stripped out using NCO commands.** Use the 'pynco' module for NCO python bindings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nco = Nco()\n",
    "\n",
    "sstConfig = {\n",
    "            'scale_factor': 0.005,     \n",
    "            'add_offset': 273.15,\n",
    "            'valid_max': 10000,\n",
    "            'valid_min': -1000,\n",
    "            'long_name': 'sea surface temperature',\n",
    "            'standard_name': 'sea_surface_skin_temperature',\n",
    "            'coverage_content_type': 'physicalMeasurement'\n",
    "            }\n",
    "\n",
    "qualityConfig = {\n",
    "            'valid_max': 0,\n",
    "            'valid_min': 5, \n",
    "            'flag_values': '0b, 1b, 2b, 3b, 4b, 5b',\n",
    "            'flag_meanings': 'no_data bad_data worst_quality low_quality acceptable_quality best_quality',\n",
    "            'long_name': 'quality level of SST pixel',\n",
    "            'coverage_content_type': 'qualityInformation'\n",
    "            }\n",
    "            \n",
    "nc_vars = ['sea_surface_temperature', 'quality_level']\n",
    "os.chdir(download_dir)\n",
    "\n",
    "for i in range(len(nc_vars)): \n",
    "    variable = nc_vars[i]\n",
    "    reg_ex = \"subset*\" + variable + \"*\"\n",
    "    for file in glob.glob(reg_ex):\n",
    "        if variable == 'sea_surface_temperature':\n",
    "          print(\" -> Updating \"+ file)\n",
    "          nco.ncrename(input=file, output=file, options=['-v .Band1,sea_surface_temperature'])\n",
    "        \n",
    "\n",
    "          # update the SST variable attributes\n",
    "          nco.ncatted(input=file, output=file, options=['-a scale_factor,'+ variable + ',o,f,{scale_factor}'.format(**sstConfig)])\n",
    "          nco.ncatted(input=file, output=file, options=['-a add_offset,'+ variable + ',o,f,{add_offset}'.format(**sstConfig)])\n",
    "          nco.ncatted(input=file, output=file, options=['-a valid_min,'+ variable + ',o,s,{valid_min}'.format(**sstConfig)])\n",
    "          nco.ncatted(input=file, output=file, options=['-a valid_max,'+ variable + ',o,s,{valid_max}'.format(**sstConfig)])\n",
    "          nco.ncatted(input=file, output=file, options=['-a long_name,'+ variable + ',o,c,\"{long_name}\"'.format(**sstConfig)])\n",
    "          nco.ncatted(input=file, output=file, options=['-a standard_name,'+ variable + ',o,c,{standard_name}'.format(**sstConfig)])\n",
    "          nco.ncatted(input=file, output=file, options=['-a coverage_content_type,'+ variable + ',o,c,{coverage_content_type}'.format(**sstConfig)])\n",
    "            \n",
    "        elif variable == 'quality_level':\n",
    "          print(\" -> Updating \"+ file)\n",
    "        \n",
    "          nco.ncrename(input=file, output=file, options=['-v .Band1,' + variable])\n",
    "\n",
    "          # update the quality variable attributes\n",
    "          nco.ncatted(input=file, output=file, options=['-a valid_min,'+ variable + ',o,s,{valid_min}'.format(**qualityConfig)])\n",
    "          nco.ncatted(input=file, output=file, options=['-a valid_max,'+ variable + ',o,s,{valid_max}'.format(**qualityConfig)])\n",
    "          nco.ncatted(input=file, output=file, options=['-a long_name,'+ variable + ',o,c,\"{long_name}\"'.format(**qualityConfig)])\n",
    "          nco.ncatted(input=file, output=file, options=['-a flag_values,'+ variable + ',o,c,\"{flag_values}\"'.format(**qualityConfig)])\n",
    "          nco.ncatted(input=file, output=file, options=['-a flag_meanings,'+ variable + ',o,c,\"{flag_meanings}\"'.format(**qualityConfig)])\n",
    "          nco.ncatted(input=file, output=file, options=['-a coverage_content_type,'+ variable + ',o,c,{coverage_content_type}'.format(**qualityConfig)])\n",
    "     \n",
    "    \n",
    "print( \"-Done-\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create the MODIS SST Data Cube.** Copy variable(s) to SST target files and catenate all of them into a final output  netCDF file using NCO."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "nco = Nco()\n",
    "os.chdir(download_dir)\n",
    "\n",
    "\n",
    "# Loop through the SST files and add variable quality_level to SST file (target)\n",
    "reg_ex = \"subset_reproject-sea_surface_temperature*\"\n",
    "\n",
    "print( \"Copying quality_level variables to target sst files . . .\")\n",
    "for sst_file in glob.glob(reg_ex):    \n",
    "    quality_file = sst_file.replace(\"sea_surface_temperature\", \"quality_level\")  \n",
    "    nco.ncks(input=quality_file, output=sst_file, options=['-v quality_level -A'])\n",
    "    \n",
    "    \n",
    "# Create the data cube using NCO ncrcat command\n",
    "print(\". . . -Done- \\n\\nCreating MODIS SST data cube . . . \")\n",
    "nco.ncrcat(input=glob.glob(reg_ex), output='MODIS_SST.data-cube.nc')\n",
    "print(\". . . -Done- \\n\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the data with xarray and perform some plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-08-11T17:34:22.362167Z",
     "start_time": "2020-08-11T17:34:22.359386Z"
    }
   },
   "outputs": [],
   "source": [
    "dataCube = download_dir + 'MODIS_SST.data-cube.nc'\n",
    "xds = xr.open_dataset(dataCube)\n",
    "\n",
    "# create objects for subplots\n",
    "fig, axes = plt.subplots(ncols=2)\n",
    "\n",
    "# plot a time series at a specific location\n",
    "xds.sea_surface_temperature.isel(lat = 200, lon = 300).plot(ax=axes[0], marker = '.', linestyle = 'None')\n",
    "\n",
    "# a histogram of all points in region of interest\n",
    "xds.sea_surface_temperature.plot(ax=axes[1])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.draw()\n",
    "\n",
    "# filter the dataset using quality information (quality_level value 4 and 5 are best data)\n",
    "qc_dataset = xds.where((xds['sea_surface_temperature'] < 310) & (xds['quality_level'] >= 4))\n",
    "\n",
    "# re-plot the time series at a specific location and the regional histogram\n",
    "fig, axes = plt.subplots(ncols=2)\n",
    "qc_dataset.sea_surface_temperature.isel(lat = 200, lon = 300).plot(ax=axes[0], marker = '.', linestyle = 'None')\n",
    "qc_dataset.sea_surface_temperature.plot(ax=axes[1])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.draw()"
   ]
  },
  {
   "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.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "241.528px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}