Units issue with isentropic_interpolation_as_dataset and CONUS404 data #3315
Comments
|
What is the output if you print the variable The |
|
Thanks for the quick response. The output of printing the pressure variable is as follows: |
|
Ok, the subtle problem here is that you manually are calling The long answer is that because So I would write your code above as: import metpy
import metpy.calc as mpcalc
from metpy.units import units
import xarray as xr
import xwrf
import numpy as np
#------------------------------------------------------------------------------
#--- Open the NetCDF file -----------------------------------------------------
filename = "/glade/work/kshourd/29_June_2012/wrfout_d01_2012-06-29_12:00:00"
#--- xWRF postprocess & grid destaggering
ds = xr.open_dataset(filename).xwrf.postprocess()
ds = ds.xwrf.destagger()
lat = ds["XLAT"].drop_vars("CLAT")
lon = ds["XLONG"].drop_vars("CLAT")
#--- Extract variables
z = ds['geopotential_height'].drop_vars("CLAT")
print(z)
ua = ds['wind_east'].drop_vars("CLAT")
va = ds['wind_north'].drop_vars("CLAT")
pressure = ds["air_pressure"].drop_vars("CLAT")
theta = ds["air_potential_temperature"].drop_vars("CLAT")
temp = mpcalc.temperature_from_potential_temperature(pressure, theta)
# Compute the geostrophic wind
print("Geostrophic wind calculation...")
ugeo, vgeo = mpcalc.geostrophic_wind(z)
#--- Isentropic interpolation w/ MetPy
thlv = np.arange(298,350,1)
isen_lev = thlv*units.kelvin
isentropic = mpcalc.isentropic_interpolation_as_dataset(isen_lev.squeeze(), temp.squeeze(), ugeo.squeeze(), vgeo.squeeze(), ua.squeeze(), va.squeeze()) |
dopplershift
added
Area: Units
|
Thanks for the tip! Unfortunately, I am still receiving the same error: |
|
Hrmmm. What does this show? print(temp.metpy.vertical)Also, any chance you could share the data file? That would make it easier to figure out what exactly is going on. |
|
Output from is as follows: As far as the file, it is quite large (2.54 GB) and I do not currently have it downloaded. I am accessing the file using NCAR GLADE. The files are publicly available here: https://rda.ucar.edu/datasets/ds559.0/dataaccess/ and I am using domain 01 from 29 June 2012 at 12UTC. If you would still like me to provide a file, please let me know and I can do so. |
|
I was helping with this over on xarray-contrib/xwrf#152, so I had the file downloaded. Here (should be) a quick dropbox link: https://www.dropbox.com/scl/fi/bs8bh22pj7nal0diwjlfg/wrfout_d01_2012-06-29_12-00-00?rlkey=bsjb8osu73heutbhdia7kj5gp&dl=0 |
|
@spacekace So the root problem here is the fact that the WRF output you're working with is on sigma vertical coordinates. Now, in theory you should be able to use isen_press, isen_ug, isen_vg, isen_ua, isen_va = mpcalc.isentropic_interpolation(isen_lev.squeeze(), pressure.squeeze(), temp.squeeze(),
ugeo.squeeze(), vgeo.squeeze(), ua.squeeze(), va.squeeze())Unfortunately, there's currently a bug due to some assumptions we make that prevent that from working. I'm hoping we can quickly fix that for inclusion in the 1.6 release that should be out "any day now". EDIT: A work-around in the meanwhile is to first interpolate to a fixed set of isobaric levels, though I'm not sure how good those results would be in comparison with interpolation on the native sigma levels. |
|
I have a PR ready to go to fix all this up, including adding a |
|
Thanks, all! Really appreciate the help here. I have also been trying to accomplish this with wrf-python to no avail. Looking forward to the 1.6 release! |
Internally we were already using 3D pressure, so just avoid needless broadcasting outside the 1D pressure case. This paves the way for direct interpolation from e.g. sigma coords. See Unidata#3315.
isentropic_interpolation_as_dataset was blindly assuming that the vertical coordinate was pressure, resulting in a confusing unit error when calling isentropic_interpolation. Check the dimensionality and raise a clearer ValueError when this is not the case. Also allow users to manually pass (as a keyword) an appropriate pressure array to use instead.
isentropic_interpolation_as_dataset was blindly assuming that the vertical coordinate was pressure, resulting in a confusing unit error when calling isentropic_interpolation. Check the dimensionality and raise a clearer ValueError when this is not the case. Also allow users to manually pass (as a keyword) an appropriate pressure array to use instead.
What went wrong?
I am trying to use metpy.calc.isentropic_interpolation_as_dataset to interpolate CONUS404 data to isentropic levels. The CONUS404 data is produced using WRF v3.9.1.1. As such, I am using xarray and xWRF to read and process the data for use with MetPy. Unfortunately, I am receiving an error from this function stating:
isentropic_interpolation given arguments with incorrect units: pressure requires "[pressure]" but given "dimensionless"Pressure is not an input for isentropic_interpolation_as_dataset and therefore I have no control over the pressure units referred to here. It looks like this function calls on the isentropic_interpolation function of metpy.calc, which does require pressure as an input, so perhaps however isentropic_interpolation_as_dataset is calculating pressure from the temperature values provided does not work with WRF/CONUS404 data. I have previously used isentropic_interpolation_as_dataset AND isentropic_interpolation without issue, but that was with ERA5 atmospheric reanalysis GRIB data, which is CF compliant and allows for the use of metpy.parse_cf.
Any insight into how I can fix this issue and complete the interpolation would be greatly appreciated. Thank you!!
Note: You will notice that the "CLAT" coordinate is dropped when reading in variables. This is because isentropic_interpolation_as_dataset would not use variables with two latitude dimensions (i.e., XLAT and CLAT).
Operating System
Linux
Version
1.5.1
Python Version
3.9.16
Code to Reproduce
Errors, Traceback, and Logs
The text was updated successfully, but these errors were encountered: