Ch6_Ocean_Example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 6 - Example: Ocean Data \n",
    "###  Days with sea surface temperature above a threshold\n",
    "\n",
    "In this chapter we exemplify the use of Sea Surface Temperature (SST) data in the cloud. \n",
    "\n",
    "This example analyzes a time series from an area of the ocean or a point. If an area, it averages SST values into a single value. Then it analyze the time series to assess when SST is above a given threshold. This could be used to study marine heatwaves, or use a SST threshold relevant to a marine species of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings \n",
    "warnings.simplefilter('ignore') \n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt \n",
    "import hvplot.pandas # this library helps to make interactive plots\n",
    "import hvplot.xarray\n",
    "import fsspec # these libraries help reading cloud data\n",
    "import s3fs\n",
    "import dask\n",
    "from dask.distributed import performance_report, Client, progress"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# input parameters\n",
    "\n",
    "# select either a range of lat/lon or a point. \n",
    "# If a point, set both entries to the same value\n",
    "latr = [19, 20] # make sure lat1 < lat2 since no test is done below to simplify the code\n",
    "lonr = [-158, -157] # lon1 < lon2, range -180:180. resolution daily 1km!\n",
    "\n",
    "# time range. data range available: 2002-06-01 to 2020-01-20. [start with a short period]\n",
    "dater = ['2012-01-01','2016-12-31'] # dates on the format 'YYYY-MM-DD' as string"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## We are going to use the Multi-Scale Ultra High Resolution (MUR) Sea Surface Temperature (SST) data set\n",
    "### This dataset is stored the Amazon (AWS) Cloud. For more info and links to the data detail and examples, see: https://registry.opendata.aws/mur/\n",
    "\n",
    "This dataset is stored in `zarr` format, which is an optimized format for the large datasets and the cloud. It is not stored as one 'image' at a time or a gigantic netcdf file, but in 'chunks', so it is perfect for extracting time series.\n",
    "\n",
    "First, we open the dataset and explore it, but we are not downloading anything yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first determine the file name using, in the format:\n",
    "# the s3 bucket [mur-sst], and the region [us-west-2], and the folder if applicable [zarr-v1] \n",
    "file_location = 'https://mur-sst.s3.us-west-2.amazonaws.com/zarr-v1'\n",
    "\n",
    "ds_sst = xr.open_zarr(file_location,consolidated=True) # open a zarr file using xarray\n",
    "# it is similar to open_dataset but it only reads the metadata\n",
    "\n",
    "ds_sst # we can treat it as a dataset!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now that we know what the file contains, we select our data (region and time), operate on it if needed (if a region, average), and download only the selected data \n",
    "It takes a while given the high resolution of the data. So, be patient.... and if you're only testing, might want to choose a small region and a short time period first. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# decide if a point or a region was given.\n",
    "if (latr[0]==latr[1]) | (lonr[0]==lonr[1]): # if we give it only one point\n",
    "    sst = ds_sst['analysed_sst'].sel(time = slice(dater[0],dater[1]),\n",
    "                                            lat  = latr[0], \n",
    "                                            lon  = lonr[0]\n",
    "                                           ).load()\n",
    "else: # if we give it an area, it extract the area and average SST over the area and returns a time series of SST\n",
    "    sst = ds_sst['analysed_sst'].sel(time = slice(dater[0],dater[1]),\n",
    "                                            lat  = slice(latr[0], latr[1]), \n",
    "                                            lon  = slice(lonr[0], lonr[1])\n",
    "                                           ).mean(dim={'lat','lon'}, skipna=True, keep_attrs=True).load() # skip 'not a number' (NaN) values and keep attributes\n",
    "\n",
    "sst = sst-273.15 # transform units from Kelvin to  Celsius\n",
    "sst.attrs['units']='deg C' # update units in metadata\n",
    "sst.to_netcdf('data/sst_example.nc') # saving the data, incase we want to come back to analyze the same data, but don't want to acquire it again from the cloud.\n",
    "sst # take a peak"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "### *Execute the next cell only if your reading the data from a file - either no access to cloud, or not want to keep reading from it. Skip otherwise. (No problem if you executed it by mistake).*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sst = xr.open_dataset('data/sst_example.nc') \n",
    "sst.close()\n",
    "sst = sst.analysed_sst # select only one variable\n",
    "sst"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Let's plot the data using two different libraries.\n",
    "#### - `matplotlib` that we already learn.\n",
    "#### - `hovplot` is a more interactive library for web display. It provides you with the data details when you hover your cursor over the figure. Very nice for inspecting the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# matplotlib method #\n",
    "print('matplotlib') \n",
    "sst.plot() # this is all you need\n",
    "\n",
    "# all the stuff here to make it look nice. \n",
    "plt.ylabel('SST ($^\\circ$C)')\n",
    "plt.xlabel('Year')\n",
    "plt.title('Location: '+str(latr)+'$^\\circ$N, '+str(lonr)+'$^\\circ$W')\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.show()\n",
    "\n",
    "# hovplot method #\n",
    "print('hovplot')\n",
    "df = pd.DataFrame(data=sst.data, index=sst.time.data,columns=['SST (C)'])\n",
    "df.index.name = 'Date'\n",
    "df.hvplot(grid=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Now, let's analyze our data.\n",
    "#### First, the basics: climatology and anomalies. Also plotting using `hovplot`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the climatology\n",
    "sst_climatology = sst.groupby('time.dayofyear').mean('time',keep_attrs=True,skipna=False) # Group by day, all years. skipna ignore missing (NaN) values \n",
    "sst_climstd = sst.groupby('time.dayofyear').std('time',keep_attrs=True,skipna=False) # Calculate standard deviation. Keep data attributes.\n",
    "\n",
    "# creates a dataset with climatology and standard deviaton for easy plotting with hvplot\n",
    "ds = xr.Dataset({'clim':sst_climatology,'+Std':sst_climatology+sst_climstd,'-Std':sst_climatology-sst_climstd}) # add standard deviation time series +/-\n",
    "ds.hvplot(color=['k','grey','grey'], grid=True, title='SST Climatology') # plot the climatology (black, and the standard deviation in grey)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the anomalies\n",
    "sst_anomaly = sst.groupby('time.dayofyear')-sst_climatology \n",
    "sst_anomaly_monthly = sst_anomaly.resample(time='1MS', loffset='15D').mean(keep_attrs=True,skipna=False) # calculate monthly anomalies/smoothing\n",
    "\n",
    "# make a plot \n",
    "plt.plot(sst_anomaly.time,sst_anomaly)\n",
    "plt.plot(sst_anomaly_monthly.time,sst_anomaly_monthly, 'r')\n",
    "\n",
    "plt.grid()\n",
    "plt.ylabel('SSTa (C)')\n",
    "plt.title('SST Anomalies')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## We analyze the data further by dividing it by a threshold.\n",
    "\n",
    "- One way is to set a threshold that has some relevance.  For example, a thermal threshold for a marine species we are studying. \n",
    "\n",
    "- Another way is choosing the maximum value in the climatology (mean value + 1 standard deviation), which we can calculate or read by hovering our cursor over the climatology plot above.\n",
    "\n",
    "### Once the threshold is choosen, we identify when SST is over that threshold, and count how many days that occurred each year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we define a function that take a threshold value, and analyze and plot our data\n",
    "def SST_above(thr):\n",
    "    \n",
    "    fig, axs = plt.subplots(1,2,figsize=(16,4)) # creates a figure with two panels\n",
    "    \n",
    "    # first part - values above threshold - timeseries\n",
    "    plt.subplot(1,2,1) # plot on the first panel (last number)\n",
    "    plt.plot(sst.time,sst.data, lw=1)\n",
    "    a=sst>=thr # test when data is equal or greater than the threshold. a is a logical array (True/False values)\n",
    "    plt.plot(sst.time[a], sst.data[a],'.r', markersize=3) # plot only the values equal or above threshold\n",
    "    # all stuff here to make it look good\n",
    "    plt.ylabel('SST ($^\\circ$C)')\n",
    "    plt.xlabel('Year')\n",
    "    plt.title('Location: '+str(latr)+'$^\\circ$N, '+str(lonr)+'$^\\circ$W')\n",
    "    plt.grid(True, alpha=0.3)\n",
    "    \n",
    "\n",
    "    # second part - days per year above threshold\n",
    "    plt.subplot(1,2,2) # plot on the second panel\n",
    "    dts = sst[sst>=thr].time # select dates when SST is equal or greater than the threshold. note that this time is not a logical array, but the time values\n",
    "    hot_days = dts.groupby('time.year').count() # agregate by year, by counting  \n",
    "    plt.bar(hot_days.year, hot_days) # bar plot of days per year\n",
    "    plt.xlim(int(dater[0][:4]), int(dater[1][:4])+1) # make it nice\n",
    "    plt.ylabel('No. days above '+str(np.round(thr,1))+'C')\n",
    "    plt.grid(True, alpha=0.3)\n",
    "    plt.show() # display and finaiize this figure, so the next is not overwritten\n",
    "\n",
    "## the actual analysis: two examples ##\n",
    "\n",
    "### Maximum climatology threshold\n",
    "thr = ds['+Std'].max() # setting threshold as maximum climatological value: mean + 1 standard deviation\n",
    "print('Max climatological SST = ',np.round(thr,1),'C')\n",
    "SST_above(thr) # Call function we defined\n",
    "\n",
    "### A relevant threshold. \n",
    "# For example, for hawaii (the select region), 28C is a relevant threshold for coral bleaching (https://coralreefwatch.noaa.gov/product/5km/tutorial/crw08a_bleaching_threshold.php)\n",
    "thr = 28\n",
    "print('\\n\\nBiologically relevant SST = ',thr,'C')\n",
    "SST_above(thr) # Call function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "#### Now, a different analsys of anomalously warm SST days. \n",
    "## Marine heatwaves\n",
    "Defined as any period with SST anomalies above the threshold determined by the 90th percentile value of a given period - in this case our data time period."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first, calculate the threshold: 90th percentile\n",
    "thr = np.percentile(sst_anomaly, 90)\n",
    "\n",
    "fig, axs = plt.subplots(3,1,figsize=(16,16)) # make a figure of 3 vertical panels\n",
    "\n",
    "# same plot as in our function above, but this time we are plotting the anomalies.\n",
    "plt.subplot(3,1,1) \n",
    "plt.plot(sst_anomaly.time,sst_anomaly.data, lw=1)\n",
    "plt.axhline(y=0, c='k', zorder=0, alpha=0.5) # add a line to highlight the x axis \n",
    "a=sst_anomaly>=thr # select data above the threshold\n",
    "plt.plot(sst_anomaly.time[a], sst_anomaly.data[a],'.r', markersize=3)\n",
    "# all stuff here to make it look good\n",
    "plt.ylabel('SST anomalies ($^\\circ$C)')\n",
    "plt.xlabel('Year')\n",
    "plt.title('Location: '+str(latr)+'$^\\circ$N, '+str(lonr)+'$^\\circ$W')\n",
    "plt.grid(True, alpha=0.3)\n",
    "\n",
    "# Now plot on the original data (not anomalies)\n",
    "plt.subplot(3,1,2) # second panel\n",
    "plt.plot(sst.time,sst.data, lw=1)\n",
    "plt.plot(sst.time[a], sst.data[a],'.r', markersize=3) # plot only the values equal or above threshold\n",
    "# all stuff here to make it look good\n",
    "plt.ylabel('SST ($^\\circ$C)')\n",
    "plt.xlabel('Year')\n",
    "plt.title('Location: '+str(latr)+'$^\\circ$N, '+str(lonr)+'$^\\circ$W')\n",
    "plt.grid(True, alpha=0.3)\n",
    "\n",
    "# plot of marine heatwave days  per year\n",
    "dts = sst_anomaly[sst_anomaly>=thr].time\n",
    "mhw = dts.groupby('time.year').count()\n",
    "plt.subplot(3,1,3) # third panel\n",
    "plt.bar(mhw.year,mhw)\n",
    "plt.ylabel('No. days SSTa > '+str(np.round(thr,1))+'C')\n",
    "plt.grid(True, alpha=0.3)\n",
    "plt.show()\n",
    "\n",
    "mhw # print the numbers of days"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resources\n",
    "\n",
    "For the cloud and data in the cloud, see resources listed in Chapter 5.\n",
    "\n",
    "### Resources specifically for this chapter:\n",
    "\n",
    "- [MUR SST Data](https://registry.opendata.aws/mur/). SST data in the cloud, with references the official datta website, examples and other resources.\n",
    "\n",
    "- [Pangeo OSM2020 Tutorial](https://github.com/pangeo-gallery/osm2020tutorial). This is a very good tutorial for ocean application and cloud computing. Plenty of examples. Many of the commands here are from this tutorial.\n",
    "\n",
    "### About MHW\n",
    "\n",
    "- [Marine heatwaves](http://www.marineheatwaves.org/all-about-mhws.html). A good place to begin to get info about the subject.\n",
    "\n",
    "- [Marine heatwaves code](https://github.com/ecjoliver/marineHeatWaves). Marine heatwaves code from E. Oliver.\n",
    "\n",
    "### If you want to learn more:\n",
    "\n",
    "- [Methods for accessing a AWS bucket](https://docs.aws.amazon.com/AmazonS3/latest/userguide/access-bucket-intro.html). Bucket is the name of the cloud storage object. S3 stands for Amazon's Simple Storage Service.\n",
    "\n",
    "- [hvplot site](https://hvplot.holoviz.org/index.html). Plotting tool used here.\n",
    "\n",
    "- [zarr](https://zarr.readthedocs.io/en/stable/). Learn more about this big data storage format."
   ]
  },
  {
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}