tidal_hltc_plots.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tidal plots for HLTC National Map polygons\n",
    "\n",
    "**Summary**: Iterates through ITEMv2 polygons, extracts Landsat observations and produces plots of matching data\n",
    "\n",
    "**Issues**: \n",
    "\n",
    "**Notes**: Tidal model code based on original by Claire Phillips."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 7,
        "height": 5,
        "hidden": false,
        "row": 42,
        "width": 4
       },
       "report_default": {}
      }
     }
    }
   },
   "outputs": [],
   "source": [
    "# Import all modules\n",
    "import os\n",
    "import glob\n",
    "import datacube\n",
    "import fiona\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import datetime as dt\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from datacube.utils import geometry\n",
    "from datacube.api.query import query_group_by\n",
    "from otps import TimePoint, predict_tide\n",
    "from dateutil.relativedelta import relativedelta\n",
    "from matplotlib.ticker import FormatStrFormatter\n",
    "\n",
    "dc = datacube.Datacube(app='Tidal HLTC plots')\n",
    "\n",
    "\n",
    "def date_range(start_date, end_date, increment, period):\n",
    "    \n",
    "    \"\"\"Generate dates seperated by given time increment/period\"\"\"\n",
    "    \n",
    "    result = []\n",
    "    nxt = start_date\n",
    "    delta = relativedelta(**{period:increment})\n",
    "    while nxt <= end_date:\n",
    "        result.append(nxt)\n",
    "        nxt += delta\n",
    "    return result\n",
    "\n",
    "\n",
    "# Setup   \n",
    "filepath='/g/data/r78/intertidal/GA_native_tidal_model.shp'\n",
    "products = ['ls5_pq_albers', 'ls7_pq_albers', 'ls8_pq_albers'] \n",
    "time_period = ('1986-01-01', '2018-01-01')   # Global time range\n",
    "ls7_slc_period = ('1986-01-01', '2003-05-01') # Removes SLC failure\n",
    "\n",
    "# List of files \n",
    "hltc_files = glob.glob(\"/g/data/fk4/datacube/002/HLTC/HLTC_2_0/geotiff/COMPOSITE*.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "extensions": {
     "jupyter_dashboards": {
      "version": 1,
      "views": {
       "grid_default": {
        "col": 0,
        "height": 4,
        "hidden": false,
        "row": 71,
        "width": 4
       },
       "report_default": {}
      }
     }
    }
   },
   "source": [
    "## Produce plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing COMPOSITE_HIGH_171_115.11_-21.16_20000101_20170101_PER_20\n",
      "   Model tide heights\n",
      "   Observed tide heights\n"
     ]
    }
   ],
   "source": [
    "for hltc_file in hltc_files:\n",
    "    \n",
    "    try:\n",
    "\n",
    "        # Extract metadata from file\n",
    "        basename = os.path.basename(hltc_file)[0:-4]\n",
    "        print(\"Processing {}\".format(basename))\n",
    "        stage, polygon, lon, lat, time_from, time_to = basename.split(\"_\")[1:-2]\n",
    "\n",
    "        # Extract year, month, days\n",
    "        year_from, month_from, day_from = time_from[:4], time_from[4:6], time_from[6:]\n",
    "        year_to, month_to, day_to = time_to[:4], time_to[4:6], time_to[6:]\n",
    "\n",
    "        # Convert to datetime; take 1 second off \"to\" date to that the previous day is used for printing\n",
    "        dt_from = dt.datetime(int(year_from), int(month_from), int(day_from))\n",
    "        dt_to = dt.datetime(int(year_to), int(month_to), int(day_to)) - dt.timedelta(seconds = 1)\n",
    "\n",
    "        # Convert to string\n",
    "        string_from = \"{}-{}-{}\".format(dt_from.year, dt_from.month, dt_from.day)\n",
    "        string_to = \"{}-{}-{}\".format(dt_to.year, dt_to.month, dt_to.day)\n",
    "\n",
    "\n",
    "        ################################\n",
    "        print(\"   Model tide heights\") #\n",
    "        ################################\n",
    "\n",
    "        # For each hour between start and end of timeperiod, create list of datetimes\n",
    "        start = dt.datetime.strptime(time_period[0], \"%Y-%m-%d\")\n",
    "        end = dt.datetime.strptime(time_period[1], \"%Y-%m-%d\")\n",
    "        all_times = date_range(start, end, 1, 'hours')\n",
    "\n",
    "        # For each timestep in list, predict tide and add to list\n",
    "        times_model = [dt.strftime('%Y-%m-%d') for dt in all_times]\n",
    "        tp_model = [TimePoint(float(lon), float(lat), dt) for dt in all_times]\n",
    "        tides_model = [tide.tide_m for tide in predict_tide(tp_model)]\n",
    "\n",
    "        # Covert to dataframe of observed dates and tidal heights\n",
    "        df2_model = pd.DataFrame({'Model_height': tides_model}, index=pd.DatetimeIndex(times_model))\n",
    "\n",
    "\n",
    "        ###################################\n",
    "        print(\"   Observed tide heights\") #\n",
    "        ###################################\n",
    "\n",
    "        all_times_obs = list()\n",
    "\n",
    "        # Open ITEM polygons\n",
    "        with fiona.open(filepath) as Input:\n",
    "            crs = geometry.CRS(str(Input.crs_wkt))\n",
    "\n",
    "            # For each polygon, extract ID\n",
    "            for feature in Input:\n",
    "\n",
    "                Id = feature['properties']['ID']            \n",
    "\n",
    "                # If polygon is selected polygon\n",
    "                if Id == int(polygon):\n",
    "\n",
    "                    # Take first geometry and and extract area covered for input to dc call\n",
    "                    first_geometry = feature['geometry']\n",
    "                    geom = geometry.Geometry(first_geometry, crs=crs)\n",
    "\n",
    "                    # For each product:                                \n",
    "                    for source in products:   \n",
    "\n",
    "                        # Use entire time range unless LS7                \n",
    "                        time_range = ls7_slc_period if source == 'ls7_pq_albers' else time_period\n",
    "\n",
    "                        # Determine matching datasets for geom area and group into solar day\n",
    "                        ds = dc.find_datasets(product=source, time=time_range, geopolygon=geom)\n",
    "                        group_by = query_group_by(group_by='solar_day')\n",
    "                        sources = dc.group_datasets(ds, group_by)\n",
    "\n",
    "                        # If data is found, add time to list then sort\n",
    "                        if len(ds) > 0:\n",
    "                            all_times_obs.extend(sources.time.data.astype('M8[s]').astype('O').tolist()) \n",
    "\n",
    "                    # Discontinue loop if polygon is found\n",
    "                    break\n",
    "\n",
    "\n",
    "        # Calculate tide data from X-Y-time location \n",
    "        all_times_obs = sorted(all_times_obs)  \n",
    "        times_obs = [dt.strftime('%Y-%m-%d') for dt in all_times_obs]\n",
    "        tp_obs = [TimePoint(float(lon), float(lat), dt) for dt in all_times_obs]\n",
    "        tides_obs = [tide.tide_m for tide in predict_tide(tp_obs)]\n",
    "\n",
    "        # Covert to dataframe of observed dates and tidal heights\n",
    "        df1_obs = pd.DataFrame({'Tide_height': tides_obs}, index=pd.DatetimeIndex(times_obs)) \n",
    "\n",
    "\n",
    "        #############################\n",
    "        print(\"   Plotting output\") #\n",
    "        #############################\n",
    "\n",
    "        # Filter observed data by time\n",
    "        df1_obs_subset = df1_obs.loc[string_from:string_to]\n",
    "\n",
    "        # Compute percentage tide height\n",
    "        min_height=df1_obs_subset.Tide_height.min()\n",
    "        max_height=df1_obs_subset.Tide_height.max()\n",
    "        dr = max_height - min_height\n",
    "\n",
    "        # Create dict of percentile values\n",
    "        per10_dict = {perc:min_height + dr * perc * 0.01 for perc in range(0, 110, 10)}\n",
    "\n",
    "        # Filter observed data by tidal stage, using 0-20 for low composites and 80-100 for high\n",
    "        tide_from, tide_to = {'LOW': (0, 20), 'HIGH': (80, 100)}[stage]\n",
    "        df1_obs_subset = df1_obs_subset[df1_obs_subset.Tide_height.between(per10_dict[tide_from], \n",
    "                                                                           per10_dict[tide_to])]\n",
    "\n",
    "        # Plot setup\n",
    "        fig = plt.figure(figsize=(12,4))\n",
    "        plt.margins(0) \n",
    "        fig.axes[0].spines['right'].set_visible(False)\n",
    "        fig.axes[0].spines['top'].set_visible(False)\n",
    "        fig.axes[0].yaxis.set_major_formatter(FormatStrFormatter('%.2f'))\n",
    "        plt.ylabel('Tide height (m)')\n",
    "#         plt.title(\"{0}\\nTidal range between {1:g} and {2:g}m, period {3} / {4}\".format(basename, \n",
    "#                                                                                        per10_dict[tide_from],\n",
    "#                                                                                        per10_dict[tide_to],\n",
    "#                                                                                        string_from,\n",
    "#                                                                                        string_to))\n",
    "\n",
    "        # Plot observations and modelled values\n",
    "        plt.plot(df2_model.index, df2_model.Model_height, \n",
    "                 color='gainsboro', linewidth=0.5, zorder=1, label = 'Modelled')\n",
    "\n",
    "        # Plot output\n",
    "        plt.scatter(df1_obs.index, df1_obs.Tide_height, \n",
    "                    s=4, color='silver', marker='o', zorder=2)\n",
    "        plt.scatter(df1_obs_subset.index, df1_obs_subset.Tide_height, \n",
    "                    s=10, color='black', marker='o', zorder=3, label = 'Observed')\n",
    "\n",
    "        # Save and plot\n",
    "        plt.savefig('figures/hltc_plots/{}.jpg'.format(basename), \n",
    "                    bbox_inches=\"tight\", pad_inches=0.05, dpi = 150)\n",
    "        plt.close()\n",
    "    \n",
    "    except:\n",
    "        print(\"Skipping\")\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/g/data/fk4/datacube/002/HLTC/HLTC_2_0/geotiff/COMPOSITE_HIGH_209_122.2_-18.07_20000101_20170101_PER_20.tif'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hltc_files[50]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "extensions": {
   "jupyter_dashboards": {
    "activeView": "grid_default",
    "version": 1,
    "views": {
     "grid_default": {
      "cellMargin": 10,
      "defaultCellHeight": 20,
      "maxColumns": 12,
      "name": "grid",
      "type": "grid"
     },
     "report_default": {
      "name": "report",
      "type": "report"
     }
    }
   }
  },
  "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.3"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": false,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {
    "height": "40px",
    "left": "24px",
    "right": "1029px",
    "top": "139px",
    "width": "227px"
   },
   "toc_section_display": "none",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}