Skip to content

Commit

Permalink
Fixes to get docs passing (full update deferred to Unidata#1385)
Browse files Browse the repository at this point in the history
  • Loading branch information
@jthielen
jthielen committed Sep 23, 2020
1 parent 57bb1e7 commit 277512c
Show file tree
Hide file tree
Showing 4 changed files with 71 additions and 83 deletions.
32 changes: 13 additions & 19 deletions examples/cross_section.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,27 +50,21 @@
# For this example, we will be plotting potential temperature, relative humidity, and
# tangential/normal winds. And so, we need to calculate those, and add them to the dataset:

temperature, pressure, specific_humidity = xr.broadcast(cross['Temperature'],
cross['isobaric'],
cross['Specific_humidity'])

theta = mpcalc.potential_temperature(pressure, temperature)
rh = mpcalc.relative_humidity_from_specific_humidity(pressure, temperature, specific_humidity)

# These calculations return unit arrays, so put those back into DataArrays in our Dataset
cross['Potential_temperature'] = xr.DataArray(theta,
coords=temperature.coords,
dims=temperature.dims,
attrs={'units': theta.units})
cross['Relative_humidity'] = xr.DataArray(rh,
coords=specific_humidity.coords,
dims=specific_humidity.dims,
attrs={'units': rh.units})

cross['Potential_temperature'] = mpcalc.potential_temperature(
cross['isobaric'],
cross['Temperature']
)
cross['Relative_humidity'] = mpcalc.relative_humidity_from_specific_humidity(
cross['isobaric'],
cross['Temperature'],
cross['Specific_humidity']
)
cross['u_wind'] = cross['u_wind'].metpy.convert_units('knots')
cross['v_wind'] = cross['v_wind'].metpy.convert_units('knots')
cross['t_wind'], cross['n_wind'] = mpcalc.cross_section_components(cross['u_wind'],
cross['v_wind'])
cross['t_wind'], cross['n_wind'] = mpcalc.cross_section_components(
cross['u_wind'],
cross['v_wind']
)

print(cross)

Expand Down
103 changes: 48 additions & 55 deletions examples/isentropic_example.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,22 +37,9 @@

#############################
# We will reduce the dimensionality of the data as it is pulled in to remove an empty time
# dimension.
# dimension, as well as add longitude and latitude as coordinates (instead of data variables).

# Assign data to variable names
lat = data['lat']
lon = data['lon']
lev = data['isobaric']
times = data['time']

tmp = data['Temperature'][0]
uwnd = data['u_wind'][0]
vwnd = data['v_wind'][0]
spech = data['Specific_humidity'][0]

# pint doesn't understand gpm
data['Geopotential_height'].attrs['units'] = 'meter'
hgt = data['Geopotential_height'][0]
data = data.squeeze().set_coords(['lon', 'lat'])

#############################
# To properly interpolate to isentropic coordinates, the function must know the desired output
Expand All @@ -65,55 +52,55 @@
#
# Once three dimensional data in isobaric coordinates has been pulled and the desired
# isentropic levels created, the conversion to isentropic coordinates can begin. Data will be
# passed to the function as below. The function requires that isentropic levels, isobaric
# levels, and temperature be input. Any additional inputs (in this case relative humidity, u,
# and v wind components) will be linearly interpolated to isentropic space.

isent_ana = mpcalc.isentropic_interpolation(isentlevs,
lev,
tmp,
spech,
uwnd,
vwnd,
hgt,
temperature_out=True)
# passed to the function as below. The function requires that isentropic levels, as well as a
# DataArray of temperature on isobaric coordinates be input. Any additional inputs (in this
# case specific humidity, geopotential height, and u and v wind components) will be
# logarithmicaly interpolated to isentropic space.

isent_data = mpcalc.isentropic_interpolation_as_dataset(
isentlevs,
data['Temperature'],
data['u_wind'],
data['v_wind'],
data['Specific_humidity'],
data['Geopotential_height']
)

#####################################
# The output is a list, so now we will separate the variables to different names before
# plotting.
# The output is an xarray Dataset:

isentprs, isenttmp, isentspech, isentu, isentv, isenthgt = isent_ana
isentu.ito('kt')
isentv.ito('kt')
isent_data

########################################
# A quick look at the shape of these variables will show that the data is now in isentropic
# coordinates, with the number of vertical levels as specified above.
# Note that the units on our wind variables are not ideal for plotting. Instead, let us
# convert them to more appropriate values.

print(isentprs.shape)
print(isentspech.shape)
print(isentu.shape)
print(isentv.shape)
print(isenttmp.shape)
print(isenthgt.shape)
isent_data['u_wind'] = isent_data['u_wind'].metpy.convert_units('kt')
isent_data['v_wind'] = isent_data['v_wind'].metpy.convert_units('kt')

#################################
# **Converting to Relative Humidity**
#
# The NARR only gives specific humidity on isobaric vertical levels, so relative humidity will
# have to be calculated after the interpolation to isentropic space.

isentrh = 100 * mpcalc.relative_humidity_from_specific_humidity(isentprs, isenttmp, isentspech)
isent_data['Relative_humidity'] = mpcalc.relative_humidity_from_specific_humidity(
isent_data['pressure'],
isent_data['temperature'],
isent_data['Specific_humidity']
).metpy.convert_units('percent')

#######################################
# **Plotting the Isentropic Analysis**

# Set up our projection
# Set up our projection and coordinates
crs = ccrs.LambertConformal(central_longitude=-100.0, central_latitude=45.0)
lon = isent_data['pressure'].metpy.longitude
lat = isent_data['pressure'].metpy.latitude

# Coordinates to limit map area
bounds = [(-122., -75., 25., 50.)]
# Choose a level to plot, in this case 296 K
# Choose a level to plot, in this case 296 K (our sole level in this example)
level = 0

fig = plt.figure(figsize=(17., 12.))
Expand All @@ -125,26 +112,28 @@

# Plot the surface
clevisent = np.arange(0, 1000, 25)
cs = ax.contour(lon, lat, isentprs[level, :, :], clevisent,
colors='k', linewidths=1.0, linestyles='solid', transform=ccrs.PlateCarree())
cs = ax.contour(lon, lat, isent_data['pressure'].isel(isentropic_level=level),
clevisent, colors='k', linewidths=1.0, linestyles='solid',
transform=ccrs.PlateCarree())
cs.clabel(fontsize=10, inline=1, inline_spacing=7, fmt='%i', rightside_up=True,
use_clabeltext=True)

# Plot RH
cf = ax.contourf(lon, lat, isentrh[level, :, :], range(10, 106, 5),
cmap=plt.cm.gist_earth_r, transform=ccrs.PlateCarree())
cf = ax.contourf(lon, lat, isent_data['Relative_humidity'].isel(isentropic_level=level),
range(10, 106, 5), cmap=plt.cm.gist_earth_r, transform=ccrs.PlateCarree())
cb = fig.colorbar(cf, orientation='horizontal', aspect=65, shrink=0.5, pad=0.05,
extendrect='True')
cb.set_label('Relative Humidity', size='x-large')

# Plot wind barbs
ax.barbs(lon.values, lat.values, isentu[level, :, :].m, isentv[level, :, :].m, length=6,
ax.barbs(lon.values, lat.values, isent_data['u_wind'].isel(isentropic_level=level).values,
isent_data['v_wind'].isel(isentropic_level=level).values, length=6,
regrid_shape=20, transform=ccrs.PlateCarree())

# Make some titles
ax.set_title(f'{isentlevs[level]:~.0f} Isentropic Pressure (hPa), Wind (kt), '
'Relative Humidity (percent)', loc='left')
add_timestamp(ax, times[0].values.astype('datetime64[ms]').astype('O'),
add_timestamp(ax, isent_data['time'].values.astype('datetime64[ms]').astype('O'),
y=0.02, high_contrast=True)
fig.tight_layout()

Expand All @@ -157,7 +146,10 @@


# Calculate Montgomery Streamfunction and scale by 10^-2 for plotting
msf = mpcalc.montgomery_streamfunction(isenthgt, isenttmp) / 100.
msf = mpcalc.montgomery_streamfunction(
isent_data['Geopotential_height'],
isent_data['temperature']
).values / 100.

# Choose a level to plot, in this case 296 K
level = 0
Expand All @@ -177,20 +169,21 @@
use_clabeltext=True)

# Plot RH
cf = ax.contourf(lon, lat, isentrh[level, :, :], range(10, 106, 5),
cmap=plt.cm.gist_earth_r, transform=ccrs.PlateCarree())
cf = ax.contourf(lon, lat, isent_data['Relative_humidity'].isel(isentropic_level=level),
range(10, 106, 5), cmap=plt.cm.gist_earth_r, transform=ccrs.PlateCarree())
cb = fig.colorbar(cf, orientation='horizontal', aspect=65, shrink=0.5, pad=0.05,
extendrect='True')
cb.set_label('Relative Humidity', size='x-large')

# Plot wind barbs.
ax.barbs(lon.values, lat.values, isentu[level, :, :].m, isentv[level, :, :].m, length=6,
# Plot wind barbs
ax.barbs(lon.values, lat.values, isent_data['u_wind'].isel(isentropic_level=level).values,
isent_data['v_wind'].isel(isentropic_level=level).values, length=6,
regrid_shape=20, transform=ccrs.PlateCarree())

# Make some titles
ax.set_title(f'{isentlevs[level]:~.0f} Montgomery Streamfunction '
r'($10^{-2} m^2 s^{-2}$), Wind (kt), Relative Humidity (percent)', loc='left')
add_timestamp(ax, times[0].values.astype('datetime64[ms]').astype('O'),
add_timestamp(ax, isent_data['time'].values.astype('datetime64[ms]').astype('O'),
y=0.02, pretext='Valid: ', high_contrast=True)

fig.tight_layout()
Expand Down
10 changes: 5 additions & 5 deletions src/metpy/calc/thermo.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1824,8 +1824,8 @@ def most_unstable_parcel(pressure, temperature, dewpoint, height=None,


@exporter.export
@preprocess_and_wrap()
@add_vertical_dim_from_xarray
@preprocess_and_wrap()
@check_units('[temperature]', '[pressure]', '[temperature]')
def isentropic_interpolation(levels, pressure, temperature, *args, vertical_dim=0,
temperature_out=False, max_iters=50, eps=1e-6,
Expand Down Expand Up @@ -2047,7 +2047,7 @@ def isentropic_interpolation_as_dataset(
'isentropic_level': xr.DataArray(
levels.m,
dims=('isentropic_level',),
coords=levels.m,
coords={'isentropic_level': levels.m},
name='isentropic_level',
attrs={
'units': str(levels.units),
Expand All @@ -2056,7 +2056,7 @@ def isentropic_interpolation_as_dataset(
),
**{
key: value
for key, value in all_args[0].coords
for key, value in all_args[0].coords.items()
if key != vertical_dim
}
}
Expand All @@ -2074,11 +2074,11 @@ def isentropic_interpolation_as_dataset(
),
'temperature': (
new_dims,
ret[-1],
ret[1],
{'standard_name': 'air_temperature'}
),
**{
all_args[i].name: (new_dims, ret[i], all_args[i].attrs)
all_args[i].name: (new_dims, ret[i + 1], all_args[i].attrs)
for i in range(1, len(all_args))
}
},
Expand Down
9 changes: 5 additions & 4 deletions tutorials/xarray_tutorial.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,6 +164,9 @@
# Calculations
# ------------
#
# TODO: The below section is out of date as of PR #1490. Updates are needed as a part of the
# 1.0 Upgrade Guide and Doc Update (Issue #1385).
#
# Most of the calculations in `metpy.calc` will accept DataArrays by converting them
# into their corresponding unit arrays. While this may often work without any issues, we must
# keep in mind that because the calculations are working with unit arrays and not DataArrays:
Expand All @@ -175,10 +178,9 @@
# As an example, we calculate geostropic wind at 500 hPa below:

lat, lon = xr.broadcast(y, x)
f = mpcalc.coriolis_parameter(lat)
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat, initstring=data_crs.proj4_init)
heights = data['height'].metpy.loc[{'time': time[0], 'vertical': 500. * units.hPa}]
u_geo, v_geo = mpcalc.geostrophic_wind(heights, f, dx, dy)
u_geo, v_geo = mpcalc.geostrophic_wind(heights, dx, dy, lat)
print(u_geo)
print(v_geo)

Expand Down Expand Up @@ -213,9 +215,8 @@

heights = data['height'].metpy.loc[{'time': time[0], 'vertical': 500. * units.hPa}]
lat, lon = xr.broadcast(y, x)
f = mpcalc.coriolis_parameter(lat)
dx, dy = mpcalc.grid_deltas_from_dataarray(heights)
u_geo, v_geo = mpcalc.geostrophic_wind(heights, f, dx, dy)
u_geo, v_geo = mpcalc.geostrophic_wind(heights, dx, dy, lat)
print(u_geo)
print(v_geo)

Expand Down

0 comments on commit 277512c

Please sign in to comment.