Example_8 Do some science!.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from dask_jobqueue import SLURMCluster\n",
    "# cluster = SLURMCluster(scheduler_port=8789)\n",
    "# cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from distributed import LocalCluster\n",
    "# cluster = LocalCluster(n_workers=3,threads_per_worker=1)\n",
    "# cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style=\"border: 2px solid white;\">\n",
       "<tr>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3>Client</h3>\n",
       "<ul>\n",
       "  <li><b>Scheduler: </b>tcp://localhost:8789\n",
       "  <li><b>Dashboard: </b><a href='http://localhost:8787/status' target='_blank'>http://localhost:8787/status</a>\n",
       "</ul>\n",
       "</td>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3>Cluster</h3>\n",
       "<ul>\n",
       "  <li><b>Workers: </b>15</li>\n",
       "  <li><b>Cores: </b>30</li>\n",
       "  <li><b>Memory: </b>360.00 GB</li>\n",
       "</ul>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Client: scheduler='tcp://146.118.38.43:8789' processes=15 cores=30>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from distributed import Client\n",
    "# client=Client(cluster)\n",
    "client=Client('tcp://localhost:8789')\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Against the Zarr on local (Lustre - fast!) file system "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how each time the number of lines of code is reducing, better chunking in the zarr probably helps in this. Most of the time is now in building the dask task graph and copying the dataset back to the memory of the notebook! The calculation only takes < 5 s (!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:    (latitude: 7001, longitude: 10001, time: 4882)\n",
       "Coordinates:\n",
       "  * latitude   (latitude) float64 10.0 9.99 9.98 9.97 ... -59.98 -59.99 -60.0\n",
       "  * longitude  (longitude) float64 80.0 80.01 80.02 80.03 ... 180.0 180.0 180.0\n",
       "  * time       (time) datetime64[ns] 2002-07-06 2002-07-07 ... 2016-12-31\n",
       "Data variables:\n",
       "    K_490      (time, latitude, longitude) float32 dask.array<shape=(4882, 7001, 10001), chunksize=(7, 856, 1223)>\n",
       "Attributes:\n",
       "    Conventions:  CF-1.6\n",
       "    history:      File initialised at 2015-12-16T15:54:44.125153\\nInitialised...\n",
       "    source_path:  imos-srs/archive/oc/aqua/1d/v201508/2002/07/A20020706.L2OC_..."
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import os\n",
    "zarr_path = os.environ['MYSCRATCH'] + '/aqua.P1D.aust.K_490.Full.zarr'\n",
    "ds_zarr = xr.open_zarr(zarr_path) # Already chunked on filesystem\n",
    "ds_zarr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14 years of daily analysis ready satellite data from IMOS-SRS team with complete australian coverage\n",
      "Dataset size 341.8 billion array elements, 1367.29 GB\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "numel = np.prod(ds_zarr.K_490.shape)\n",
    "print('14 years of daily analysis ready satellite data from IMOS-SRS team with complete australian coverage')\n",
    "print('Dataset size %.1f billion array elements, %.2f GB' % (numel/1E9,ds_zarr.nbytes/1E9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate anomalies from a climatology for a region of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "cropped = ds_zarr.sel(longitude=slice(113,117),latitude=slice(-20,-23))\n",
    "\n",
    "month_group = cropped.groupby('time.month')\n",
    "month_mean = month_group.mean('time')\n",
    "month_std = month_group.std('time')\n",
    "\n",
    "# persist these calculations to the cluster\n",
    "# kicks off the calculation and then retains the result spread across the workers for later use\n",
    "month_mean, month_std = client.persist([month_mean, month_std])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume that the light attenuation is gaussian distributed, define an anomoly as 2.5 times the standard deviation above the mean (possibly a poor assumption!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "month_group_2014 = cropped.sel(time=slice('2013-01-01', '2014-01-01')).groupby('time.month')\n",
    "month_mean_2014 = month_group_2014.mean('time')\n",
    "excess_2014 = month_mean_2014 - month_mean\n",
    "threshold = month_mean_2014 - (month_mean + 2.5*month_std)\n",
    "\n",
    "excess_2014, threshold = client.persist([excess_2014, threshold])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the results in the browser with holoviews+geoviews"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import holoviews as hv\n",
    "import hvplot.xarray\n",
    "import hvplot.pandas\n",
    "import geoviews as gv\n",
    "import numpy as np\n",
    "import cartopy\n",
    "from cartopy import crs as ccrs\n",
    "\n",
    "hv.notebook_extension('bokeh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "plate=ccrs.PlateCarree(central_longitude=133.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%opts Overlay [width=700 height=400]\n",
    "%%opts Image [projection=plate]\n",
    "\n",
    "thres_dataset = gv.Dataset(threshold, kdims=['month', 'longitude', 'latitude'], vdims=['K_490'])\n",
    "\n",
    "(excess_2014.K_490.hvplot('longitude','latitude',groupby='month',dynamic=True, rasterize=True, cmap='Oranges',crs=ccrs.PlateCarree()).redim(K_490=dict(range=(0, 0.2))))*\\\n",
    "(gv.feature.coastline.geoms('10m')).opts(color='black')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:c3dis2] *",
   "language": "python",
   "name": "conda-env-c3dis2-py"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}