-
|
I don't have time at the moment to simplify my code to something reproducible, but given a variable this code: from metpy.units import units
from metpy.interpolate import interpolate_to_isosurface
pvu = 1e-6 * units('degK') / units('kg') * units('m2') / units('s')
theta_dyn_trop = interpolate_to_isosurface(ds['pv'], ds['theta'], 2*pvu)generates: I am not sure if this is a shortcoming of the MetPy function, or if I need to massage the xarray Dataset in some way. I don't see anything in the Ideas? |
Beta Was this translation helpful? Give feedback.
Replies: 3 comments 3 replies
-
|
Try using the following: from metpy.units import units
from metpy.interpolate import interpolate_to_isosurface
theta_dyn_trop = interpolate_to_isosurface(ds['pv'], ds['theta'], 2)The choice of level is not unit-aware, although it probably could/should be. I believe I would have to take git blame on that one. |
Beta Was this translation helpful? Give feedback.
-
|
This change gives me: |
Beta Was this translation helpful? Give feedback.
-
|
Well, that error is because you have units for the data variables, but not your threshold. Units can work fine, it's |
Beta Was this translation helpful? Give feedback.
-
|
If you cast PV and theta to be Numpy arrays (without units), then use 2e-6 for the level should work. But we definitely need to look at what is going on here to determine if any changes are needed to work with DataArrays...I haven't tried to interpolate with a DataArray and have just cast them to Numpy arrays in the last (because we used to just have arrays with units attached and I ignored the units!). |
Beta Was this translation helpful? Give feedback.
-
|
For posterity (and a future bug report perhaps), here is the complete code corresponding to my original attempt: from datetime import datetime, timedelta
import xarray as xr
from xarray.backends import NetCDF4DataStore
import metpy.calc as mpcalc
from metpy.interpolate import interpolate_to_isosurface
from metpy.units import units
from siphon.catalog import TDSCatalog
def get_nam212(init_time, valid_time):
"""Obtain NAM data on the 212 grid."""
ymd = init_time.strftime('%Y%m%d')
hr = init_time.strftime('%H')
filename = f'{ymd}_{hr}00.grib2'
ds_name = 'NAM_CONUS_40km_conduit_' + filename
cat_name = ('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/NAM/'
'CONUS_40km/conduit/' + ds_name + '/catalog.xml')
cat = TDSCatalog(cat_name)
ds = cat.datasets[ds_name]
ncss = ds.subset()
query = ncss.query()
query.time(valid_time).variables('all')
nc = ncss.get_data(query)
data = xr.open_dataset(NetCDF4DataStore(nc)).metpy.parse_cf()
return data
def extract_and_normalize(xda, name):
var = xda.squeeze()
var = var.rename({var.metpy.vertical.name: 'pres'})
var = var.rename(name).metpy.assign_latitude_longitude()
return var
init_time = datetime(2021, 12, 6, 12)
plot_time = init_time + timedelta(hours=6)
nam = get_nam212(init_time, plot_time)
tmpk = extract_and_normalize(nam['Temperature_isobaric'], 'tmpk')
urel = extract_and_normalize(nam['u-component_of_wind_isobaric'], 'urel')
vrel = extract_and_normalize(nam['v-component_of_wind_isobaric'], 'vrel')
ds = xr.merge((tmpk, urel, vrel))
ds['theta'] = mpcalc.potential_temperature(ds['pres'], ds['tmpk'])
ds['pv'] = mpcalc.potential_vorticity_baroclinic(ds['theta'], ds['pres'], ds['urel'], ds['vrel'])
pvu = 1e-6 * units('degK') / units('kg') * units('m2') / units('s')
theta_dyn_trop = interpolate_to_isosurface(ds['pv'], ds['theta'], 2*pvu)I will try dropping down to NumPy arrays next and see what happens... |
Beta Was this translation helpful? Give feedback.
-
|
Here's the working version: from datetime import datetime, timedelta
import xarray as xr
from xarray.backends import NetCDF4DataStore
import metpy.calc as mpcalc
from metpy.interpolate import interpolate_to_isosurface
from metpy.plots import FilledContourPlot, MapPanel, PanelContainer
from metpy.units import units
from siphon.catalog import TDSCatalog
def get_nam212(init_time, valid_time):
"""Obtain NAM data on the 212 grid."""
ymd = init_time.strftime('%Y%m%d')
hr = init_time.strftime('%H')
filename = f'{ymd}_{hr}00.grib2'
ds_name = 'NAM_CONUS_40km_conduit_' + filename
cat_name = ('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/NAM/'
'CONUS_40km/conduit/' + ds_name + '/catalog.xml')
cat = TDSCatalog(cat_name)
ds = cat.datasets[ds_name]
ncss = ds.subset()
query = ncss.query()
query.time(valid_time).variables('all')
nc = ncss.get_data(query)
data = xr.open_dataset(NetCDF4DataStore(nc)).metpy.parse_cf()
return data
def extract_and_normalize(xda, name):
var = xda.squeeze()
var = mpcalc.smooth_gaussian(var, 16)
var = var.rename({var.metpy.vertical.name: 'pres'})
var = var.rename(name).metpy.assign_latitude_longitude()
return var
init_time = datetime(2021, 12, 6, 12)
plot_time = init_time + timedelta(hours=6)
nam = get_nam212(init_time, plot_time)
tmpk = extract_and_normalize(nam['Temperature_isobaric'], 'tmpk')
urel = extract_and_normalize(nam['u-component_of_wind_isobaric'], 'urel')
vrel = extract_and_normalize(nam['v-component_of_wind_isobaric'], 'vrel')
ds = xr.merge((tmpk, urel, vrel), combine_attrs='drop_conflicts')
ds['theta'] = mpcalc.potential_temperature(ds['pres'], ds['tmpk'])
ds['pv'] = mpcalc.potential_vorticity_baroclinic(ds['theta'], ds['pres'], ds['urel'], ds['vrel'])
# Feature request: Allow this
#pvu = 1e-6 * units('degK') / units('kg') * units('m2') / units('s')
#theta_dt = interpolate_to_isosurface(ds['pv'], ds['theta'], 2*pvu)
# Workaround
pv = ds['pv'].values
theta = ds['theta'].values
theta_dyn_trop = interpolate_to_isosurface(pv, theta, 2e-6)
theta_dt = xr.DataArray(theta_dyn_trop, dims=['y', 'x'],
coords={'y': ds['theta'].y, 'x': ds['theta'].x, 'metpy_crs': ds['theta'].metpy_crs},
attrs=ds.attrs)
fcp = FilledContourPlot()
fcp.data = theta_dt
fcp.colorbar = 'horizontal'
mp = MapPanel()
mp.area = [-97, -60, 33, 49]
mp.layers = ['coastline', 'borders', 'states']
mp.plots = [fcp]
pc = PanelContainer()
pc.size = (12, 8)
pc.panels = [mp]
pc.show() |
Beta Was this translation helpful? Give feedback.
If you cast PV and theta to be Numpy arrays (without units), then use 2e-6 for the level should work. But we definitely need to look at what is going on here to determine if any changes are needed to work with DataArrays...I haven't tried to interpolate with a DataArray and have just cast them to Numpy arrays in the last (because we used to just have arrays with units attached and I ignored the units!).