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) catalog. From there we can drill down to find satellite data products.
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. 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!
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.
image_date, region, and channel. Assign them to today's date, the mesoscale-1 region, and ABI channel 8. data_url from these variables and the URL we navigated to above.catalog.html to catalog.xml. What do you see?
In [ ]:
from datetime import datetime
# Create variables for URL generation
# YOUR CODE GOES HERE
# Construct the data_url string
# YOUR CODE GOES HERE
# Print out your URL and verify it works!
# YOUR CODE GOES HERE
In [ ]:
# %load solutions/data_url.py
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.
We can use 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.
We want to create a TDSCatalog object called cat that we can examine and use to get handles to work with the data.
In [ ]:
from siphon.catalog import TDSCatalog
In [ ]:
cat = TDSCatalog(data_url)
To find the latest file, we can look at the cat.datasets attribute. Let’s look at the first five datasets:
In [ ]:
cat.datasets[:5]
In [ ]:
cat.datasets[0]
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.
In [ ]:
dataset = cat.datasets[1]
In [ ]:
print(dataset)
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).
In [ ]:
ds = dataset.remote_access(use_xarray=True)
Now that we've got some data - let's see what we actually got our hands on.
In [ ]:
ds
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.
To plot our data we'll be using MetPy's new declarative plotting functionality. You can write lots of matplotlib based code, but this interface greatly reduces the number of lines you need to write to get a great starting plot and then lets you customize it. The declarative plotting interface consists of three fundamental objects/concepts:
ImagePlot, ContourPlot, or Plot2D.MapPanel is the only panel type available.So containers have panels which have plots. It takes a second to get that straight in your mind, but it makes setting up complex figures very simple.
For this plot we need a single panel and we want to plot the satellite image, so we'll use the ImagePlot.
In [ ]:
from metpy.plots import ImagePlot, MapPanel, PanelContainer
%matplotlib inline
Let's start out with the smalles element, the plot, and build up to the largest, the panel container.
First, we'll make the ImagePlot:
In [ ]:
img = ImagePlot()
img.data = ds
img.field = 'Sectorized_CMI'
Next, we'll make the panel that our image will go into, the MapPanel objecet and add the image to the plots on the panel.
In [ ]:
panel = MapPanel()
panel.plots = [img]
Finally, we make the PanelContainer and add the panel to it's panels. Remember that since we can have multiple plots on a panel and multiple panels on a plot, we use lists. In this case is just happens to be a list of length 1.
In [ ]:
pc = PanelContainer()
pc.panels = [panel]
Unlike working with matplotlib directly in the notebooks, this figure hasn't actually been rendered yet. Calling the show method of the panel container builds up everything, renders, and shows it to us.
In [ ]:
pc.show()
ImagePlot here and figure out how to set the colormap of the image. For this image, let's go with the WVCIMMS_r colormap as this is a mid-level water vapor image.add_timestamp method from metpy.plots to add a timestamp to the plot. You can get the axes object to plot on from the ImagePlot. The call will look something like img.ax. This needs to happen after the panels have been added to the PanelContainer.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.
In [ ]:
# Import for the bonus exercise
from metpy.plots import add_timestamp
# Make the image plot
# YOUR CODE GOES HERE
# Make the map panel and add the image to it
# YOUR CODE GOES HERE
# Make the panel container and add the panel to it
# YOUR CODE GOES HERE
# Show the plot
# YOUR CODE GOES HERE
In [ ]:
# %load solutions/sat_map.py
Colormapping in matplotlib (which backs CartoPy) is handled through two pieces:
Let's try to determine a good range of values for the normalization. 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.
We use compressed to remove any masked elements before making our histogram.
In [ ]:
import matplotlib.pyplot as plt
plt.hist(ds['Sectorized_CMI'].to_masked_array().compressed(), bins=255);
In [ ]:
# Make the image plot
img = ImagePlot()
img.data = ds
img.field = 'Sectorized_CMI'
img.colormap = 'WVCIMSS_r'
img.image_range = (200, 255)
# Make the map panel and add the image to it
panel = MapPanel()
panel.plots = [img]
# Make the panel container and add the panel to it
pc = PanelContainer()
pc.panels = [panel]
pc.size = (10, 8)
# Bonus
start_time = datetime.strptime(ds.start_date_time, '%Y%j%H%M%S')
add_timestamp(panel.ax, time=start_time)
# Show the plot
pc.show()
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. The WVCIMSS colormap is a direct conversion of a GEMPAK colormap. Let's remake that image with a range of 195 to 265. This was empirically determined to closely match other displays of water vapor data.
We'll also dress up the plot a bit more to show a few other features.
In [ ]:
# Make the image plot
img = ImagePlot()
img.data = ds
img.field = 'Sectorized_CMI'
img.colormap = 'WVCIMSS_r'
img.image_range = (195, 265)
# Make the map panel and add the image to it
panel = MapPanel()
panel.plots = [img]
panel.title = f'GOES-16 Ch.{channel}'
# Make the panel container and add the panel to it
pc = PanelContainer()
pc.panels = [panel]
pc.size = (10, 8)
# Bonus
start_time = datetime.strptime(ds.start_date_time, '%Y%j%H%M%S')
add_timestamp(panel.ax, time=start_time,
pretext='Mid-level Water Vapor ',
high_contrast=True, fontsize=16, y=0.01)
# Show the plot
pc.show()
NOTE: This is just a quick taste of producing an animation using matplotlib. The animation support in matplotlib is robust, but sometimes installation of the underlying tool (ffmpeg) can be a little tricky. In order to make sure we get don't get bogged down, this is really more of a demo than something expected to work out of the box.
Conda-forge has packages, so it may be as easy as:
In [ ]:
#!conda install -y -n unidata-workshop -c conda-forge ffmpeg
First, we'll import the animation support from matplotlib. We also tell it that we want it to render the animations to HTML using the HTML5 video tag:
In [ ]:
import os.path
import sys
from IPython.display import HTML
from matplotlib.animation import ArtistAnimation
import cartopy.feature as cfeature
import cartopy.crs as ccrs
from metpy.plots import colortables
We create the base figure, then we loop over a bunch of the datasets to create an animation. For each one we pull out the data and plot both the timestamp and the image. The ArtistAnimation class takes the Figure instance and a list as required arguments. The contents of this list are a collection of matplotlib artists for each frame of the animation. In the loop below, we populate this list with the Text instance created when adding the timestamp as well as the image that results from plotting the data.
In [ ]:
import matplotlib as mpl
mpl.rcParams['animation.embed_limit'] = 50
# List used to store the contents of all frames. Each item in the list is a tuple of
# (image, text)
artists = []
case_date = datetime(2017, 9, 9)
# Get the IRMA case study catalog
cat = TDSCatalog('http://thredds.ucar.edu/thredds/catalog/casestudies/irma'
f'/goes16/Mesoscale-1/Channel{channel:02d}/{case_date:%Y%m%d}/'
'catalog.xml')
datasets = cat.datasets.filter_time_range(datetime(2017, 9, 9), datetime(2017, 9, 9, 6))
# Grab the first dataset and make the figure using its projection information
ds = datasets[0]
ds = ds.remote_access(use_xarray=True)
dat = ds.metpy.parse_cf('Sectorized_CMI')
proj = dat.metpy.cartopy_crs
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1, 1, 1, projection=proj)
plt.subplots_adjust(left=0.005, bottom=0.005, right=0.995, top=0.995, wspace=0, hspace=0)
ax.coastlines(resolution='50m', color='black')
ax.add_feature(cfeature.BORDERS, linewidth=2)
wv_norm, wv_cmap = colortables.get_with_range('WVCIMSS_r', 195, 265)
# Loop over the datasets and make the animation
for ds in datasets[::-6]:
# Open the data
ds = ds.remote_access(use_xarray=True)
dat = ds.metpy.parse_cf('Sectorized_CMI')
# Pull out the image data, x and y coordinates, and the time. Also go ahead and
# convert the time to a python datetime
x = dat['x']
y = dat['y']
timestamp = datetime.strptime(ds.start_date_time, '%Y%j%H%M%S')
img_data = ds['Sectorized_CMI']
# Plot the image and the timestamp. We save the results of these plotting functions
# so that we can tell the animation that these two things should be drawn as one
# frame in the animation
im = ax.imshow(dat, extent=(x.min(), x.max(), y.min(), y.max()), origin='upper',
cmap=wv_cmap, norm=wv_norm)
text_time = add_timestamp(ax, timestamp, pretext=f'GOES-16 Ch.{channel} ',
high_contrast=True, fontsize=16, y=0.01)
# Stuff them in a tuple and add to the list of things to animate
artists.append((im, text_time))
# Create the animation--in addition to the required args, we also state that each
# frame should last 200 milliseconds
anim = ArtistAnimation(fig, artists, interval=200., blit=False)
anim.save('GOES_Animation.mp4')
HTML(anim.to_jshtml())