Pressure levels variables

In this tutorial, we will guide you through the process of analyzing multiple output dataset from the WRF (Weather Research and Forecasting) model and creating Cross section graphs featuring specific humidity and wind variables using the scahpy package. To begin, let’s import the glob package to manage files and paths, as well as the scahpy package.

from scahpy import *
import glob

Step 1: Reading WRF Data

We initiate the process by listing all the files we want to read using the glob package and assigning them to the variable list_files.

list_files = sorted(glob.glob('/data/datos/COW/OUT_DIAG_WRF/wrfouts/wrfout_d01_*'))

Given the capability to specify excluded variables when reading netCDF files using the drop_variables argument (refer to xarray functions open_dataset and open_mfdataset), we utilize the _drop_vars function from the module in_out. This function takes the list of variables we require and generates a list containing all variables present in the output file, subsequently removing those we are not interested in (such as ‘P’, ‘PB’, ‘U’, ‘W’, ‘QVAPOR’). For this purpose, we use the first file from our list of files, assigning it to the variable dvars.

dvars = in_out._drop_vars(list_files[0], ['P', 'PB', 'U', 'W', 'QVAPOR'], model='wrf')

Subsequently, we utilize the read_wrf_multi function to selectively read the variables of interest. This function accepts the input path (file_name in this case), the list of variables to be excluded (dvars), any required time difference (e.g., '5 hours'), and the corresponding sign of the time difference (-1 for negative, 1 for positive). The outcome is an xarray.Dataset containing longitude, latitude, bottom_top, time, and the specified variables. Optionally, you can designate a save path to export the netCDF.

ds = in_out.read_wrf_multi(list_files, dvars, '5 hours', -1)

Step 2: Calculating Specific Humidity and Total pressure

In this step, we will utilize the met_diag module to calculate specific humidity (calc_qe), and total pressure (calc_pres). Whe have the option of the parameter elimthat can be set to True or False in order to remove some variables used to calculate the final variable.

ds_lvl = met_vars.calc_qe(ds, elim=True)
ds_lvl = met_vars.calc_pres(ds_lvl, elim=True)

By running these commands, we ensure that our dataset ds_lvl now contains calculated specific humidity and total pressure, ready for further analysis or visualization.

Step 3: Aggregating the Data

Now, we’ll aggregate the data to operate on a daily time scale instead of hourly. To achieve this, we’ll utilize the dmy_var function from the temp_scales module. This function takes an xarray.Dataset as input, where we specify the desired time scale (e.g., ‘1D’ for daily, ‘ME’ for monthly, ‘YE’ for yearly). Additionally, we can specify which variables should be aggregated by sum, average, or median by providing lists for each aggregation method.

dd = temp_scales.dmy_var(ds_lvl, tiempo='1D', accum=None, avg=['Presion','U','W','QE'], mediana=None)

By executing this code, we’ll have our data aggregated to a daily time scale, with variables averaged according to our specifications.

Step 4: Interpolation to vertical levels

We use the functionvert_levs from spatial_scalesmodule to interpolate the data to same pressure levels, when we do not specify the levels, by default the interpolation is to: 1000,975,950,925,900,850,800,700,600,500,400,300,200 hPa. The dataset dd should contain the total pressure and the variables we are interested in.

dd2=vert_levs(dd,['U','W','QE'],lvls=None)

Step 5: Plotting Precipitation Maps

Next, we’ll generate cross section plots with specific humidity contours and wind vectors using the cross_section_xz function from the map_plots module. This function takes the dataset with specific humidity, total pressure and wind components, humidity levels, exportation settings, output path, temporal scale (‘H’ for hourly, ‘D’ for daily, ‘M’ for monthly, ‘Y’ for yearly) and vector speed.

# Example usage
levs=[0,0.2,0.4,0.6,0.8,1,1.5,2,2.5,5,7.5,10,12,15,18]
cmaps=cmocean.tools.lighten(matplotlib.colormaps['rainbow'],0.90)# 1d
df=dd2.sel(lat=-5,method='nearest').sel(lon=slice(-90,-80),time=slice('2023-03-10','2023-03-13'))
df['QE']=df['QE']*1000
cross_section_xz(df,'QE',levs,cmaps,'QE',quiverkey_speed=8, output_path=None, freq='D',
                           save_maps=False)