Adding an example for CERA workflow

examples/CERA_workflow.ipynb

@@ -0,0 +1,227 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "51b8f9b4-c76a-4a78-97fc-652a0ad8f9d7",
+   "metadata": {},
+   "source": [
+    "The CERA workflow for stations involved the following station filtering and adjustments:\n",
+    "- Get by parameters\n",
+    "- Get by region (coastal)\n",
+    "- Adjust by vertical datum\n",
+    "- Adjust time zone\n",
+    "- Station active vs inactive"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f8dfa70e-f62f-453d-b107-b0eadb6e1b40",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "from datetime import datetime, timedelta\n",
+    "from matplotlib import pyplot\n",
+    "\n",
+    "import geopandas as gpd\n",
+    "\n",
+    "from searvey import usgs, stations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0e336dca-652b-484c-b0c4-c986954d4553",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
+    "us = world[world.name.isin(['United States of America', 'Puerto Rico'])]\n",
+    "us_coast = us.boundary.intersection(world.unary_union.boundary)\n",
+    "ax = world.plot(color='k', alpha=0.1)\n",
+    "us.plot(ax=ax, color='b', alpha=0.2)\n",
+    "us_coast.plot(ax=ax, color='r')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1f676ca-3b2f-47fa-90b6-dfcae780202a",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "params_of_interest = ['62620', '62615']\n",
+    "region_of_interest = us_coast.unary_union.buffer(0.5)  # Buffer coast lines to overlap with some stations."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "58d786be-3560-499e-b787-a61436f2b4c0",
+   "metadata": {},
+   "source": [
+    "Note that currently USGS implemented all parameters of interest by CERA workflow, for further filtering one needs to fetch all and then filter. Also note that currently `stations.get_stations` API doesn't have paramter information."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ac8a7ad7-07e0-4e4f-b674-0e55bd70388a",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "#usgs_stations = stations.get_stations(providers='USGS', region=region_of_interest)\n",
+    "usgs_stations = usgs.get_usgs_stations(region=region_of_interest)\n",
+    "usgs_stations"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7c43b139-7bbe-4a2f-9129-6f80b6663356",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "ax2 = us.plot()\n",
+    "usgs_stations.plot(ax=ax2, color='r')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "986e3b94-9d38-44e4-9f51-54d4c86dac63",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "usgs_stations_w_param = usgs_stations[usgs_stations.parm_cd.isin(params_of_interest)]\n",
+    "is_active = (datetime.now() - usgs_stations_w_param.end_date) < timedelta(days=3)\n",
+    "usgs_stations_of_interest = usgs_stations_w_param[is_active]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f8c83ab7-2034-4cd1-9636-c655a1fb2d72",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "ax3 = us.plot()\n",
+    "usgs_stations_of_interest.plot(ax=ax3, color='r')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9c256fc7-cc4f-4b59-b0ee-57abb64d303c",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "usgs_stations_of_interest.columns"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ae5da10b-6718-4ff6-9683-82477f397945",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "data = usgs.get_usgs_data(usgs_metadata=usgs_stations_of_interest)\n",
+    "data = data.assign(vdatum=('site_no', usgs_stations_of_interest.drop_duplicates(subset='site_no').set_index('site_no').alt_datum_cd.loc[data.site_no]))\n",
+    "data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3a30559a-0619-4200-a2b5-b11874a48281",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "def adjust_vdatum(ds):\n",
+    "    # TODO: Adjust \"value\"s based on the \"vdatum\" for the \"site_no\"\n",
+    "    return ds"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c44ea17a-2031-4dcb-89ae-d397086b312a",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "adjusted_data = adjust_vdatum(data)\n",
+    "adjusted_data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e32605ef-28d7-46ba-a49d-3ed98e205fe6",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "plot_data = []\n",
+    "for site_no in adjusted_data.site_no.values:\n",
+    "    for code in adjusted_data.code.values:\n",
+    "        for option in adjusted_data.option.values:\n",
+    "            ds = adjusted_data.sel(site_no=site_no, code=code, option=option).squeeze()\n",
+    "            if ds.value.isnull().all():\n",
+    "                continue\n",
+    "            plot_data.append(ds.to_dataframe())\n",
+    "            \n",
+    "\n",
+    "ncols = 3\n",
+    "fig, axes = pyplot.subplots(ncols=ncols, nrows=len(plot_data) // ncols + 1, figsize=(15, 150))\n",
+    "\n",
+    "for ds, ax in zip(plot_data, axes.ravel()):\n",
+    "    ds.value.dropna().plot(ax=ax, style='x', xlabel='datetime', ylabel=f'{ds.code.iloc[0]} ({ds.unit.iloc[0]})')\n",
+    "fig.tight_layout()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "searvey",
+   "language": "python",
+   "name": "searvey"
+  },
+  "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.10.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}