diff --git a/examples/cross_section.py b/examples/cross_section.py index d3ce9144f27..10a4dce8747 100644 --- a/examples/cross_section.py +++ b/examples/cross_section.py @@ -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) diff --git a/examples/isentropic_example.py b/examples/isentropic_example.py index a14038191a8..fee7532ac66 100644 --- a/examples/isentropic_example.py +++ b/examples/isentropic_example.py @@ -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 @@ -65,37 +52,31 @@ # # 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** @@ -103,17 +84,23 @@ # 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.)) @@ -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() @@ -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 @@ -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() diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py index 1a083bc8e7c..3c7391d0237 100644 --- a/src/metpy/calc/thermo.py +++ b/src/metpy/calc/thermo.py @@ -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, @@ -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), @@ -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 } } @@ -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)) } }, diff --git a/tutorials/xarray_tutorial.py b/tutorials/xarray_tutorial.py index d64fcf7094d..7bef7aa307f 100644 --- a/tutorials/xarray_tutorial.py +++ b/tutorials/xarray_tutorial.py @@ -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: @@ -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) @@ -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)