Declarative Satellite Data

Unidata Logo

Plotting Satellite Data

Unidata Python Workshop


Example Satellite Image

Questions

  1. Where are current GOES data available?
  2. How can satellite data be obtained with Siphon?
  3. How can maps be made with the declarative plotting interface?

Table of Contents

  1. Finding GOES data
  2. Accessing data with Siphon
  3. Plotting the data
  4. How Subplots Work
  5. Adding Annotations

Finding GOES Data

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 Satellite Data directory, then GOES East Products and CloudAndMoistureImagery. 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.

EXERCISE:
  • Create variables named image_date, region, and channel. Assign them to today's date, the mesoscale-1 region, and ABI channel 8.
  • Construct a string data_url from these variables and the URL we navigated to above.
  • Verify that following your link will take you where you think it should.
  • Change the extension from catalog.html to catalog.xml. What do you see?
In [1]:
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
SOLUTION
In [2]:
# %load solutions/data_url.py

# Cell content replaced by load magic replacement.
from datetime import datetime

# Create variables for URL generation
image_date = datetime.utcnow().date()
region = 'CONUS'
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.ucar.edu/thredds/catalog/satellite/goes/east/products/'
            f'CloudAndMoistureImagery/{region}/Channel{channel:02d}/'
            f'{image_date:%Y%m%d}/catalog.xml')

# Print out your URL and verify it works!
print(data_url)
https://thredds.ucar.edu/thredds/catalog/satellite/goes/east/products/CloudAndMoistureImagery/CONUS/Channel08/20200117/catalog.xml

Top


Accessing data with Siphon

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 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, but siphon will issue a warning.

We want to create a TDSCatalog object called cat that we can examine and use to get handles to work with the data.

In [3]:
from siphon.catalog import TDSCatalog
In [4]:
image_date = datetime.utcnow().date()
region = 'CONUS'
channel = 8

data_url = ('https://thredds.ucar.edu/thredds/catalog/satellite/goes/east/products/'
            f'CloudAndMoistureImagery/{region}/Channel{channel:02d}/'
            f'{image_date:%Y%m%d}/catalog.xml')

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 [5]:
cat.datasets[:5]
Out[5]:
[OR_ABI-L2-CMIPC-M6C08_G16_s20200171511160_e20200171511160_c20200171511160.nc,
 OR_ABI-L2-CMIPC-M6C08_G16_s20200171506160_e20200171506160_c20200171506160.nc,
 OR_ABI-L2-CMIPC-M6C08_G16_s20200171501160_e20200171501160_c20200171501160.nc,
 OR_ABI-L2-CMIPC-M6C08_G16_s20200171456160_e20200171456160_c20200171456160.nc,
 OR_ABI-L2-CMIPC-M6C08_G16_s20200171451160_e20200171451160_c20200171451160.nc]
In [6]:
cat.datasets[0]
Out[6]:
OR_ABI-L2-CMIPC-M6C08_G16_s20200171511160_e20200171511160_c20200171511160.nc

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 [7]:
dataset = cat.datasets[1]
In [8]:
print(dataset)
OR_ABI-L2-CMIPC-M6C08_G16_s20200171506160_e20200171506160_c20200171506160.nc

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 [9]:
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 [10]:
ds
Out[10]:
<xarray.Dataset>
Dimensions:               (x: 2500, y: 1500)
Coordinates:
    time                  datetime64[ns] ...
  * y                     (y) float32 128212.0 128156.0 ... 44324.0 44268.0
  * x                     (x) float32 -101332.0 -101276.0 ... 38556.0 38612.0
Data variables:
    Sectorized_CMI        (y, x) float32 ...
    fixedgrid_projection  int32 ...
Attributes:
    title:                       Sectorized Cloud and Moisture Imagery for th...
    ICD_version:                 GROUND SEGMENT (GS) TO ADVANCED WEATHER INTE...
    Conventions:                 CF-1.6
    channel_id:                  8
    central_wavelength:          6.19
    abi_mode:                    6
    source_scene:                CONUS
    periodicity:                 5.0
    production_location:         WCDAS
    product_name:                ECONUS-020-B12-M6C08
    satellite_id:                GOES-16
    product_center_latitude:     30.08300267372796
    product_center_longitude:    -87.09695844824527
    projection:                  Fixed Grid
    bit_depth:                   12
    source_spatial_resolution:   2.0
    request_spatial_resolution:  2.0
    start_date_time:             2020017150616
    number_product_tiles:        15
    product_tile_width:          512
    product_tile_height:         512
    product_rows:                1500
    product_columns:             2500
    pixel_x_size:                2.0
    pixel_y_size:                2.0
    satellite_latitude:          0.0
    satellite_longitude:         -75.0
    satellite_altitude:          35786023.0
    created_by:                  ldm-alchemy
    product_tiles_received:      15

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. MetPy has tools to make parsing all of this much easier than it used to be, but when using the delcarative plotting syntax, everything happens "automagically" for you and you don't need to worry about coordinate systems!

Top


Plotting the Data

When using the declarative plotting syntax, we need to import XArray and the plotting tools we'll use from MetPy. If you're used to working with Matplotlib - here's how the declarative objects should map in your thinking:

Matplotlib Declarative Plotting Description
Figure Panel Container Place to draw one or more plots
Axes Panel (multiple kinds) A set of x-y places to draw data inside a figure/panel container
Plotting Methods Plotting Methods Methods that draw data on an axes/panel in various ways (line plot, contours, images, etc).

For satellite data, we'd use imshow in matplotlib, but with the declarative syntax we want to make an ImagePlot.

In [11]:
import xarray as xr

from metpy.plots.declarative import *

