Reference: http://scitools.org.uk/cartopy/docs/latest/gallery.html
# Importing libraries we will need.
import netCDF4
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
# Open the netCDF file
f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')
# Getting the n-dimensional data
tempv = f.variables['temperature']
depth = f.variables['Depth']
print("The temperature variable...\n")
# Note the temperature variable has time, depth, y and x dimensions
print(tempv)
print("The dimensions...\n")
print(tempv.dimensions)
# Continued from previous cell..
# Our goal is temperature as a function of depth so slicing along the depth axis
# at a specific time and place
temp = tempv[0, :, 123, 486]
# Masked arrays are arrays that have bad values idenitifed by the mask array.
print("The masked array containing the temperature data...")
print(repr(temp))
# trick for filtering out good values
x = temp[~temp.mask]
y = depth[~temp.mask]
print("Plotting...")
# plot and show data
plt.plot(y, x)
# close netCDF
f.close()
Peruse matplotlib gallery and see and emulate what you like.
# Adding text, adjusting borders, and figure size
# Get figure hook to manipulate our plot
fig = plt.figure()
desc ='Figure 1. Temperature as a function of ocean depth as\n'\
'predicted by the RTOFS model'
# Adding our description
plt.figtext(.5,.15,desc,fontsize=12,ha='center')
#adjust margin
fig.subplots_adjust(bottom=0.35)
#adjust figure size
fig.set_size_inches(7,7)
# Improve axes
# Get axis hook to manipulate our plot
ax = fig.add_subplot(111)
# Add axis labels
ax.set_xlabel('Depth (m)', fontweight='bold')
# \u00b0 : degree symbol
ax.set_ylabel(u'Temperature (\u00b0C)', fontweight='bold')
# Don't show top and right axis
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)
# Define ticks
ax.tick_params(axis='both', direction='out')
ax.get_xaxis().tick_bottom() # remove unneeded ticks
ax.get_yaxis().tick_left()
# Getting the data as we did before
f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')
tempv = f.variables['temperature']
depth = f.variables['Depth']
temp = tempv[0,:,123,486]
x = temp[~temp.mask] #trick for getting data
y = depth[~temp.mask]
# Plotting line with triangle markers, and red line color.
plt.plot(y, x, marker=r'^', color='r', markersize=10, clip_on=False, linewidth=2.0)
plt.show()
f.close()
# Importing CartoPy
import cartopy
# Works with matplotlib's built-in transform support.
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.Mercator())
# Sets the extent to cover the whole globe
ax.set_global()
# Adds standard background map
ax.stock_img()
Cartopy also has a lot of built-in support for a variety of map features:
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.PlateCarree())
# Add variety of features
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.COASTLINE)
# Can also supply matplotlib kwargs
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.add_feature(cartopy.feature.LAKES, alpha=0.5)
ax.add_feature(cartopy.feature.RIVERS)
# Set limits in lat/lon space
ax.set_extent([-140, -60, 25, 50])
You can also grab other features from the Natural Earth project: http://www.naturalearthdata.com/
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.PlateCarree())
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
# Grab state borders
state_borders = cartopy.feature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m',
facecolor='none')
ax.add_feature(state_borders, linestyle="--", edgecolor='blue')
# Set limits in lat/lon space
ax.set_extent([-140, -60, 25, 50])
Plotting data works with the set projection by specifying the transform for the data when plotting. CartoPy will take care of any needed interpolation and transformation. The Plate Carree projection (Equi-rectangular projection) specifies the data in lat/lon, but not on the earth sphere. The Geodetic projection uses lat/lon, but are considered on the sphere. This example shows the difference.
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.Robinson())
ax.set_global()
ax.stock_img()
ax.coastlines()
# Plot circle with and without specifying the transform
ax.plot(-105.295175, 40.013176, 'o')
ax.plot(-105.295175, 40.013176, 'ro', transform=cartopy.crs.Geodetic())
# Plot line between two lat/lon points. One using equi-rectangular and one geodetic
plt.plot([-0.08, 132.], [51.53, 43.17], 'r', transform=cartopy.crs.PlateCarree())
plt.plot([-0.08, 132.], [51.53, 43.17], 'g', transform=cartopy.crs.Geodetic())
Try it yourself: Let's look at those curved flight paths on a map
# Open and read netCDF variables
nc = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')
tempv = nc.variables['temperature']
lat = nc.variables['Latitude'][:]
lon = nc.variables['Longitude'][:]
data = tempv[0, 0, :, :]
# Set up a stereographic projection
proj = cartopy.crs.Stereographic(central_latitude=60., central_longitude=-120.)
# Construct figure
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(1, 1, 1, projection=proj)
# Setting the plot size and text
plt.figtext(.5, .15, 'Figure 1. Sea surface temperature as predicted by the RTOFS model', fontsize=12, ha='center')
# define color map
cmap = plt.cm.hsv
# Nice high-level, human-readable abstractions for dealing with maps.
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines(zorder=3)
ax.gridlines()
#Color-filled contour plot
cs = ax.contourf(lon, lat, data, 50, cmap=cmap, transform=cartopy.crs.PlateCarree(), zorder=2)
# Color bar
cbar = plt.colorbar(cs)
cbar.set_label('$^{o}C$')
nc.close()
from datetime import datetime
date = datetime(2014, 9, 15, 0) # date to plot.
# open dataset.
dataset = netCDF4.Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/OISST-V2-AVHRR_agg')
timevar = dataset.variables['time']
timeindex = netCDF4.date2index(date, timevar) # find time index for desired date.
# read sst. Will automatically create a masked array using
# missing_value variable attribute. 'squeeze out' singleton dimensions.
sst = dataset.variables['sst'][timeindex, :].squeeze()
# read lats and lons (representing centers of grid boxes).
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]
# create figure, axes instances.
fig = plt.figure(figsize=(20, 10))
proj = cartopy.crs.Robinson(central_longitude=-105.)
ax = fig.add_axes([0.05, 0.05, 0.9, 0.9], projection=proj)
# plot sst, then ice with imshow--this relies on evenly spaced points
bounds = (lons[0], lons[-1], lats[0], lats[-1])
im1 = ax.imshow(sst, extent=bounds, interpolation='nearest',
cmap=plt.cm.jet, transform=cartopy.crs.PlateCarree())
# Add gridlines and coastlines
ax.gridlines()
ax.coastlines('50m')
# add colorbar
plt.colorbar(im1)
# add a title.
ax.set_title('Analysis for %s' % date)
# specify date to plot.
yyyy=1993; mm=3; dd=14; hh=0
# set OpenDAP server URL.
URLbase = "http://nomads.ncdc.noaa.gov/thredds/dodsC/modeldata/cmd_pgbh/"
URL = URLbase + "{year}/{year}{month:02}/{year}{month:02}{day:02}/pgbh00.gdas.{year}{month:02}{day:02}{hour:02}.grb2".format(
year=yyyy, month=mm, day=dd, hour=hh)
data = netCDF4.Dataset(URL)
# read lats,lons
latitudes = data.variables['lat'][:]
longitudes = data.variables['lon'][:]
# get sea level pressure and 10-m wind data.
# mult slp by 0.01 to put in units of hPa.
slpin = 0.01 * data.variables['Pressure_msl'][:].squeeze()
uin = data.variables['U-component_of_wind_height_above_ground'][:].squeeze()
vin = data.variables['V-component_of_wind_height_above_ground'][:].squeeze()
# make 2-d grid of lons, lats
lons, lats = np.meshgrid(longitudes, latitudes)
# make orthographic basemap.
proj = cartopy.crs.Orthographic(central_longitude=-60., central_latitude=60.)
# create figure, add axes
fig = plt.figure(figsize=(20, 10))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=proj)
ax.set_global()
# set desired contour levels.
clevs = np.arange(960., 1061., 8.)
# plot SLP contours.
ax.contour(lons, lats, slpin, clevs, linewidths=0.5, colors='k', transform=cartopy.crs.PlateCarree())
ax.contourf(lons, lats, slpin, clevs, cmap=plt.cm.RdBu_r, transform=cartopy.crs.PlateCarree())
stride = 15
ax.quiver(lons[::stride, ::stride], lats[::stride, ::stride],
uin[::stride, ::stride], vin[::stride, ::stride],
transform=cartopy.crs.PlateCarree())
# draw coastlines and gridlines
ax.coastlines(linewidth=1.5)
ax.gridlines()
date = datetime(yyyy, mm, dd, hh)
ax.set_title('SLP and Wind Vectors ' + str(date))