This example is based on a notebook posted by Scott Collis. This notebook shows how to do the plotting part of Scott's example using cartopy.
First step is to load the required modules/objects. We'll be using the cartopy.crs
module to handle coordinate reference systems, and the feature module to plot the US states and boundaries from a shapefile:
import urllib2
import cartopy.crs as ccrs
from cartopy.feature import NaturalEarthFeature
import matplotlib.pyplot as plt
from netCDF4 import Dataset
This next function gets the path to a GOES IR data file (direct from Scott's notebook):
def fetch_gini(dir_string=('http://motherlode.ucar.edu:8080/'
'thredds/dodsC/satellite/IR/EAST-CONUS_4km/current/'),
pattern_match = 'EAST-CONUS_4km_IR_'):
dirlisting = urllib2.urlopen(dir_string)
all_lines = dirlisting.read().splitlines()
has_gini = []
for line in all_lines:
if '.gini' in line:
happy_place = line.find(pattern_match)
end_place = line.find('.gini')
my_upper_name = line[happy_place:end_place]+'.gini'
url = dir_string + my_upper_name
has_gini.append(url)
return has_gini[0]
Reading the data from OpenDAP is the same as in Scott's notebook. We just need to be careful that the IR data is read as the correct type and that the units of the x and y coordinates are converted from kilometers to meters:
url = fetch_gini()
print 'retreiving data from: {}'.format(url)
ds = Dataset(url)
ir = ds.variables['IR'][0].astype(np.uint8)
x = ds.variables['x'][:] * 1000
y = ds.variables['y'][:] * 1000
retreiving data from: http://motherlode.ucar.edu:8080/thredds/dodsC/satellite/IR/EAST-CONUS_4km/current/EAST-CONUS_4km_IR_20140126_1945.gini
We need to create a coordinate reference system (CRS) that describes the input data, we can do that with cartopy:
grid = ds.variables['LambertConformal']
lon0 = grid.longitude_of_central_meridian
lat0 = grid.latitude_of_projection_origin
data_crs = ccrs.LambertConformal(central_latitude=lat0, central_longitude=lon0,
secant_latitudes=(lat0, lat0))
For the plot we need to create a set of axes with the desired projection, in this case a Miller projection. Cartopy can do this for us using familiar matplotlib syntax. We'll then add some detail to the map and plot the IR data:
fig = plt.figure(figsize=(15, 11))
# Miller projection map over Lake Michigan:
ax = plt.axes(projection=ccrs.Miller())
ax.set_extent([-90, -85, 40, 45], crs=ccrs.Geodetic())
# Add state boundaries etc.:
us_states = NaturalEarthFeature('cultural',
'admin_1_states_provinces_lakes_shp',
'50m',
facecolor='none')
ax.add_feature(us_states)
# Draw some grid lines:
ax.gridlines()
# Plot the IR data, making sure to specify the data CRS:
pc = ax.pcolormesh(x, y, ir, transform=data_crs,
cmap=plt.get_cmap('gray'), vmin=105)
# Add a colorbar:
cb = plt.colorbar(pc, orientation='vertical')
Finally, let's plot again but with a zoomed out view to verify that we have defined our data CRS correctly. The corners of the GOES Lambert Conformal Conic grid are available at http://www.nws.noaa.gov/noaaport/html/icdtb48e.html, so we plot those points on top to visually verify the grid location:
fig = plt.figure(figsize=(15, 11))
# Orthographic projection:
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-95,
central_latitude=25))
ax.set_global()
# Draw coastlines on the map:
ax.coastlines()
# Draw some grid lines:
ax.gridlines()
# Plot the IR data, making sure to specify the data CRS:
pc = ax.pcolormesh(x, y, ir, transform=data_crs,
cmap=plt.get_cmap('gray'))
# Add a colorbar:
cb = plt.colorbar(pc, orientation='vertical')
# Plot corner points accoriding to http://www.nws.noaa.gov/noaaport/html/icdtb48e.html:
for point in [(-113.133, 16.369), (-65.091, 14.335),
(-49.385, 57.289), (-123.044, 59.844)]:
ax.plot(point[0], point[1], transform=ccrs.Geodetic(),
markersize=10, marker='o', markerfacecolor='y')