RH at PBL height

[pgamez]

Hi,

I was trying to interpolate a certain variable at the PBL height. According to the docs, the way to proceed is:

https://wrf-python.readthedocs.io/en/latest/user_api/generated/wrf.interplevel.html

from netCDF4 import Dataset
from wrf import getvar, interplevel

wrfin = Dataset("wrfout_d02_2010-06-13_21:00:00")

rh = getvar(wrfin, "rh")
z = getvar(wrfin, "z")
pblh = getvar(wrfin, "PBLH")

rh_pblh = interplevel(rh, z, pblh)

AFAIK the variable z has units of "meters above mean sea level", whilst variable PBLH has units of "meters above ground". I suppose that both variables must share the same units, so... Could you check if I am missing something or otherwise there is a bug in the docs?

Thanks in advance for your clarification, regards

Pedro

[michaelavs]

Do you mind clarifying your question for me a bit? Are you confused about the units used in the variables, or are you getting an error when trying to do this function locally?

[pgamez]

Hi Michaela,

I suspect the docs have an error at https://wrf-python.readthedocs.io/en/latest/user_api/generated/wrf.interplevel.html (2nd example: "Interpolate Relative Humidity to Boundary Layer Heights")

In this example from the docs, the vertical level z is given as meters above msl, but the desired level PBLH is the thickness above ground level. I think that both variables must be given in the same reference system to be properly interpolated

I am asking for confirmation. Thanks!

Pedro

[pgamez]

(...indeed, I am confused about the units)

[pgamez]

I think that z should be computed by using the argument msl=False at the function getvar, so:

from netCDF4 import Dataset
from wrf import getvar, interplevel

wrfin = Dataset("wrfout_d02_2010-06-13_21:00:00")

rh = getvar(wrfin, "rh")
z = getvar(wrfin, "z", msl=False)  # Height AGL
pblh = getvar(wrfin, "PBLH")       # Height AGL

rh_pblh = interplevel(rh, z, pblh)

[michaelavs]

I see what you are saying and why this was maybe confusing in set up, I had to do a bit of research before seeing what is meant by each variable input in interplevel(). The variable z is our reference height. This variable will be in meters above sea level, which is an elevation. The variable PBLH is our desired level. This variable will be in meters above ground, which is altitude. While I see why it would make sense to set msl=False, the way that the PBLH variable is recorded will already see itself as meters above z for the area you are looking at. So it isn't necessary to set msl=False for this case as PBLH should already be recorded in a way that would circumvent that. If, however, you are finding there are drastically different outputs to your models between when you set msl=False and when you don't, that would be a different case and we would need to do some debugging and figure out why that would be the case.

[pgamez]

Dear Michaela

Thank you for your explanation. I see that both z and rh are 3D vars, so the sentence: "the way that the PBLH variable is recorded will already see itself as meters above z" is not clear for me...

I have made a test in my wrfouf file. This corresponds to a region with an average elevation of 3000 m. I have run this script to check the differences between calling getvar with msl=True or msl=False:

#!/usr/bin/env python3

from netCDF4 import Dataset
from wrf import getvar, interplevel
import numpy as np

wrfin = Dataset("wrfout_d03_2020-03-19_12:00:00")

rh = getvar(wrfin, "rh")
pblh = getvar(wrfin, "PBLH")

for msl in [False, True]:
    print('> MSL =', msl)
    z = getvar(wrfin, "z", msl=msl)
    rh_pblh = interplevel(rh, z, pblh)

    values = np.ma.fix_invalid(rh_pblh.values)

    print('N:', values.count(), '/', values.size)
    print('min:', values.min())
    print('max:', values.max())
    print('mean:', values.mean())
    print()

The output of my script shows that there are really big differences depending on the msl argument:

> MSL = False
N: 53280 / 53280
min: 29.981607
max: 100.0
mean: 77.959045

> MSL = True
N: 251 / 53280
min: 31.27118
max: 79.26576
mean: 70.36978647908367

Setting msl=True results in a lot of missing data (the output has only 251 valid values). It make sense since you are interpolating from 3000+ meters (the 3D z coordinates of rh) to the pbl height (2D height coordinates, lets say about 1000 m) so the function interplevel thinks that it is below the surface! (there is no rh data below 3000m)

[pgamez]

(Using WRF v3.9.1.1)

[michaelavs]

Ok, now I think I have a better understanding of what the problem is. I think you are correct in saying/using msl=False. When msl=False, you are looking at above ground levels (AGL) and those are the levels you are interpolating with, which is where my confusion was. When msl=True, this would 'allow' for interpolating to sea level which is most likely why you were getting below ground results. I'm going to investigate the wrfout file that is used on the wrf readthedocs website that you linked to earlier and see if there is a reason this method wasn't used earlier. In the mean time, if you would like to create a pull request for this update so that the team can go over these changes and potentially merge them, please feel free to do so.

[michaelavs]

I was wondering if you could tell me the lat/lon range your wrfout file is using and also the lat/lon it is centered on? I want to double check my assumption that this may be due to elevation of the location of interest and ensure we document correctly. The example wrf file on the redthedocs page is centered in the American plains and has an elevation of 213 m, so it is relatively close to sea level and doesn't have many differences between msl=True and msl=False.

[pgamez]

Hi Michaela,

Good, the pull request is open!

The domain is 240x222 with this projection:

{'MAP_PROJ': 3,
 'CEN_LAT': 5.284218,
 'CEN_LON': -74.18291,
 'TRUELAT1': 4.3,
 'TRUELAT2': 4.3,
 'MOAD_CEN_LAT': 4.2999954,
 'STAND_LON': -74.9,
 'POLE_LAT': 90.0,
 'POLE_LON': 0.0,
 'DX': 1000.0,
 'DY': 1000.0}

And the lat/lon range is:

• latitude: from 4.2909927 to 6.2758484
• longitude: from -75.260796 to -73.10502

Thank you very much, Pedro

[ruratotoye2021]

Dear all, I have wrf output data file, I would like to plot wind speed at a given height using python but I failed. Someone to help me to interpolate wind speed for example at 60m agl.

[michaelavs]

Hi @ruratotoye2021, You will most likely want to do something like the code found on the interp1d function page. The changes to that code I would recommend would look something like:

import numpy as np
from wrf import getvar, interp1d
from netCDF4 import Dataset

wrfnc = Dataset("wrfout_d02_2010-06-13_21:00:00")

# Get a 1D vertical column for height 
ht_1d = wrf.getvar(wrfnc, "z", units="m")

# Get a 1D vertical column for wind speed 
wspd_1d = wrf.getvar(wrfnc, "wspd", units="knots")

# Want the speeds (in decameters) at 0, 60m
levels = np.asarray([0., 60.])

# Get the 0m and 60m wind speed values
interp_vals = interp1d(ht_1d, wspd_1d, levels)

This should give you data that can then be plotted with matplotlib and cartopy.

[marcelooyaneder]

Hi, its possible to use interplevel() function to interpoleta rh at a certain pressure level ?,

p = getvar(ncfile, "pressure",timeidx=idx) z = getvar(ncfile, "z", units="dm",timeidx=idx) hr = getvar(ncfile, "rh",timeidx=idx)[0,:]

isoLevel=900 ht_isoLevel = interplevel(z, p, isoLevel) hr_isoLevel = interplevel(hr, p, isoLevel)

Thanks in advance.