First, let's make the image plot. We create an instance of it and name it, then set the data we want to use (our XArray data array) and the field from that data array we want to plot (Sectorized_CMI in this case).

In [12]:
img = ImagePlot()
img.data = ds
img.field = 'Sectorized_CMI'

Now we want to build a panel for that data to go in. Since the data are georeferenced, we'll use a MapPanel.

In [13]:
panel = MapPanel()
panel.plots = [img]

Finally, we want to build a panel container (figure) for the panel(s) to go in. Here we just have one panel.

In [14]:
pc = PanelContainer()
pc.panels = [panel]

Show the results of our "hard" work!

In [15]:
pc.show()
EXERCISE:
  • Now it's your turn! Set the colormap attribute of the ImagePlot to an appropriate colortable from those in metpy.plots.colortables here.
  • Using the title attribute of the MapPanel, give the plot a meaningful title.
  • Set the size attribute on the PanelContainer to a tuple representing the image size in the format (x, y).
In [16]:
# Create the ImagePlot
# YOUR CODE GOES HERE

# Create the MapPanel
# YOUR CODE GOES HERE

# Create the PanelContainer
# YOUR CODE GOES HERE

# Show it!
# YOUR CODE GOES HERE
SOLUTION
In [17]:
# %load solutions/dec_map_modifications.py

# Cell content replaced by load magic replacement.
# Create the ImagePlot
img = ImagePlot()
img.data = ds
img.field = 'Sectorized_CMI'
img.colormap = 'WVCIMSS'

# Create the MapPanel
panel = MapPanel()
panel.plots = [img]
panel.title = "CONUS GOES-East"

# Create the PanelContainer
pc = PanelContainer()
pc.panels = [panel]
pc.size = (10, 8)

# Show it!
pc.show()

How Subplots Work

EXERCISE: Now that you're familiar with the basic workflow of getting some satellite data and creating a plot, this is a challenge exercise to extend and combine those ideas.
  • Create a two panel plot (one PanelContainer, two MapPanels, two ImagePlots) with band 7 (IR) and band 9 (Mid-Level Water Vapor).
  • Set the layout attribute on the left panel to (2, 1, 1) and on the right to (2, 1, 2)
  • Add both panels to the PanelContainer
  • Label the plots
  • Apply appropriate color maps
Remember: build from the smallest component (the ImagePlot) to the biggest (the PanelContainer). BONUS: Functionalize this to eliminate repeated code!
In [18]:
#
# Get the data
#

# Get the data for band 7
# YOUR CODE GOES HERE

# Get the data for band 10
# YOUR CODE GOES HERE

#
# Make the plot
#
# YOUR CODE GOES HERE
SOLUTION
In [19]:
# %load solutions/dec_map_multipanel.py

# Cell content replaced by load magic replacement.
#
# Get the data
#

# Get the data for band 7
url = ('https://thredds.ucar.edu/thredds/catalog/satellite/goes/east/products/'
            f'CloudAndMoistureImagery/{region}/Channel07/'
            f'{datetime.utcnow():%Y%m%d}/catalog.xml')
cat = TDSCatalog(url)
band_7 = cat.datasets[-2].remote_access(use_xarray=True)

# Get the data for band 10
url = ('https://thredds.ucar.edu/thredds/catalog/satellite/goes/east/products/'
            f'CloudAndMoistureImagery/{region}/Channel10/'
            f'{datetime.utcnow():%Y%m%d}/catalog.xml')
cat = TDSCatalog(url)
band_10 = cat.datasets[-2].remote_access(use_xarray=True)

#
# Make the plot
#
band_7_img = ImagePlot()
band_7_img.data = band_7
band_7_img.field = 'Sectorized_CMI'
band_7_img.colormap = 'ir_rgbv'

band_10_img = ImagePlot()
band_10_img.data = band_10
band_10_img.field = 'Sectorized_CMI'
band_10_img.colormap = 'WVCIMSS'

left_panel = MapPanel()
left_panel.plots = [band_7_img]
left_panel.title = 'IR - GOES EAST'
left_panel.layout = (1, 2, 1)

right_panel = MapPanel()
right_panel.plots = [band_10_img]
right_panel.title = 'Mid-Level Water Vapor - GOES EAST'
right_panel.layout = (1, 2, 2)

pc = PanelContainer()
pc.panels = [left_panel, right_panel]
pc.size = (12, 10)

pc.show()
SOLUTION
In [20]:
# %load solutions/dec_map_multipanel_bonus.py

# Cell content replaced by load magic replacement.
def make_channel_image_plot(channel):
    url = ('https://thredds.ucar.edu/thredds/catalog/satellite/goes/east/products/'
                f'CloudAndMoistureImagery/{region}/Channel{channel:02d}/'
                f'{datetime.utcnow():%Y%m%d}/catalog.xml')
    cat = TDSCatalog(url)
    band_data = cat.datasets[-2].remote_access(use_xarray=True)
    
    img = ImagePlot()
    img.data = band_data
    img.field = 'Sectorized_CMI'
    return img

#
# Make the plot
#
band_7_img = make_channel_image_plot(7)
band_10_img = make_channel_image_plot(10)

band_7_img.colormap = 'ir_rgbv'
band_10_img.colormap = 'WVCIMSS'

left_panel = MapPanel()
left_panel.plots = [band_7_img]
left_panel.title = 'IR - GOES EAST'
left_panel.layout = (1, 2, 1)

right_panel = MapPanel()
right_panel.plots = [band_10_img]
right_panel.title = 'Mid-Level Water Vapor - GOES EAST'
right_panel.layout = (1, 2, 2)

pc = PanelContainer()
pc.panels = [left_panel, right_panel]
pc.size = (12, 10)

pc.show()