Converting/Interpolating Eta Model Levels to Height

[mariandob]

Hello,

I am new to using xWRF, and have been following the vertical interpolation example provided in the documentation to learn how to convert the native WRF sigma or Eta levels into a different coordinate system. Although the example provided in the documentation does this for the pressure coordinate, I would like to do this for the height (z) coordinate.

When I try to apply the same interpolation workflow (with minor adjustments) shown in the documentation for pressure, to this interpolation for height, I unfortunately find that the transformed variable (wind speed computed from MetPy in this case) is full of only NaN values. I have provided a minimal working example below for reference:

import xwrf
import xgcm
import numpy as np
from metpy.calc import wind_speed

#import sample data and run xWRF post-processing 
ds = xwrf.tutorial.open_dataset("wrfout", chunks='auto').xwrf.postprocess()

#destagger dataset
destaggered = ds.xwrf.destagger().metpy.quantify()

#Compute wind speed using metpy and add to destaggered dataset
destaggered['wind_speed'] = wind_speed(destaggered.U, destaggered.V)

#Define target levels
#In this case we want to convert eta/sigma coorindates to height (z)
target_levels = np.arange(0, 3000, 50) # in meters

#Define xgcm grid using destaggered dataset
grid = xgcm.Grid(destaggered, periodic=False)

#apply linear transformation using xgcm
wind_speed_transformed = grid.transform(destaggered.wind_speed.metpy.dequantify(), 'Z', target_levels, method='linear')

#print the first value of the transformed variable to see if we have any values
print(wind_speed_transformed.isel(Time=0, y=0, x=0, z=0).values)

nan 

Would anyone know what might be the piece I am missing here? Thank you for your time and attention, and for this wonderful package.

Best,

Marian

[lpilz]

Hey Marian,

thanks for the question! I think the problem with your approach is that you don't tell XGCM what to interpolate the height levels to. It needs the target_data kwarg containing the heights of the individual cells to do this properly, however the xwrf.tutorial example you used is so small that it doesn't actually contain any height information.

I tried to get an example working with the intake integration which we have at https://github.com/xarray-contrib/xwrf-data however it seems to me that it might be broken, since I'm getting KeyError: 'wrf-cmip6-noversioning/downscaled_products/gcm/cesm2_r11i1p1f1_ssp245/6hourly/2099/d02/wrfout_d01_2099-10-01_00'. @jthielen @mgrover1 Could you maybe have a look at this and check whether we have to update the .jsons? (And maybe add some CI so that this doesn't happen again?)

At any rate, this should work with any wrf output file, you just need the HGT, PH and PHB fields in there (xwrf computes geopotential_height from the latter two).

import xwrf
import xgcm
import numpy as np
from metpy.calc import wind_speed

#import sample data and run xWRF post-processing
import intake
import intake_xarray
import s3fs

# Read in the data catalog from github
cat = intake.open_catalog('https://raw.githubusercontent.com/xarray-contrib/xwrf-data/main/catalogs/catalog.yml')

# Read one of the sample datasets
data_var_list = ["U", "V", "geopotential_height", "HGT"]
ds = cat["xwrf-sample-ssp245"].to_dask().xwrf.postprocess()[data_var_list]

#destagger dataset
destaggered = ds.xwrf.destagger().metpy.quantify()

#Compute wind speed using metpy and add to destaggered dataset
destaggered['wind_speed'] = wind_speed(destaggered.U, destaggered.V)

#Define target levels
#In this case we want to convert eta/sigma coorindates to height (z)
target_levels = np.arange(0, 3000, 50) # in meters

#Define xgcm grid using destaggered dataset
grid = xgcm.Grid(destaggered, periodic=False)
destaggered['height_agl'] = (destaggered['geopotential_height']-destaggered['HGT']).metpy.convert_units('meters').metpy.dequantify()

#apply linear transformation using xgcm
wind_speed_transformed = grid.transform(destaggered.wind_speed, 'Z', target_levels, target_data=destaggered['height_agl'], method='linear')

#print the first value of the transformed variable to see if we have any values
print(wind_speed_transformed.isel(Time=0, y=0, x=0, height_agl=0).values)

[mariandob]

Thank you very much for this informative and concise answer, @lpilz. I think I can get things to work from here.

If I could ask a (non-urgent) follow-up question: Perhaps this is simply Meteorology 101, but what is the meaning of subtracting the geopotential height from the model terrain height (HGT)? Feel free to reference some text if you do not have the time to explain here.

Thank you again for the assistance. I appreciate it.

[lpilz]

No problem at all. The variable HGT stores the terrain height and subtracting it from the geopotential height yields the height above ground level (which is what I think you wanted to interpolate to).