python-awips/examples/notebooks/Precip_Accumulation_Region_of_Interest.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a name=\"top\"></a>\n",
    "<div style=\"width:1000 px\">\n",
    "\n",
    "<div style=\"float:right; width:98 px; height:98px;\">\n",
    "<img src=\"https://docs.unidata.ucar.edu/images/logos/unidata_logo_vertical_150x150.png\" alt=\"Unidata Logo\" style=\"height: 98px;\">\n",
    "</div>\n",
    "\n",
    "# Precipitation Accumulation Region of Interest\n",
    "**Python-AWIPS Tutorial Notebook**\n",
    "\n",
    "<div style=\"clear:both\"></div>\n",
    "</div>\n",
    "\n",
    "---\n",
    "\n",
    "<div style=\"float:right; width:250 px\"><img src=\"../images/precip_roi_preview.png\" alt=\"Colorized Precipation Accumulation with Identified Region of Interest\" style=\"height: 300px;\"></div>\n",
    "\n",
    "\n",
    "# Objectives\n",
    "\n",
    "* Access the model data from an EDEX server and limit the data returned by using model specific parameters\n",
    "* Calculate the total precipitation over several model runs\n",
    "* Create a colorized plot for the continental US of the accumulated precipitation data\n",
    "* Calculate and identify area of highest of precipitation\n",
    "* Use higher resolution data to draw region of interest\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Imports\" data-toc-modified-id=\"Imports-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Imports</a></span></li><li><span><a href=\"#Initial-Setup\" data-toc-modified-id=\"Initial-Setup-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Initial Setup</a></span><ul class=\"toc-item\"><li><span><a href=\"#Geographic-Filter\" data-toc-modified-id=\"Geographic-Filter-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Geographic Filter</a></span></li><li><span><a href=\"#EDEX-Connection\" data-toc-modified-id=\"EDEX-Connection-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>EDEX Connection</a></span></li><li><span><a href=\"#Refine-the-Request\" data-toc-modified-id=\"Refine-the-Request-2.3\"><span class=\"toc-item-num\">2.3&nbsp;&nbsp;</span>Refine the Request</a></span></li><li><span><a href=\"#Get-Times\" data-toc-modified-id=\"Get-Times-2.4\"><span class=\"toc-item-num\">2.4&nbsp;&nbsp;</span>Get Times</a></span></li></ul></li><li><span><a href=\"#Function:-calculate_accumulated_precip()\" data-toc-modified-id=\"Function:-calculate_accumulated_precip()-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Function: calculate_accumulated_precip()</a></span></li><li><span><a href=\"#Fuction:-make_map()\" data-toc-modified-id=\"Fuction:-make_map()-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Fuction: make_map()</a></span></li><li><span><a href=\"#Get-the-Data!\" data-toc-modified-id=\"Get-the-Data!-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Get the Data!</a></span></li><li><span><a href=\"#Plot-the-Data!\" data-toc-modified-id=\"Plot-the-Data!-6\"><span class=\"toc-item-num\">6&nbsp;&nbsp;</span>Plot the Data!</a></span><ul class=\"toc-item\"><li><span><a href=\"#Create-CONUS-Image\" data-toc-modified-id=\"Create-CONUS-Image-6.1\"><span class=\"toc-item-num\">6.1&nbsp;&nbsp;</span>Create CONUS Image</a></span></li><li><span><a href=\"#Create-Region-of-Interest-Image\" data-toc-modified-id=\"Create-Region-of-Interest-Image-6.2\"><span class=\"toc-item-num\">6.2&nbsp;&nbsp;</span>Create Region of Interest Image</a></span></li></ul></li><li><span><a href=\"#High-Resolution-ROI\" data-toc-modified-id=\"High-Resolution-ROI-7\"><span class=\"toc-item-num\">7&nbsp;&nbsp;</span>High Resolution ROI</a></span><ul class=\"toc-item\"><li><span><a href=\"#New-Data-Request\" data-toc-modified-id=\"New-Data-Request-7.1\"><span class=\"toc-item-num\">7.1&nbsp;&nbsp;</span>New Data Request</a></span></li><li><span><a href=\"#Calculate-Data\" data-toc-modified-id=\"Calculate-Data-7.2\"><span class=\"toc-item-num\">7.2&nbsp;&nbsp;</span>Calculate Data</a></span></li><li><span><a href=\"#Plot-ROI\" data-toc-modified-id=\"Plot-ROI-7.3\"><span class=\"toc-item-num\">7.3&nbsp;&nbsp;</span>Plot ROI</a></span></li></ul></li><li><span><a href=\"#See-Also\" data-toc-modified-id=\"See-Also-8\"><span class=\"toc-item-num\">8&nbsp;&nbsp;</span>See Also</a></span><ul class=\"toc-item\"><li><span><a href=\"#Related-Notebooks\" data-toc-modified-id=\"Related-Notebooks-8.1\"><span class=\"toc-item-num\">8.1&nbsp;&nbsp;</span>Related Notebooks</a></span></li><li><span><a href=\"#Additional-Documentation\" data-toc-modified-id=\"Additional-Documentation-8.2\"><span class=\"toc-item-num\">8.2&nbsp;&nbsp;</span>Additional Documentation</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports\n",
    "\n",
    "The imports below are used throughout the notebook. Note the first import is coming directly from python-awips and allows us to connect to an EDEX server. The subsequent imports are for data manipulation and visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from awips.dataaccess import DataAccessLayer\n",
    "import cartopy.crs as ccrs\n",
    "import matplotlib.pyplot as plt\n",
    "from metpy.units import units\n",
    "import numpy as np\n",
    "from shapely.geometry import Point, Polygon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"#top\">Top</a>\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geographic Filter\n",
    "\n",
    "By defining a bounding box for the Continental US (CONUS), we’re able to optimize the data request sent to the EDEX server."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "conus=[-125, -65, 25, 55]\n",
    "conus_envelope = Polygon([(conus[0],conus[2]),(conus[0],conus[3]),\n",
    "                          (conus[1],conus[3]),(conus[1],conus[2]),\n",
    "                          (conus[0],conus[2])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EDEX Connection\n",
    "\n",
    "First we establish a connection to Unidata’s public EDEX server. With that connection made, we can create a [new data request object](http://unidata.github.io/python-awips/api/IDataRequest.html) and set the data type to **grid**, and use the geographic envelope we just created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "DataAccessLayer.changeEDEXHost(\"edex-beta.unidata.ucar.edu\")\n",
    "request = DataAccessLayer.newDataRequest(\"grid\", envelope=conus_envelope)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Refine the Request\n",
    "\n",
    "Here we specify which model we're interested in by setting the *LocationNames*, and the specific data we're interested in by setting the *Levels* and *Parameters*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "request.setLocationNames(\"GFS1p0\")\n",
    "request.setLevels(\"0.0SFC\")\n",
    "request.setParameters(\"TP\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get Times\n",
    "\n",
    "We need to get the available times and cycles for our model data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[5], line 3\u001b[0m\n\u001b[1;32m      1\u001b[0m cycles \u001b[38;5;241m=\u001b[39m DataAccessLayer\u001b[38;5;241m.\u001b[39mgetAvailableTimes(request, \u001b[38;5;28;01mTrue\u001b[39;00m)\n\u001b[1;32m      2\u001b[0m times \u001b[38;5;241m=\u001b[39m DataAccessLayer\u001b[38;5;241m.\u001b[39mgetAvailableTimes(request)\n\u001b[0;32m----> 3\u001b[0m fcstRun \u001b[38;5;241m=\u001b[39m DataAccessLayer\u001b[38;5;241m.\u001b[39mgetForecastRun(\u001b[43mcycles\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m]\u001b[49m, times)\n",
      "\u001b[0;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "cycles = DataAccessLayer.getAvailableTimes(request, True)\n",
    "times = DataAccessLayer.getAvailableTimes(request)\n",
    "fcstRun = DataAccessLayer.getForecastRun(cycles[-1], times)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"#top\">Top</a>\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Function: calculate_accumulated_precip()\n",
    "\n",
    "Since we'll want to calculate the accumulated precipitation of our data more than once, it makes sense to create a function that we can call instead of duplicating the logic.\n",
    "\n",
    "This function cycles through all the grid data responses and adds up all of the rainfall to produce a numpy array with the total ammount of rainfall for the given data request.  It also finds the maximum rainfall point in x and y coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_accumulated_precip(dataRequest, forecastRun):\n",
    "\n",
    "    for i, tt in enumerate(forecastRun):\n",
    "        response = DataAccessLayer.getGridData(dataRequest, [tt])\n",
    "        grid = response[0]\n",
    "        if i>0:\n",
    "            data += grid.getRawData()\n",
    "        else:\n",
    "            data = grid.getRawData()\n",
    "        data[data <= -9999] = 0\n",
    "        print(data.min(), data.max(), grid.getDataTime().getFcstTime()/3600)\n",
    "\n",
    "    # Convert from mm to inches\n",
    "    result = (data * units.mm).to(units.inch)\n",
    "    \n",
    "    ii,jj = np.where(result==result.max())\n",
    "    i=ii[0]\n",
    "    j=jj[0]\n",
    "\n",
    "    return result, i, j"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"#top\">Top</a>\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fuction: make_map()\n",
    "\n",
    "This function creates the basics of the map we're going to plot our data on. It takes in a bounding box to determine the extent and then adds coastlines for easy frame of reference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_map(bbox, projection=ccrs.PlateCarree()):\n",
    "    fig, ax = plt.subplots(figsize=(20, 14),\n",
    "            subplot_kw=dict(projection=projection))\n",
    "    ax.set_extent(bbox)\n",
    "    ax.coastlines(resolution='50m')\n",
    "    return fig, ax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"#top\">Top</a>\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the Data!\n",
    "\n",
    "Access the data from the DataAccessLayer interface using the *getGridData* function.  Use that data to calculate the accumulated rainfall, the maximum rainfall point, and the region of interest bounding box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## get the grid response from edex\n",
    "response = DataAccessLayer.getGridData(request, [fcstRun[-1]])    \n",
    "## take the first result to get the location information from\n",
    "grid = response[0]\n",
    "\n",
    "## get the location coordinates and create a bounding box for our map\n",
    "lons, lats = grid.getLatLonCoords()\n",
    "bbox = [lons.min(), lons.max(), lats.min(), lats.max()]\n",
    "fcstHr = int(grid.getDataTime().getFcstTime()/3600)\n",
    "\n",
    "## calculate the total precipitation\n",
    "tp_inch, i, j = calculate_accumulated_precip(request, fcstRun)\n",
    "print(tp_inch.min(), tp_inch.max())\n",
    "\n",
    "## use the max points coordinates to get the max point in lat/lon coords\n",
    "maxPoint = Point(lons[i][j], lats[i][j])\n",
    "inc = 3.5\n",
    "## create a region of interest bounding box\n",
    "roi_box=[maxPoint.x-inc, maxPoint.x+inc, maxPoint.y-inc, maxPoint.y+inc]\n",
    "roi_polygon = Polygon([(roi_box[0],roi_box[2]),(roi_box[0],roi_box[3]), \n",
    "                (roi_box[1],roi_box[3]),(roi_box[1],roi_box[2]),(roi_box[0],roi_box[2])])\n",
    "\n",
    "print(maxPoint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"#top\">Top</a>\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the Data!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create CONUS Image\n",
    "\n",
    "Plot our data on our CONUS map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = plt.get_cmap('rainbow')\n",
    "fig, ax = make_map(bbox=bbox)\n",
    "cs = ax.pcolormesh(lons, lats, tp_inch, cmap=cmap)\n",
    "cbar = fig.colorbar(cs, shrink=0.7, orientation='horizontal')\n",
    "cbar.set_label(grid.getLocationName() + \" Total accumulated precipitation in inches, \" \\\n",
    "               + str(fcstHr) + \"-hr fcst valid \" + str(grid.getDataTime().getRefTime()))\n",
    "\n",
    "ax.scatter(maxPoint.x, maxPoint.y, s=300,\n",
    "           transform=ccrs.PlateCarree(),marker=\"+\",facecolor='black')\n",
    "\n",
    "ax.add_geometries([roi_polygon], ccrs.PlateCarree(), facecolor='none', edgecolor='white', linewidth=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create Region of Interest Image\n",
    "\n",
    "Now crop the data and zoom in on the region of interest (ROI) to create a new plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cmap = plt.get_cmap('rainbow')\n",
    "fig, ax = make_map(bbox=roi_box)\n",
    "\n",
    "cs = ax.pcolormesh(lons, lats, tp_inch, cmap=cmap)\n",
    "cbar = fig.colorbar(cs, shrink=0.7, orientation='horizontal')\n",
    "cbar.set_label(grid.getLocationName() + \" Total accumulated precipitation in inches, \" \\\n",
    "               + str(fcstHr) + \"-hr fcst valid \" + str(grid.getDataTime().getRefTime()))\n",
    "\n",
    "ax.scatter(maxPoint.x, maxPoint.y, s=300,\n",
    "           transform=ccrs.PlateCarree(),marker=\"+\",facecolor='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"#top\">Top</a>\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## High Resolution ROI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### New Data Request\n",
    "\n",
    "To see the region of interest more clearly, we can redo the process with a higher resolution model (GFS20 vs. GFS1.0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roiRequest = DataAccessLayer.newDataRequest(\"grid\", envelope=conus_envelope)\n",
    "roiRequest.setLocationNames(\"GFS20\")\n",
    "roiRequest.setLevels(\"0.0SFC\")\n",
    "roiRequest.setParameters(\"TP\")\n",
    "\n",
    "roiCycles = DataAccessLayer.getAvailableTimes(roiRequest, True)\n",
    "roiTimes = DataAccessLayer.getAvailableTimes(roiRequest)\n",
    "roiFcstRun = DataAccessLayer.getForecastRun(roiCycles[-1], roiTimes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roiResponse = DataAccessLayer.getGridData(roiRequest, [roiFcstRun[-1]])  \n",
    "print(roiResponse)\n",
    "roiGrid = roiResponse[0]\n",
    "\n",
    "roiLons, roiLats = roiGrid.getLatLonCoords()\n",
    "\n",
    "roi_data, i, j = calculate_accumulated_precip(roiRequest, roiFcstRun)\n",
    "\n",
    "roiFcstHr = int(roiGrid.getDataTime().getFcstTime()/3600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot ROI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cmap = plt.get_cmap('rainbow')\n",
    "fig, ax = make_map(bbox=roi_box)\n",
    "\n",
    "cs = ax.pcolormesh(roiLons, roiLats, roi_data, cmap=cmap)\n",
    "cbar = fig.colorbar(cs, shrink=0.7, orientation='horizontal')\n",
    "cbar.set_label(roiGrid.getLocationName() + \" Total accumulated precipitation in inches, \" \\\n",
    "               + str(roiFcstHr) + \"-hr fcst valid \" + str(roiGrid.getDataTime().getRefTime()))\n",
    "\n",
    "ax.scatter(maxPoint.x, maxPoint.y, s=300,\n",
    "           transform=ccrs.PlateCarree(),marker=\"+\",facecolor='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"#top\">Top</a>\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## See Also\n",
    "\n",
    "### Related Notebooks\n",
    "\n",
    "* [Colorized Grid Data](https://unidata.github.io/python-awips/examples/generated/Colorized_Grid_Data.html)\n",
    "* [Grid Levels and Parameters](https://unidata.github.io/python-awips/examples/generated/Grid_Levels_and_Parameters.html)\n",
    "\n",
    "### Additional Documentation\n",
    "\n",
    "**python-awips:**\n",
    "* [awips.DataAccessLayer](http://unidata.github.io/python-awips/api/DataAccessLayer.html)\n",
    "* [awips.PyGridData](http://unidata.github.io/python-awips/api/PyGridData.html)\n",
    "\n",
    "**matplotlib:**\n",
    "* [matplotlib.pyplot](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.html)\n",
    "* [matplotlib.pyplot.subplot](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.subplot.html)\n",
    "* [matplotlib.pyplot.pcolormesh](https://matplotlib.org/3.3.3/api/_as_gen/matplotlib.pyplot.pcolormesh.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"#top\">Top</a>\n",
    "\n",
    "---"
   ]
  }
 ],
 "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.11.11"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}