In [3]:
# Import the necessary libraries
%matplotlib inline
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.units import units
from netCDF4 import num2date
import numpy as np
import xarray as xr
import scipy.ndimage as ndimage
In [7]:
ds = xr.open_dataset('hgt.199601.nc')
ds
/home/jon/dev/xarray/xarray/conventions.py:494: SerializationWarning: variable 'hgt' has multiple fill values {9.96921e+36, -9.96921e+36}, decoding all values to NaN.
use_cftime=use_cftime,
Out[7]:
<xarray.Dataset>
Dimensions: (level: 29, time: 248, x: 349, y: 277)
Coordinates:
* time (time) datetime64[ns] 1996-01-01 ... 1996-01-31T21:00:00
* level (level) float32 1000.0 975.0 950.0 ... 150.0 125.0 100.0
lat (y, x) float32 ...
lon (y, x) float32 ...
* y (y) float32 0.0 32463.0 64926.0 ... 8927325.0 8959788.0
* x (x) float32 0.0 32463.0 64926.0 ... 11264660.0 11297120.0
Data variables:
Lambert_Conformal int32 ...
hgt (time, level, y, x) float32 ...
Attributes:
Conventions: CF-1.2
centerlat: 50.0
centerlon: -107.0
comments:
institution: National Centers for Environmental Prediction
latcorners: [ 1.000001 0.897945 46.3544 46.63433 ]
loncorners: [-145.5 -68.32005 -2.569891 148.6418 ]
platform: Model
standardpar1: 50.0
standardpar2: 50.000001
title: 8x Daily NARR
history: created Fri Feb 12 16:59:48 MST 2016 by NOAA/ESRL/PSD
dataset_title: NCEP North American Regional Reanalysis (NARR)
references: https://www.esrl.noaa.gov/psd/data/gridded/data.narr.html
source: http://www.emc.ncep.noaa.gov/mmb/rreanl/index.html
In [8]:
# Make all longitudes positive for plotting purposes
ds['lon'].data[ds['lon'].data < 0] += 360
In [11]:
hgt = ds.metpy.parse_cf('hgt')
hgt
Out[11]:
<xarray.DataArray 'hgt' (time: 248, level: 29, y: 277, x: 349)>
[695272216 values with dtype=float32]
Coordinates:
* time (time) datetime64[ns] 1996-01-01 ... 1996-01-31T21:00:00
* level (level) float32 1000.0 975.0 950.0 925.0 ... 150.0 125.0 100.0
lat (y, x) float32 ...
lon (y, x) float32 214.5 214.685 214.8705 ... 357.35638 357.4301
* y (y) float32 0.0 32463.0 64926.0 ... 8894862.0 8927325.0 8959788.0
* x (x) float32 0.0 32463.0 64926.0 ... 11264660.0 11297120.0
crs object Projection: lambert_conformal_conic
Attributes:
GRIB_id: 7
GRIB_name: HGT
dataset: NARR 3-hourly
grid_mapping: Lambert_Conformal
level_desc: Pressure Levels
long_name: 3-hourly Geopotential Heights on Pressure Levels
parent_stat: Other
standard_name: geopotential_height
statistic: Individual Obs
units: m
valid_range: [-1500. 35800.]
var_desc: Geopotential height
actual_range: [ -452.56787 16661.324 ]
In [12]:
# Calculate grid parameters
dx, dy = mpcalc.lat_lon_grid_deltas(hgt.metpy.longitude, hgt.metpy.latitude)
f = mpcalc.coriolis_parameter(hgt.metpy.latitude).to_base_units()
In [16]:
hgt_500 = hgt.metpy.sel(vertical=500 * units.hPa, time='1996-01-25 06:00')
uwnd, vwnd = mpcalc.geostrophic_wind(hgt_500, f, dx, dy)
In [23]:
# Set projection of data and grab data for plotting state boundaries
data_proj = hgt.metpy.cartopy_crs
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lakes',
scale='50m',
facecolor='none')
Plotting in native projection¶
In [25]:
# Make figure and set title
fig = plt.figure(figsize=(17., 11.))
ax = plt.subplot(111, projection=data_proj)
plt.title('500-hPa heights (m) and geostrophic wind vectors, Jan 25, 1996 06z', {'fontsize':18})
# Set extent and plot map lines
ax.set_extent([-90., -60, 35., 55.], ccrs.PlateCarree())
ax.coastlines('50m', edgecolor='gray', linewidth=0.75)
ax.add_feature(states_provinces, edgecolor='gray', linewidth=0.5)
# Plot 500-hPa heights
clevs_500_hght = np.arange(0, 8000, 60)
cs = ax.contour(hgt.metpy.x, hgt.metpy.y, hgt_500, clevs_500_hght,
colors='black', linewidths=1.5, linestyles='solid')
plt.clabel(cs, fmt='%d')
# Plot Wind Barbs
ax.quiver(hgt.metpy.x[0::4], hgt.metpy.y[0::4],
uwnd[0::4, 0::4], vwnd[0::4, 0::4], scale=2e3,
width=0.0015, pivot='middle')
plt.show()
Plotting in PlateCarree¶
In [29]:
# Make figure and set title
fig = plt.figure(figsize=(17., 11.))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
plt.title('500-hPa heights (m) and geostrophic wind vectors, Jan 25, 1996 06z', {'fontsize':18})
# Set extent and plot map lines
ax.set_extent([-90., -60, 35., 55.], ccrs.PlateCarree())
ax.coastlines('50m', edgecolor='gray', linewidth=0.75)
ax.add_feature(states_provinces, edgecolor='gray', linewidth=0.5)
# Plot 500-hPa heights
clevs_500_hght = np.arange(0, 8000, 60)
cs = ax.contour(hgt.metpy.x, hgt.metpy.y, hgt_500, clevs_500_hght,
colors='black', linewidths=1.5, linestyles='solid',
transform=data_proj)
plt.clabel(cs, fmt='%d')
# Plot Wind Barbs
ax.quiver(hgt.metpy.x[0::4], hgt.metpy.y[0::4],
uwnd[0::4, 0::4].m, vwnd[0::4, 0::4].m, scale=2e3,
width=0.0015, pivot='middle',
transform=data_proj)
plt.show()
Where the error arises¶
In [31]:
# Make figure and set title
fig = plt.figure(figsize=(17., 11.))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
plt.title('500-hPa heights (m) and geostrophic wind vectors, Jan 25, 1996 06z', {'fontsize':18})
# Set extent and plot map lines
ax.set_extent([-90., -60, 35., 55.], ccrs.PlateCarree())
ax.coastlines('50m', edgecolor='gray', linewidth=0.75)
ax.add_feature(states_provinces, edgecolor='gray', linewidth=0.5)
# Plot 500-hPa heights
# Using lat/lon rather than native CRS is fine for a scalar
clevs_500_hght = np.arange(0, 8000, 60)
cs = ax.contour(hgt.metpy.longitude, hgt.metpy.latitude, hgt_500, clevs_500_hght,
colors='black', linewidths=1.5, linestyles='solid',
transform=ccrs.PlateCarree())
plt.clabel(cs, fmt='%d')
# Plot Wind Barbs
# Can't use lat/lon for vector when not on lat/lon CRS!
ax.quiver(hgt.metpy.longitude[0::4, 0::4], hgt.metpy.latitude[0::4, 0::4],
uwnd[0::4, 0::4].m, vwnd[0::4, 0::4].m, scale=2e3,
width=0.0015, pivot='middle',
transform=ccrs.PlateCarree())
plt.show()
In [ ]: