{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Surface map of METAR data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Table of Contents \n", "\n", "* [Objective](#objective)\n", "* [Strategy](#strategy)\n", "* [Step 0: Import required packages](#step0)\n", "* [Step 1: Browse a THREDDS Data Server (TDS)](#step1)\n", "* [Step 2: Determine the date and hour to gather observations](#step2)\n", "* [Step 3: Surface observations](#step3)\n", " * [Step 3a: Determine catalog location for Surface observations](#step3a)\n", " * [Step 3b: Subset Surface observations](#step3b) \n", " * [Step 3c: Retrieve Surface observations](#step3c)\n", " * [Step 3d: Process Surface observations](#step3d)\n", " * [Step 3e: Visualize Surface observations](#step3e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Objective\n", "> In this notebook, we will make a surface map based on current observations from METAR sites across North America." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 0: Import required packages \n", "[Top](#top)\n", "> But first, we need to import all our required packages. Today we're working with:\n", "> - datetime\n", "> - numpy\n", "> - pandas\n", "> - matplotlib\n", "> - cartopy\n", "> - metpy\n", "> - siphon\n", "\n", "> We will also import a couple of functions that convert weather and cloud cover symbols from METAR files encoded in GEMPAK format." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from datetime import datetime,timedelta\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "\n", "import cartopy.crs as ccrs\n", "import cartopy.feature as cfeature\n", "\n", "from metpy.calc import wind_components, reduce_point_density\n", "from metpy.units import units\n", "from metpy.plots import StationPlot\n", "from metpy.plots.wx_symbols import current_weather, sky_cover, wx_code_map\n", "\n", "from siphon.catalog import TDSCatalog\n", "from siphon.ncss import NCSS" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Load in a collection of functions that process GEMPAK weather conditions and cloud cover data.\n", "%run /kt11/ktyle/python/metargem.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 1: Browse a THREDDS Data Server (TDS) \n", "[Top](#top)\n", "\n", "> A **THREDDS Data Server** provides us with coherent access to a large collection of real-time and archived datasets from a variety of environmental data sources at a number of distributed server sites. \n", "> Various institutions serve data via THREDDS, including our department. You can browse our department's TDS in your web browser using this link: \n", "\n", "http://thredds.atmos.albany.edu:8080/thredds\n", "\n", "\n", "> Let's take a few moments to browse the catalog in a new tab of your web browser." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 2: Determine the date and hour to gather observations " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Use the current time, or set your own for a past time.\n", "# Set current to False if you want to specify a past time.\n", "\n", "nowTime = datetime.utcnow()\n", "current = True\n", "#current = False\n", "if (current):\n", " validTime = datetime.utcnow()\n", " year = validTime.year\n", " month = validTime.month\n", " day = validTime.day\n", " hour = validTime.hour\n", "else:\n", " year = 2024\n", " month = 1\n", " day = 11\n", " hour = 22\n", " \n", "validTime = datetime(year, month, day, hour)\n", "deltaTime = nowTime - validTime\n", "deltaDays = deltaTime.days\n", "\n", "timeStr = validTime.strftime(\"%Y-%m-%d %H UTC\")\n", "timeStr2 = validTime.strftime(\"%Y%m%d%H\")\n", "print(timeStr)\n", "print(validTime)\n", "print(deltaTime)\n", "print(deltaDays)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 3: Surface observations " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 3a: Determine catalog location for Surface observations \n", "[Top](#top)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the date is within the past week, use the current METAR data catalog.\n", "Otherwise, use the archive METAR data catalog." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if (deltaDays <= 7):\n", " # First URL has a subset of stations, just over N. America\n", " NAmCurrent = 'http://thredds.atmos.albany.edu:8080/thredds/catalog/metar/ncdecodedNAm/catalog.xml?dataset=metar/ncdecodedNAm/Metar_NAm_Station_Data_fc.cdmr'\n", " # Second URL has all stations worldwide ... will take longer to load and you definitely need to adjust the point density in step 3d\n", " WorldCurrent = 'http://thredds.atmos.albany.edu:8080/thredds/catalog/metar/ncdecoded/catalog.xml?dataset=metar/ncdecoded/Metar_Station_Data_fc.cdmr'\n", " \n", " metar_cat_url = NAmCurrent\n", "# metar_cat_url = WorldCurrent\n", "else:\n", " metar_cat_url = 'http://thredds.atmos.albany.edu:8080/thredds/catalog/metarArchive/ncdecoded/catalog.xml?dataset=metarArchive/ncdecoded/Archived_Metar_Station_Data_fc.cdmr'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "

Note:

\n", " In general, our METAR archive goes back two years. Requests for older dates may fail. If you need data for an earlier date, email Ross or Kevin.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Create a TDS catalog object from the URL we created above, and examine it.\n", "# The .cdmr extension represents a dataset that works especially well with irregularly-spaced point data sources, such as METAR sites.\n", "\n", "catalog = TDSCatalog(metar_cat_url)\n", "catalog" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# This tells us that there is only one dataset type within this catalog, namely a Feature Collection. We now create an object to store this dataset.\n", "metar_dataset = catalog.datasets['Feature Collection']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### This dataset consists of multiple hours and locations. We are only interested in a subset of this dataset; namely, the most recent hour, and sites within a geographic domain that we will define below. For this, we will exploit a feature of THREDDS that allows for easy subsetting ... the NetCDF Subset Service (NCSS). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ncss_url = metar_dataset.access_urls['NetcdfSubset']\n", "print(ncss_url)\n", "# Import ncss client\n", "ncss = NCSS(ncss_url)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One of the `ncss` object's attributes is `variables`. Let's look at that:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ncss.variables" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are abbreviations that the GEMPAK data analysis/visualization package uses for all the meteorological data potentially encoded into a METAR observation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 3b: Subset Surface Observations \n", "[Top](#top)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Now, let's request all stations within a bounding box for a given time, select certain variables, and create a surface station plot\n", " * Make new NCSS query that specifies a regional domain and a particular time.\n", " * Request data closest to \"now\", or specify a specific YYMMDDHH yourself.\n", " * Specify what variables to retreive\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Set the domain to gather data from and for defining the plot region." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "latN = 55\n", "latS = 20\n", "lonW = -125\n", "lonE = -60" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cLat = (latN + latS)/2\n", "cLon = (lonW + lonE )/2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 3c: Retrieve Surface Observations \n", "[Top](#top)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create an object to build our data query" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "query = ncss.query()\n", "\n", "query.lonlat_box(north=latN-.25, south=latS+0.5, east=lonE-.25, west=lonW+.25)\n", "\n", "# We actually specify a time one minute earlier than the specified hour; this is a quirk of how the data is stored on the THREDDS server and ensures we will get data\n", "# corresponding to the hour we specified\n", "query.time(validTime - timedelta(minutes = 1))\n", "\n", "# Select the variables to query. Note that the variable names depend on the source of the METAR data.\n", "# The 'GEMPAK-like' 4-character names are from the UAlbany THREDDS.\n", "\n", "query.variables('TMPC', 'DWPC', 'PMSL',\n", " 'SKNT', 'DRCT','ALTI','WNUM','VSBY','CHC1', 'CHC2', 'CHC3','CTYH', 'CTYM', 'CTYL' )\n", "query.accept('csv')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pass the query to the THREDDS server's NCSS service. We can then create a Pandas `dataframe` from the returned object." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "data = ncss.get_data(query)\n", "df = pd.DataFrame(data)\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 3d: Process Surface Observations \n", "[Top](#top)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Read in several of the columns; assign / convert units as necessary; convert GEMPAK cloud cover and present wx symbols to MetPy's representation " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lats = data['latitude']\n", "lons = data['longitude']\n", "tair = (data['TMPC'] * units ('degC')).to('degF')\n", "dewp = (data['DWPC'] * units ('degC')).to('degF')\n", "altm = (data['ALTI'] * units('inHg')).to('mbar')\n", "slp = data['PMSL'] * units('hPa')\n", "\n", "# Convert wind to components\n", "u, v = wind_components(data['SKNT'] * units.knots, data['DRCT'] * units.degree)\n", "\n", "# replace missing wx codes or those >= 100 with 0, and convert to MetPy's present weather code\n", "wnum = (np.nan_to_num(data['WNUM'],True).astype(int))\n", "convert_wnum (wnum)\n", "\n", "# Need to handle missing (NaN) and convert to proper code\n", "chc1 = (np.nan_to_num(data['CHC1'],True).astype(int))\n", "chc2 = (np.nan_to_num(data['CHC2'],True).astype(int))\n", "chc3 = (np.nan_to_num(data['CHC3'],True).astype(int))\n", "cloud_cover = calc_clouds(chc1, chc2, chc3)\n", "\n", "# For some reason station id's come back as bytes instead of strings. This line converts the array of station id's into an array of strings.\n", "stid = np.array([s.decode() for s in data['station']])\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next step deals with the removal of overlapping stations, using `reduce_point_density`. This returns a mask we can apply to data to filter the points." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Project points so that we're filtering based on the way the stations are laid out on the map\n", "proj = ccrs.Stereographic(central_longitude=cLon, central_latitude=cLat)\n", "xy = proj.transform_points(ccrs.PlateCarree(), lons, lats)\n", "\n", "# Reduce point density so that there's only one point within a circle whose distance is specified in meters.\n", "# This value will need to change depending on how large of an area you are plotting.\n", "density = 150000\n", "mask = reduce_point_density(xy, density)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Step 3e: Visualize Surface Observations \n", "[Top](#top)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Simple station plotting using plot methods\n", "\n", "One way to create station plots with MetPy is to create an instance of `StationPlot` and call various plot methods, like `plot_parameter`, to plot arrays of data at locations relative to the center point.\n", "\n", "In addition to plotting values, `StationPlot` has support for plotting text strings, symbols, and plotting values using custom formatting.\n", "\n", "Plotting symbols involves mapping integer values to various custom font glyphs in our custom weather symbols font. MetPy provides mappings for converting WMO codes to their appropriate symbol. The `sky_cover` and `current_weather` functions below are two such mappings." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we just plot with `arr[mask]` for every `arr` of data we use in plotting." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Set up a plot with map features\n", "# First set dpi (\"dots per inch\") - higher values will give us a less pixelated final figure.\n", "dpi = 125\n", "\n", "fig = plt.figure(figsize=(15,10), dpi=dpi)\n", "\n", "proj = ccrs.Stereographic(central_longitude=cLon, central_latitude=cLat)\n", "ax = fig.add_subplot(1, 1, 1, projection=proj)\n", "ax.set_facecolor(cfeature.COLORS['water'])\n", "\n", "land_mask = cfeature.NaturalEarthFeature('physical', 'land', '50m',\n", " edgecolor='face',\n", " facecolor=cfeature.COLORS['land'])\n", "lake_mask = cfeature.NaturalEarthFeature('physical', 'lakes', '50m',\n", " edgecolor='face',\n", " facecolor=cfeature.COLORS['water'])\n", "state_borders = cfeature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lakes',\n", " scale='50m', facecolor='none')\n", "\n", "ax.add_feature(land_mask)\n", "ax.add_feature(lake_mask)\n", "ax.add_feature(state_borders, linestyle='solid', edgecolor='black')\n", "\n", "# Slightly reduce the extent of the map as compared to the subsetted data region; this helps eliminate data from being plotted beyond the frame of the map\n", "ax.set_extent ((lonW+2,lonE-4,latS+1,latN-2), crs=ccrs.PlateCarree())\n", "ax.set_extent ((lonW+0.2,lonE-0.2,latS+.1,latN-.1), crs=ccrs.PlateCarree())\n", "\n", "\n", "#If we wanted to add grid lines to our plot:\n", "#ax.gridlines()\n", "\n", "# Create a station plot pointing to an Axes to draw on as well as the location of points\n", "stationplot = StationPlot(ax, lons[mask], lats[mask], transform=ccrs.PlateCarree(),\n", " fontsize=8)\n", "\n", "stationplot.plot_parameter('NW', tair[mask], color='red', fontsize=10)\n", "stationplot.plot_parameter('SW', dewp[mask], color='darkgreen', fontsize=10)\n", "\n", "# Below, we are using a custom formatter to control how the sea-level pressure\n", "# values are plotted. This uses the standard trailing 3-digits of the pressure value\n", "# in tenths of millibars.\n", "stationplot.plot_parameter('NE', slp[mask], color='purple', formatter=lambda v: format(10 * v, '.0f')[-3:])\n", "\n", "\n", "\n", "stationplot.plot_symbol('C', cloud_cover[mask], sky_cover)\n", "stationplot.plot_symbol('W', wnum[mask], current_weather,color='blue',fontsize=12)\n", "stationplot.plot_text((2, 0),stid[mask], color='gray')\n", "#zorder - Higher value zorder will plot the variable on top of lower value zorder. This is necessary for wind barbs to appear. Default is 1.\n", "stationplot.plot_barb(u[mask], v[mask],zorder=2)\n", "\n", "plotTitle = (f\"Sfc Map valid at: {timeStr}\")\n", "ax.set_title (plotTitle);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# In order to see the entire figure, type the name of the figure object below.\n", "fig" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "figName = (f'{timeStr2}_sfmap.png')\n", "fig.savefig(figName)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 Jan. 2024 Environment", "language": "python", "name": "jan24" }, "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.7" } }, "nbformat": 4, "nbformat_minor": 4 }