From 35be20484036432f20c902381e7935c49db9c9bb Mon Sep 17 00:00:00 2001 From: "J.R. Leeman" Date: Fri, 16 Nov 2018 14:28:29 -0700 Subject: [PATCH 1/3] Refactor satellite notebook. --- ...Data.ipynb => PlottingSatelliteData.ipynb} | 442 +++++++++--------- .../Satellite_Data/solutions/data_url.py | 19 + notebooks/Satellite_Data/solutions/sat_map.py | 15 + 3 files changed, 251 insertions(+), 225 deletions(-) rename notebooks/Satellite_Data/{Working with Satellite Data.ipynb => PlottingSatelliteData.ipynb} (53%) create mode 100644 notebooks/Satellite_Data/solutions/data_url.py create mode 100644 notebooks/Satellite_Data/solutions/sat_map.py diff --git a/notebooks/Satellite_Data/Working with Satellite Data.ipynb b/notebooks/Satellite_Data/PlottingSatelliteData.ipynb similarity index 53% rename from notebooks/Satellite_Data/Working with Satellite Data.ipynb rename to notebooks/Satellite_Data/PlottingSatelliteData.ipynb index d83d2d7e..79dc8110 100644 --- a/notebooks/Satellite_Data/Working with Satellite Data.ipynb +++ b/notebooks/Satellite_Data/PlottingSatelliteData.ipynb @@ -4,12 +4,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "\n", + "\n", "
\n", "\n", "
\"Unidata
\n", "\n", - "

Working with Satellite Data

\n", + "

Plotting Satellite Data

\n", "

Unidata Python Workshop

\n", "\n", "
\n", @@ -26,17 +26,16 @@ "* **Exercises:** 30 minutes\n", "\n", "### Questions\n", + "1. Where are current GOES data available?\n", "1. How can satellite data be obtained with Siphon?\n", - "1. What format is satellite data distributed in?\n", - "1. How can MetPy be used to open and explore satellite data?\n", + "1. How can MetPy simplify metadata parsing?\n", "1. How can maps of satellite data be made?\n", - "1. How can other datasets be integrated with satellite data?\n", "\n", - "### Objectives\n", - "1. Download satellite data with Siphon\n", - "1. Parse out the file\n", - "1. Make a plot of the data\n", - "1. Use colortables and annotations\n", + "### Table of Contents\n", + "1. Finding GOES data\n", + "1. Accessing data with Siphon\n", + "1. Digging into and parsing the data\n", + "1. Plotting the data\n", "1. Bonus: Animations" ] }, @@ -44,43 +43,33 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "\n", - "## 1. Download satellite data with Siphon\n", - "\n", - "The first step is to find the satellite data. Normally, we would browse over to http://thredds.ucar.edu/thredds/ and find the top-level [THREDDS Data Server (TDS)](https://www.unidata.ucar.edu/software/thredds/current/tds/TDS.html) catalog. From there we can drill down to find satellite data products.\n", - "\n", - "For this tutorial, we want to utilize the new GOES-16 imagery. This is found at http://thredds-test.unidata.ucar.edu/. For current data, you could navigate to the `Test Datasets` directory, then `GOES-16 Products` and `GOES-16`. There are subfolders for the CONUS, full disk, mesoscale sector images, and other products. In each of these is a folder for each [channel of the ABI](http://www.goes-r.gov/education/ABI-bands-quick-info.html). In each channel there is a folder for every day in the month-long rolling archive. As you can see, there is a massive amount of data coming down from GOES-16! For the workshop we will be using archived case study data from Irma located at [http://thredds-test.unidata.ucar.edu/thredds/catalog/casestudies/irma/goes16/catalog.html](http://thredds-test.unidata.ucar.edu/thredds/catalog/casestudies/irma/goes16/catalog.html).\n", - "\n", - "We could download the files to our computers from here, but that can become tedious for downloading many files, requires us to store them on our computer, and does not lend itself to automation.\n", - "\n", - "We can use Unidata's [Siphon](https://github.com/Unidata/siphon) package to parse the catalog from the TDS. This provides us a nice programmatic way of accessing the data. We start by importing the `TDSCatalog` class from siphon and giving it the URL to the catalog we just surfed to manually. **Note:** Instead of giving it the link to the HTML catalog, we change the extension to XML, which asks the TDS for the XML version of the catalog. This is much better to work with in code. If you forget, the extension will be changed for you with a warning being issued from siphon.\n", + "
\n", "\n", - "Next we'll build a URL to get to the current data for a specified channel and region:" + "\n", + "## Finding GOES Data" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "from datetime import datetime\n", - "import metpy\n", - "from siphon.catalog import TDSCatalog\n", + "The first step is to find the satellite data. Normally, we would browse over to http://thredds.ucar.edu/thredds/ and find the top-level [THREDDS Data Server (TDS)](https://www.unidata.ucar.edu/software/thredds/current/tds/TDS.html) catalog. From there we can drill down to find satellite data products. For some technical reasons, the new GOES series data is currently found on the server thredds-test at http://thredds-test.unidata.ucar.edu/ - it will be moved to the main thredds server in the near future. \n", "\n", - "date = datetime.utcnow()\n", - "channel = 8\n", - "region = 'Mesoscale-2'\n", + "For current data, you could navigate to the `Test Datasets` directory, then `GOES Products` and `GOES-16`. There are subfolders for the CONUS, full disk, mesoscale sector images, and other products. In each of these is a folder for each [channel of the ABI](http://www.goes-r.gov/education/ABI-bands-quick-info.html). In each channel there is a folder for every day in the approximately month-long rolling archive. As you can see, there are a massive amount of data coming down from GOES-16!\n", "\n", - "cat = TDSCatalog('http://thredds-test.unidata.ucar.edu/thredds/catalog/satellite/goes16/GOES16/'\n", - " '{}/Channel{:02d}/{:%Y%m%d}/catalog.xml'.format(region, channel, date))" + "In the next section we'll be downloading the data in a pythonic way, so our first task is to build a URL that matches the URL we manually navigated to in the web browser. To make it as flexible as possible, we'll want to use variables for the sector name (CONUS, full-disk, mesoscale-1, etc.), the date, and the ABI channel number." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We now have a `TDSCatalog` object called `cat` that we can examine and use to get handles to work with the data. To find the latest file, we can look at the `cat.datasets` attribute. Let’s look at the last five datasets:" + "### Exercise\n", + "\n", + "* Create variables named `image_date`, `region`, and `channel`. Assign them to today's date, the mesoscale-1 region, and ABI channel 8. \n", + "* Construct a string `data_url` from these variables and the URL we navigated to above.\n", + "* Verify that following your link will take you where you think it should.\n", + "* Change the extension from `catalog.html` to `catalog.xml`. What do you see?" ] }, { @@ -89,14 +78,23 @@ "metadata": {}, "outputs": [], "source": [ - "cat.datasets[-5:]" + "from datetime import datetime\n", + "\n", + "# Create variables for URL generation\n", + "# YOUR CODE GOES HERE\n", + "\n", + "# Construct the data_url string\n", + "# YOUR CODE GOES HERE\n", + "\n", + "# Print out your URL and verify it works!\n", + "# YOUR CODE GOES HERE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We'll get the next to most recent dataset (sometimes the most recent will not have received all tiles yet) using the OPENDAP protocol." + "#### Solution" ] }, { @@ -105,69 +103,59 @@ "metadata": {}, "outputs": [], "source": [ - "ds = cat.datasets[-2]\n", - "ds = ds.remote_access(use_xarray=True)" + "# %load solutions/data_url.py\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now we'll use xarray to open the dataset as a netCDF data store." + "Top\n", + "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "
\n", - " EXERCISE:\n", - "
    \n", - "
  • Using the link to the Hurricane Irma archive above, get data from September 10, 2017 at 15Z for the mid-level water vapor channel (9) in the Mesoscale-1 sector. You could search through the datasets list manually, but instead try the `cat.datasets.filter_time_nearest` functionality! Go through opening the xarray dataset.
    • \n", - "
\n", - "
" + "\n", + "## Accessing data with Siphon" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "# Your code goes here\n" + "We could download the files to our computers from the THREDDS web interface, but that can become tedious for downloading many files, requires us to store them on our computer, and does not lend itself to automation.\n", + "\n", + "We can use [Siphon](https://github.com/Unidata/siphon) parse the catalog from the TDS. This provides us a nice programmatic way of accessing the data. We start by importing the `TDSCatalog` class from siphon and giving it the URL to the catalog we just surfed to manually. **Note:** Instead of giving it the link to the HTML catalog, we change the extension to XML, which asks the TDS for the XML version of the catalog. This is much better to work with in code. If you forget, the extension will be changed for you with a warning being issued from siphon.\n", + "\n", + "We want to create a `TDSCatalog` object called `cat` that we can examine and use to get handles to work with the data." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "\n", - "
\n", - "
\n",
-    "cat = TDSCatalog('http://thredds-test.unidata.ucar.edu/thredds/catalog/casestudies/irma/goes16/Mesoscale-1/Channel09/20170910/catalog.xml')\n",
-    "ds = cat.datasets.filter_time_nearest(datetime(2017, 9, 10, 15))\n",
-    "print(ds.name)\n",
-    "ds = ds.remote_access(use_xarray=True)\n",
-    "
\n", - "
" + "from siphon.catalog import TDSCatalog" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "Top\n", - "
" + "cat = TDSCatalog(data_url)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "\n", - "## 2. Parse out the file\n", - "\n", - "We have an xarray dataset now that contains all of the data for a given satellite image. We need to explore the variables available in the file and pull out the useful parts that we need to make a map." + "To find the latest file, we can look at the `cat.datasets` attribute. Let’s look at the last five datasets:" ] }, { @@ -176,17 +164,14 @@ "metadata": {}, "outputs": [], "source": [ - "list(ds)" + "cat.datasets[-5:]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Our goal is to plot the imagery, so we'll be using the `Sectorized_CMI` variable from the `.variables` dictionary.\n", - "Rather than just giving back the raw array of data, this gives back a `DataArray` object; from here not only\n", - "can we get the raw data values, but there are useful metadata as well. We can see just what additional information\n", - "is present by printing out the `DataArray` object:" + "We'll get the next to most recent dataset (sometimes the most recent will not have received all tiles yet) and store it as variable `dataset`. Note that we haven't actually downloaded or transferred any real data yet, just bits of metadata have been received from THREDDS and parsed by siphon." ] }, { @@ -195,58 +180,55 @@ "metadata": {}, "outputs": [], "source": [ - "print(ds['Sectorized_CMI'])" + "dataset = cat.datasets[-2]" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "This reveals several useful pieces of information, but we're going to focus on the `grid_mapping` data. This attribute is defined by the [NetCDF Climate and Forecast (CF) Metadata Conventions](http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html) and specifies a variable that contains information about the grid's projection." + "print(dataset)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The projection is generally [Lambert conformal conic](https://en.wikipedia.org/wiki/Lambert_conformal_conic_projection) or [Mercator](https://en.wikipedia.org/wiki/Mercator_projection) depending on the sector. The meta-data includes a few parameters (such as the latitude and longitude of the origin) needed to properly set up the projection to match what was used to create the image. There is also has information about the assumed shape of the Earth. We can let MetPy parse the information and build a projection for us!" + "We're finally ready to get the actual data. We could download the file, then open that, but there is no need! We can use siphon to help us only get what we need and hold it in memory. Notice that we're using the XArray accessor which will make life much nicer that dealing with the raw netCDF (like we used to back in the days of early 2018)." ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ - "dat = ds.metpy.parse_cf('Sectorized_CMI')\n", - "proj = dat.metpy.cartopy_crs" + "ds = dataset.remote_access(use_xarray=True)" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "print(dat)" + "Top\n", + "
" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "print(proj)" + "\n", + "## Digging into and parsing the data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We can then pull out the x and y coordinate information easily:" + "Now that we've got some data - let's see what we actually got our hands on." ] }, { @@ -255,26 +237,26 @@ "metadata": {}, "outputs": [], "source": [ - "x = dat['x']\n", - "y = dat['y']" + "ds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Top\n", - "
" + "Great, so we have an XArray Dataset object, something we've dealt with before! We also see that we have the coordinates `time`, `y`, and `x` as well as the data variables of `Sectorized_CMI` and the projection information.\n", + "\n", + "The first thing we can do is get an XArray DataArray with the CF compliant metadata parsed out to do things like automatically generate our cartopy crs for us!" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "\n", - "## 3. Make a plot of the data\n", - "\n", - "We are finally ready to plot the satellite data! Using `imshow()` we can get an image, but there is a lot of work to do here. The data need to be projected to a meaningful representation, map outlines added, and annotations added." + "import metpy # This gives us the metpy data accessor\n", + "dat = ds.metpy.parse_cf('Sectorized_CMI')" ] }, { @@ -283,16 +265,14 @@ "metadata": {}, "outputs": [], "source": [ - "# Make sure the notebook puts figures inline\n", - "%matplotlib inline\n", - "import matplotlib.pyplot as plt" + "dat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "CartoPy's projections are designed to interface with matplotlib, so they can just be passed as the `projection` keyword argument when creating an `Axes` using the `add_subplot` method. Since the x and y coordinates, as well as the image data, are referenced in the projection we are making the map in, we can pass all of them directly to plotting methods (such as `imshow`) with no additional information. The `extent` keyword argument to `imshow` is used to specify the bounds of the image data being plotted. It is **especially important** to specify that the `origin` is at the upper left of the image (standard practice in imagery). If your forget this, your image will be flipped. Try it!" + "We can see some useful metadata are attached to this DataArray - including a coordinate refernce system. Let's look at what it is:" ] }, { @@ -301,26 +281,14 @@ "metadata": {}, "outputs": [], "source": [ - "# Create a new figure with size 10\" by 10\"\n", - "fig = plt.figure(figsize=(10, 10))\n", - "\n", - "# Put a single axes on this figure; set the projection for the axes to be our\n", - "# Lambert conformal projection\n", - "ax = fig.add_subplot(1, 1, 1, projection=proj)\n", - "\n", - "# Plot the data with mostly the defaults\n", - "# Note, we save the image returned by imshow for later...\n", - "im = ax.imshow(dat, extent=(x.min(), x.max(), y.min(), y.max()), origin='upper')\n", - "\n", - "# Add high-resolution coastlines to the plot\n", - "ax.coastlines(resolution='50m', color='black')" + "dat.metpy.crs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "This is a nice start, but it would be nice to have better geographic references for the image. For example, where are the states? Cartopy's `feature` module has support for adding geographic features to plots, with many features are built in. The `BORDERS` built-in feature contains country borders and `STATES` contains US state boundaries. There is also support for creating \"custom\" features from the [Natural Earth](http://www.naturalearthdata.com/) set of free vector and raster map data (CartoPy will automatically download the necessary data and cache it locally)." + "Ok - so we have as `CFProjection` object - but we really want a cartopy compatible crs. Metpy provides a cartopy_crs accessor to do the translation for you:" ] }, { @@ -329,43 +297,23 @@ "metadata": {}, "outputs": [], "source": [ - "import cartopy.feature as cfeature\n", - "\n", - "# Add country borders with a thick line.\n", - "ax.add_feature(cfeature.BORDERS, linewidth=2, edgecolor='black')\n", - "\n", - "# Add state borders with dotted lines, denoted by ':'\n", - "ax.add_feature(cfeature.STATES, linestyle=':', edgecolor='black')\n", - "\n", - "# Redisplay modified figure\n", - "fig" + "proj = dat.metpy.cartopy_crs" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "Top\n", - "
" + "proj" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "\n", - "## 4. Use colortables and annotations\n", - "\n", - "The map is much improved now, but it would look much better with a different color scheme. \n", - "\n", - "Colormapping in matplotlib (which backs CartoPy) is handled through two pieces:\n", - "\n", - "- The colormap controls how values are converted from floating point values in the range [0, 1] to colors (think colortable)\n", - "- The norm (normalization) controls how data values are converted to floating point values in the range [0, 1]\n", - "\n", - "Let's start by setting the colormap to be black and white and normalizing the data to get the best contrast. We'll make a histogram to see the distribution of values in the data, then clip that range down to enhance contrast in the data visualization. **Note:** `cmap` and `norm` can also be set during the `imshow` call as keyword arguments.\n", - "\n", - "We `flatten` the 2D array down to 1D for the histogram and remove any masked elements with `compressed`." + "Before we get to plotting, let's go ahead and get a handle to the x and y coordinates that we'll need. This is easy to do with the accessor on our data variable DataArray." ] }, { @@ -374,7 +322,8 @@ "metadata": {}, "outputs": [], "source": [ - "plt.hist(dat.to_masked_array().compressed().flatten(), bins=255);" + "x = dat.metpy.x\n", + "y = dat.metpy.y" ] }, { @@ -383,39 +332,37 @@ "metadata": {}, "outputs": [], "source": [ - "# Set colormap\n", - "im.set_cmap('Greys')\n", - "\n", - "# Set norm\n", - "im.set_norm(plt.Normalize(200,255))\n", - "\n", - "# Show figure again\n", - "fig" + "x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "In meteorology, we have many ‘standard’ colortables that have been used for certain types of data. We have included these in Metpy in the `metpy.plots.ctables` module. By importing the `ColortableRegistry` we gain access to the colortables, as well as handy normalizations to go with them. We can see the colortables available by looking at the dictionary keys. The colortables ending in `_r` are the reversed versions of the named colortable." + "Metadata travels with these coordinates!" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "from metpy.plots import colortables\n", - "\n", - "[table for table in colortables.keys()]" + "Top\n", + "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Let’s use the `WVCIMSS` colormap, a direct conversion of the GEMPAK colormap. The code below asks for the colormap, as well as a normalization that covers the range 195 to 265. This was empirically determined to closely match other displays of water vapor data. We then apply it to the existing image we have been working with." + "\n", + "## Plotting the Data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To plot our data we'll use the matplotlib method `imshow`. This function maps a 2D array of values to a 2D visual representation of those values, commonly known as a picture. It can even handle arrays that are MxNx3 or MxNx4 for RGB pictures (think digital pictures) or RGBA data. Let's look at a simple example:" ] }, { @@ -424,17 +371,23 @@ "metadata": {}, "outputs": [], "source": [ - "wv_norm, wv_cmap = colortables.get_with_range('WVCIMSS_r', 195, 265)\n", - "im.set_cmap(wv_cmap)\n", - "im.set_norm(wv_norm)\n", - "fig" + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "%matplotlib inline" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "One more thing that would be nice is putting the date and time on the image, so let's do that. First grab the `start_date_time` attribute from the dataset object:" + "a = np.array([[1, 0.5, 0.25, 0, 0, 0],\n", + " [0.5, 1, 0.5, 0.25, 0, 0],\n", + " [0.25, 0.5, 1, 0.5, 0.25, 0],\n", + " [0, 0.25, 0.5, 1, 0.5, 0.25],\n", + " [0, 0, 0.25, 0.5, 1, 0.5],\n", + " [0, 0, 0, 0.25, 0.5, 1]])" ] }, { @@ -443,14 +396,18 @@ "metadata": {}, "outputs": [], "source": [ - "ds.start_date_time" + "fig, ax = plt.subplots(1, 1, figsize=(10, 8))\n", + "ax.imshow(a)" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "This timestamp contains the year, Julian day, hour, minute, and second that the image was collected. We will use `datetime.strptime` to parse this out into a Python `datetime` object." + "lats = np.array([50, 48, 46, 44, 42, 40])\n", + "lons = np.array([-136, -134, -132, -130, -128, -126])" ] }, { @@ -459,17 +416,42 @@ "metadata": {}, "outputs": [], "source": [ - "timestamp = datetime.strptime(ds.start_date_time, '%Y%j%H%M%S')\n", - "print(timestamp)" + "import cartopy.crs as ccrs\n", + "fig = plt.figure(figsize=(10, 8))\n", + "ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())\n", + "ax.imshow(a)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import cartopy.crs as ccrs\n", + "fig = plt.figure(figsize=(10, 8))\n", + "ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())\n", + "ax.imshow(a, extent=(lons[0], lons[-1], lats[0], lats[-1]),\n", + " origin='upper')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can see this is really a visual representation of the array data - it's also a handy way to quickly look at data, especially things like resolution matricies, etc." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "A Python `datetime` object is easy to work with. Let's add it to our plot.\n", + "### Exercise\n", "\n", - "MetPy has a built in `add_timestamp` method that makes this easy, as well as adding some custom text. We'll be changing this next, so we grab a handle to the created text so it can be removed later." + "* Using what you've learned about `imshow` plot the data array dat with it.\n", + "* Add the cartopy states and borders features to your map\n", + "* BONUS: Use the metpy `add_timestamp` method from `metpy.plots` to add a timestamp to the plot.\n", + "* DAILY DOUBLE: Using the `start_date_time` attribute on the dataset `ds`, change the call to `add_timestamp` to use that date and time and the pretext to say `GOES 16 Channel X`." ] }, { @@ -478,18 +460,27 @@ "metadata": {}, "outputs": [], "source": [ + "import cartopy.feature as cfeature\n", "from metpy.plots import add_timestamp\n", "\n", - "text_artist = add_timestamp(ax, timestamp, pretext='GOES-16 Ch.{} '.format(channel))\n", + "fig = plt.figure(figsize=(10, 8))\n", + "ax = fig.add_subplot(1, 1, 1, projection=proj)\n", "\n", - "fig" + "# Plot the data using imshow\n", + "# YOUR CODE GOES HERE\n", + "\n", + "# Add country borders and states (use your favorite linestyle!)\n", + "# YOUR CODE GOES HERE\n", + "\n", + "# Bonus/Daily Double\n", + "# YOUR CODE GOES HERE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "That works, but maybe we prefer the look of outlined fonts inside the image. This can be done using matplotlib's [path effects](http://matplotlib.org/users/patheffects_guide.html). More details on path effects can be found in the [documentation](http://matplotlib.org/users/patheffects_guide.html) if you want to know more. Luckily, we just need to set the `high_contrast` arguement to `True` for this to work and adjust the y position to be inside the plot." + "#### Solution" ] }, { @@ -498,16 +489,24 @@ "metadata": {}, "outputs": [], "source": [ - "add_timestamp(ax, timestamp, pretext='GOES-16 Ch.{} '.format(channel),\n", - " high_contrast=True, fontsize=16, y=0.01)\n", - "fig" + "# %load solutions/sat_map.py\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now we have two timestamps! Let's remove the old one (this is why we kept a handle to it)." + "### Using Colortables\n", + "The map is much improved now, but it would look much better with a different color scheme. \n", + "\n", + "Colormapping in matplotlib (which backs CartoPy) is handled through two pieces:\n", + "\n", + "- The norm (normalization) controls how data values are converted to floating point values in the range [0, 1]\n", + "- The colormap controls how values are converted from floating point values in the range [0, 1] to colors (think colortable)\n", + "\n", + "Let's start by setting the colormap to be black and white and normalizing the data to get the best contrast. We'll make a histogram to see the distribution of values in the data, then clip that range down to enhance contrast in the data visualization. **Note:** `cmap` and `norm` can also be set during the `imshow` call as keyword arguments.\n", + "\n", + "We use `compressed` to remove any masked elements before making our histogram." ] }, { @@ -516,15 +515,30 @@ "metadata": {}, "outputs": [], "source": [ - "text_artist.remove()\n", - "fig" + "plt.hist(dat.to_masked_array().compressed(), bins=255);" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "Finally, we'll add the MetPy (or Unidata) logo to the plot for a nice finshing touch." + "fig = plt.figure(figsize=(10, 8))\n", + "ax = fig.add_subplot(1, 1, 1, projection=proj)\n", + "\n", + "# Plot the data using imshow\n", + "im = ax.imshow(dat, origin='upper', \n", + " extent=(x.min(), x.max(), y.min(), y.max()))\n", + "\n", + "# Add country borders and states (use your favorite linestyle!)\n", + "ax.add_feature(cfeature.BORDERS, linewidth=2, edgecolor='black')\n", + "ax.add_feature(cfeature.STATES, linestyle=':', edgecolor='black')\n", + "\n", + "# Bonus/Daily Double\n", + "timestamp = datetime.strptime(ds.start_date_time, '%Y%j%H%M%S')\n", + "add_timestamp(ax, timestamp, pretext='GOES-16 Ch.{} '.format(channel),\n", + " high_contrast=True, fontsize=16, y=0.01)" ] }, { @@ -533,8 +547,13 @@ "metadata": {}, "outputs": [], "source": [ - "from metpy.plots import add_metpy_logo, add_unidata_logo\n", - "add_metpy_logo(fig)\n", + "# Set colormap\n", + "im.set_cmap('Greys')\n", + "\n", + "# Set norm\n", + "im.set_norm(plt.Normalize(200,255))\n", + "\n", + "# Show figure again\n", "fig" ] }, @@ -542,12 +561,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "
\n", - " EXERCISE:\n", - "
    \n", - "
  • In one cell write all necessary code to make a plot of the CONUS sector using a channel of your choice. You may plot data from any time on September 9, 2017. (Re-writing the import statements is not necessary.)
    • \n", - "
\n", - "
" + "In meteorology, we have many ‘standard’ colortables that have been used for certain types of data. We have included these in Metpy in the `metpy.plots.ctables` module. By importing the `ColortableRegistry` we gain access to the colortables, as well as handy normalizations to go with them. We can see the colortables available by looking at the dictionary keys. The colortables ending in `_r` are the reversed versions of the named colortable.\n", + "\n", + "Let’s use the `WVCIMSS` colormap, a direct conversion of the GEMPAK colormap. The code below asks for the colormap, as well as a normalization that covers the range 195 to 265. This was empirically determined to closely match other displays of water vapor data. We then apply it to the existing image we have been working with." ] }, { @@ -556,45 +572,26 @@ "metadata": {}, "outputs": [], "source": [ - "# Your code goes here\n" + "from metpy.plots import colortables" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "\n", - "
\n", - "
\n",
-    "cat = TDSCatalog('http://thredds-test.unidata.ucar.edu/thredds/catalog/casestudies/irma/goes16/CONUS/Channel10/20170909/catalog.xml')\n",
-    "ds = cat.datasets.filter_time_nearest(datetime(2017, 9, 9, 6))\n",
-    "print(ds.name)\n",
-    "ds = ds.remote_access(use_xarray=True)\n",
-    "timestamp = datetime.strptime(ds.start_date_time, '%Y%j%H%M%S')\n",
-    "dat = ds.metpy.parse_cf('Sectorized_CMI')\n",
-    "\n",
-    "x = dat['x']\n",
-    "y = dat['y']\n",
-    "\n",
     "wv_norm, wv_cmap = colortables.get_with_range('WVCIMSS_r', 195, 265)\n",
-    "\n",
-    "fig = plt.figure(figsize=(10, 10))\n",
-    "ax = fig.add_subplot(1, 1, 1, projection=dat.metpy.cartopy_crs)\n",
-    "\n",
-    "im = ax.imshow(dat, extent=(x.min(), x.max(), y.min(), y.max()), origin='upper',\n",
-    "               cmap=wv_cmap, norm=wv_norm)\n",
-    "ax.coastlines(resolution='50m', color='black')\n",
-    "ax.add_feature(cfeature.STATES, linestyle=':')\n",
-    "ax.add_feature(cfeature.BORDERS, linewidth=2)\n",
-    "
\n", - "
" + "im.set_cmap(wv_cmap)\n", + "im.set_norm(wv_norm)\n", + "fig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Top\n", + "Top\n", "
" ] }, @@ -648,9 +645,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "import matplotlib as mpl\n", @@ -679,9 +674,6 @@ "ax.coastlines(resolution='50m', color='black')\n", "ax.add_feature(cfeature.BORDERS, linewidth=2)\n", "\n", - "# Add the MetPy Logo\n", - "add_metpy_logo(fig, x=80, y=15)\n", - "\n", "# Loop over the datasets and make the animation\n", "for ds in datasets[::6]:\n", "\n", @@ -719,7 +711,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Top\n", + "Top\n", "
" ] } @@ -740,7 +732,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.6" + "version": "3.7.1" } }, "nbformat": 4, diff --git a/notebooks/Satellite_Data/solutions/data_url.py b/notebooks/Satellite_Data/solutions/data_url.py new file mode 100644 index 00000000..c665bde4 --- /dev/null +++ b/notebooks/Satellite_Data/solutions/data_url.py @@ -0,0 +1,19 @@ +from datetime import datetime + +# Create variables for URL generation +image_date = datetime.utcnow().date() +region = 'Mesoscale-1' +channel = 8 + +# We want to match something like: +# https://thredds-test.unidata.ucar.edu/thredds/catalog/satellite/goes16/GOES16/Mesoscale-1/Channel08/20181113/catalog.html + +# Construct the data_url string +data_url = ('https://thredds-test.unidata.ucar.edu/thredds/catalog/satellite' + '/goes16/GOES16/{region}/Channel{channel:02d}/' + '{dt:%Y%m%d}/catalog.xml'.format(region=region, + channel=channel, + dt=image_date)) + +# Print out your URL and verify it works! +print(data_url) diff --git a/notebooks/Satellite_Data/solutions/sat_map.py b/notebooks/Satellite_Data/solutions/sat_map.py new file mode 100644 index 00000000..9aaf1f87 --- /dev/null +++ b/notebooks/Satellite_Data/solutions/sat_map.py @@ -0,0 +1,15 @@ +fig = plt.figure(figsize=(10, 8)) +ax = fig.add_subplot(1, 1, 1, projection=proj) + +# Plot the data using imshow +ax.imshow(dat, origin='upper', + extent=(x.min(), x.max(), y.min(), y.max())) + +# Add country borders and states (use your favorite linestyle!) +ax.add_feature(cfeature.BORDERS, linewidth=2, edgecolor='black') +ax.add_feature(cfeature.STATES, linestyle=':', edgecolor='black') + +# Bonus/Daily Double +timestamp = datetime.strptime(ds.start_date_time, '%Y%j%H%M%S') +add_timestamp(ax, timestamp, pretext='GOES-16 Ch.{} '.format(channel), + high_contrast=True, fontsize=16, y=0.01) From 947458ec20df73cc5c1b5c4eb406ee8e6c76b639 Mon Sep 17 00:00:00 2001 From: "J.R. Leeman" Date: Mon, 26 Nov 2018 11:37:59 -0700 Subject: [PATCH 2/3] Add jupyterlab to environment.yml --- environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/environment.yml b/environment.yml index 43bc0e99..a36751f8 100644 --- a/environment.yml +++ b/environment.yml @@ -7,6 +7,7 @@ - matplotlib=2 - cartopy=0.16 - jupyter + - jupyterlab - netcdf4 - owslib - metpy From 1cf5da710041c7fd1cc1a14d95d0661d2ffe7b01 Mon Sep 17 00:00:00 2001 From: "J.R. Leeman" Date: Wed, 28 Nov 2018 08:07:26 -0700 Subject: [PATCH 3/3] Move to different radar for example --- notebooks/Siphon/Siphon Overview.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/notebooks/Siphon/Siphon Overview.ipynb b/notebooks/Siphon/Siphon Overview.ipynb index 2c1e8f9c..8e33beb2 100644 --- a/notebooks/Siphon/Siphon Overview.ipynb +++ b/notebooks/Siphon/Siphon Overview.ipynb @@ -85,7 +85,7 @@ "from siphon.catalog import TDSCatalog\n", "date = datetime.utcnow() - timedelta(days=1)\n", "cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/nexrad/level3/'\n", - " 'N0Q/DGX/{dt:%Y%m%d}/catalog.xml'.format(dt=date))" + " 'N0Q/LRX/{dt:%Y%m%d}/catalog.xml'.format(dt=date))" ] }, { @@ -336,7 +336,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.6" + "version": "3.7.1" } }, "nbformat": 4,