{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "
\n", "\"Unidata\n", "
\n", "\n", "# Precipitation Accumulation Region of Interest\n", "**Python-AWIPS Tutorial Notebook**\n", "\n", "
\n", "
\n", "\n", "---\n", "\n", "
\"Colorized
\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": [ "

Table of Contents

\n", "
" ] }, { "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": [ "Top\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": [ "Top\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": [ "Top\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": [ "Top\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": [ "Top\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": [ "Top\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": [ "Top\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": [ "Top\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 }