DL_March_Vents_v1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vents Group - Version 1\n",
    "OOI Data Labs Workshop - March 2019\n",
    "\n",
    "A possible interactive would include timeseries of 3D-Temp and Earthquakes that allows students to see the larger context, while viewing video stills from 2 different scenes.\n",
    "\n",
    "Potetial timeframe: June 2018-January 2019 (or to present)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Notebook Setup\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Earthquake Data - First Attempt\n",
    "Lets try the automatic approach first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in Earthquake data\n",
    "eq_url = 'http://axial.ocean.washington.edu/hypo71.dat'\n",
    "data = pd.read_fwf(eq_url)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks pretty good, but we need to reformat several columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert time\n",
    "data['DateTime'] = pd.to_datetime(data['yyyymmdd'].astype(str) + ' ' + data['HHMMSSS.SS'].astype(str), format='%Y%m%d %H%M %S.%f')\n",
    "\n",
    "# Convert Lat and Lon\n",
    "def dm2dd(s):\n",
    "  degrees, minutes, = s.split()\n",
    "  dd = float(degrees) + float(minutes)/60\n",
    "  return dd\n",
    "\n",
    "data['Latitude'] = data['Lat(D M)'].apply(dm2dd)\n",
    "data['Longitude'] = -data['Lon(D M)'].apply(dm2dd) #Add negative for West\n",
    "\n",
    "# Split the MW and NWR columns\n",
    "new = data['MW NWR'].str.split(' ', n=1, expand=True)\n",
    "data['MW'] = new[0].astype(float)\n",
    "data['NWR'] = new[1].astype(float)\n",
    "data.drop(columns=['yyyymmdd','HHMMSSS.SS','Lat(D M)','Lon(D M)','MW NWR'], inplace=True)\n",
    "\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quickplot\n",
    "fig, (ax1,ax2) = plt.subplots(2,1, sharex=True, figsize=(10,6))\n",
    "data.plot(ax=ax1, x='DateTime',y='Depth',marker='.',linestyle='');\n",
    "data.plot(ax=ax2, x='DateTime',y='MW',marker='.',linestyle='');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Earthquake Data - Second Attempt\n",
    "This time, we'll explicitely extract just the columns we need.  However, we'll still need to do some conversions, so this doesn't save us much code, but it does save a few lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.read_fwf?\n",
    "pd.read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_fwf(eq_url, colspecs=[(0,19),(20,28),(29,38),(40,45),(47,52)])\n",
    "\n",
    "# Convert time\n",
    "data.index = pd.to_datetime(data['yyyymmdd HHMMSSS.SS'].astype(str), format='%Y%m%d %H%M %S.%f')\n",
    "data.index.name = 'DateTime'\n",
    "\n",
    "data['Latitude'] = data['Lat(D M)'].apply(dm2dd)\n",
    "data['Longitude'] = -data['Lon(D M)'].apply(dm2dd) #Add negative for West\n",
    "\n",
    "data.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Earthqake map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's plot the quakes\n",
    "plt.plot(data.Longitude,data.Latitude,marker='.',linestyle='');\n",
    "xlim = plt.xlim()\n",
    "ylim = plt.ylim()\n",
    "\n",
    "# Add the OOI Sites\n",
    "sites = pd.read_csv('https://github.com/seagrinch/data-team-python/raw/master/infrastructure/sites.csv')\n",
    "plt.plot(sites.longitude,sites.latitude,'d',markersize=8)\n",
    "\n",
    "# Reset the limits\n",
    "plt.xlim(xlim)\n",
    "plt.ylim(ylim);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hexbin map\n",
    "First, let's create a hexbin map the regular way, but we'll set the aspect ratio of the plot so it's close to a Mercator projection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plt.hexbin(data.Longitude,data.Latitude,gridsize=200,bins='log')\n",
    "plt.colorbar();\n",
    "plt.plot(sites.longitude,sites.latitude,'d',markersize=8)\n",
    "plt.xlim(-130.1,-129.9)\n",
    "plt.ylim(45.8,46.1);\n",
    "aspect_ratio = np.cos(np.mean(plt.ylim())*np.pi/180)\n",
    "ax.set_aspect(aspect_ratio);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also get a little fancier and also use cartopy so we get the proper map projection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!apt-get -qq install python-cartopy python3-cartopy\n",
    "import cartopy.crs as ccrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.axes(projection=ccrs.PlateCarree())\n",
    "ax.set_extent([-130.2,-129.8,45.8,46.1], crs=ccrs.PlateCarree())\n",
    "plt.hexbin(data.Longitude,data.Latitude,gridsize=200,bins='log',transform=ccrs.Geodetic())\n",
    "plt.colorbar();\n",
    "plt.plot(sites.longitude,sites.latitude,'d',markersize=8)\n",
    "gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,\n",
    "                  linewidth=2, color='gray', alpha=0.5, linestyle='--')\n",
    "gl.xlabels_top = False\n",
    "gl.ylabels_right = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Daily Earthquake Averages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_sub = data.loc['2018-6-1':'2019-6-1']\n",
    "\n",
    "fig,(ax1,ax2) = plt.subplots(2,1,figsize=(8,6),sharex=True,sharey=False)\n",
    "# data_sub['MW'].resample('1H').mean().plot(ax=ax1,label='Hourly Mean',marker='.',markersize=1,linestyle='');\n",
    "data_sub['MW'].resample('D').mean().plot(ax=ax1,label='Daily Mean');\n",
    "data_sub['MW'].resample('D').median().plot(ax=ax1,label='Daily Median');\n",
    "ax1.legend();\n",
    "\n",
    "data_sub['MW'].resample('D').count().plot(ax=ax2,label='Count');\n",
    "\n",
    "ax1.set_ylabel('Magnitude')\n",
    "ax2.set_ylabel('Count');\n",
    "ax1.set_title('Daily Earthquakes at Axial');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TMPSF Data\n",
    "Next will request and process data from the Thermistor array at Axial Seamount."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Request Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# More setup\n",
    "import requests\n",
    "import os\n",
    "import re\n",
    "import xarray as xr\n",
    "! pip install netcdf4==1.5.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def request_data(reference_designator,method,stream,start_date,end_date):\n",
    "  site = reference_designator[:8]\n",
    "  node = reference_designator[9:14]\n",
    "  instrument = reference_designator[15:]\n",
    "  \n",
    "  # Create the request URL\n",
    "  api_base_url = 'https://ooinet.oceanobservatories.org/api/m2m/12576/sensor/inv'\n",
    "  data_request_url ='/'.join((api_base_url,site,node,instrument,method,stream))\n",
    "\n",
    "  # All of the following are optional, but you should specify a date range\n",
    "  params = {\n",
    "    'beginDT':start_date,\n",
    "    'endDT':end_date,\n",
    "    'format':'application/netcdf',\n",
    "    'include_provenance':'true',\n",
    "    'include_annotations':'true'\n",
    "  }\n",
    "\n",
    "  # Make the data request\n",
    "  r = requests.get(data_request_url, params=params, auth=(API_USERNAME, API_TOKEN))\n",
    "  data = r.json()\n",
    "  \n",
    "  # Return just the THREDDS URL\n",
    "  return data['allURLs'][0] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "API_USERNAME = ''\n",
    "API_TOKEN = ''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next line is used to request a dataset from the OOI Data Portal.  \n",
    "\n",
    "**Note, this line only needs to be run once** and then I recommend commenting it out so you don't accidentlally rerun it again. When you run it, simply save the outputted URL in the next line of code, so you can then grab the data.\n",
    "\n",
    "Of course, if you wish to grab a different time range, including more recent data, you will need to rerun the request line.  Also, note that it will take a few minutes for new data requests to be processed, so you will not be able to contiue on with the notebook until the data files are ready.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data Request\n",
    "# request_data('RS03ASHS-MJ03B-07-TMPSFA301','streamed','tmpsf_sample',\n",
    "#              '2018-06-01T00:00:00.000Z','2019-06-01T00:00:00.000Z')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load and Process Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is the output URL from the request_data line\n",
    "data_url = 'https://opendap.oceanobservatories.org/thredds/catalog/ooi/sage@marine.rutgers.edu/20190515T024738-RS03ASHS-MJ03B-07-TMPSFA301-streamed-tmpsf_sample/catalog.html'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_data(url,bad_inst=''):\n",
    "  tds_url = 'https://opendap.oceanobservatories.org/thredds/dodsC'\n",
    "  datasets = requests.get(url).text\n",
    "  urls = re.findall(r'href=[\\'\"]?([^\\'\" >]+)', datasets)\n",
    "  x = re.findall(r'(ooi/.*?.nc)', datasets)\n",
    "  for i in x:\n",
    "    if i.endswith('.nc') == False:\n",
    "      x.remove(i)\n",
    "  for i in x:\n",
    "    try:\n",
    "      float(i[-4])\n",
    "    except:\n",
    "      x.remove(i)\n",
    "  datasets = [os.path.join(tds_url, i) for i in x]\n",
    "  \n",
    "  # Remove extraneous files if necessary\n",
    "  selected_datasets = []\n",
    "  for d in datasets:\n",
    "    if (bad_inst) and bad_inst in d:\n",
    "      pass\n",
    "    else:\n",
    "      selected_datasets.append(d)\n",
    "#   print(selected_datasets)\n",
    "  \n",
    "  # Load in dataset\n",
    "  ds = xr.open_mfdataset(selected_datasets)\n",
    "  ds = ds.swap_dims({'obs': 'time'}) # Swap the primary dimension\n",
    "  ds = ds.chunk({'time': 1000}) # Used for optimization\n",
    "  ds = ds.sortby('time') # Data from different deployments can overlap so we want to sort all data by time stamp.\n",
    "  return ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grab the data\n",
    "tmpsf_data = get_data(data_url)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, for plotting, we will use the 3 center thermistors (1=bottom, 14=mid, and 15=top), as defined in the [TMPSF DPS](https://oceanobservatories.org/instrument-class/tmpsf/).  You can easily plot others as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's make a plot\n",
    "fig,ax = plt.subplots(4,1,figsize=(8,8),sharex=True)\n",
    "\n",
    "# And the Earthquake dataset at the top\n",
    "data.loc['2018-6-1':'2019-6-1']['MW'].resample('D').count().plot(ax=ax[0]);\n",
    "ax[0].set_ylabel('Eq Count')\n",
    "\n",
    "# Now add some selected thermistors\n",
    "thermistors = [1,14,15]\n",
    "spa = 1 # Subplot axes\n",
    "for thermistor in thermistors:\n",
    "  # Due to a quirk, I had to add .to_dataframe() to the next line so both variables are DataFrames\n",
    "  tmpsf_data[\"temperature%02d\" % (thermistor)].resample(time='1D').mean().to_dataframe().plot(ax=ax[spa],legend=False)\n",
    "  ax[spa].set_ylabel(\"temperature%02d\" % (thermistor))\n",
    "  spa = spa+1\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Export the dataset\n",
    "The next few lines will concatenate the earthquake dataset along with a few selected thermistors so we can export the full (daily averaged) dataset as a CSV file for use in Excel or an interactive widget."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First will create the averages as xarrays\n",
    "eq1 = data.loc['2018-6-1':'2019-6-1']['MW'].resample('D').count()\n",
    "eq1.index.name = 'time'\n",
    "eq1.name = 'Eq Count'\n",
    "eq1 = eq1.to_xarray()\n",
    "\n",
    "eq2 = data.loc['2018-6-1':'2019-6-1']['MW'].resample('D').mean()\n",
    "eq2.index.name = 'time'\n",
    "eq2.name = 'Eq Magnitude'\n",
    "eq2 = eq2.to_xarray()\n",
    "\n",
    "a = tmpsf_data[\"temperature01\"].resample(time='1D').mean()\n",
    "b = tmpsf_data[\"temperature14\"].resample(time='1D').mean()\n",
    "c = tmpsf_data[\"temperature15\"].resample(time='1D').mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we'll merge the datasets and convert to a pandas dataframe\n",
    "# This could potentially be done the other way around, but I found this order easier.\n",
    "x = xr.merge([a,b,c,eq1,eq2]).to_dataframe()\n",
    "\n",
    "# Print the first few rows\n",
    "x.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Quickplot\n",
    "x.plot();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export to CSV\n",
    "x.to_csv('axial_data.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "name": "DL March Vents v1.ipynb",
   "provenance": [],
   "toc_visible": true,
   "version": "0.3.2"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}