...
Illustrations of model levels, model half levels and model layers
Figure 1. An illustration of IFS model levels, showing | Figure 2. An illustration of IFS model levels, model half-levels and model layers. The pressure |
Geopotential on model levels
...
We use a Python script to download the ERA5 data from the MARS catalogue using the CDS API. The procedure is:
- Copy the script below to a text editor on your computer
- Edit the date, type, step, time, grid and area in the script to meet your requirements
- Save the script (for example with the filename as 'get_data_geopotential_on_ml.py')
- Run the script
Code Block | ||||||||
---|---|---|---|---|---|---|---|---|
| ||||||||
#!/usr/bin/env python import cdsapi c = cdsapi.Client() # data download specifications: cls = "ea" # do not change expver = "1" # do not change levtype = "ml" # do not change stream = "oper" # do not change date = "2018-01-01" # date: Specify a single date as "2018-01-01" or a period as "2018-08-01/to/2018-01-31". For periods > 1 month see https://confluence.ecmwf.int/x/l7GqB tp = "an" # type: Use "an" (analysis) unless you have a particular reason to use "fc" (forecast). time = "00:00:00" # time: ERA5 data is hourly. Specify a single time as "00:00:00", or a range as "00:00:00/01:00:00/02:00:00" or "00:00:00/to/23:00:00/by/1". c.retrieve('reanalysis-era5-complete', { 'class' : cls, 'date' : date, 'expver' : expver, 'levelist': '1/2/3/4/5/6/7/8/9/10/11/12/13/14/15/16/17/18/19/20/21/22/23/24/25/26/27/28/29/30/31/32/33/34/35/36/37/38/39/40/41/42/43/44/45/46/47/48/49/50/51/52/53/54/55/56/57/58/59/60/61/62/63/64/65/66/67/68/69/70/71/72/73/74/75/76/77/78/79/80/81/82/83/84/85/86/87/88/89/90/91/92/93/94/95/96/97/98/99/100/101/102/103/104/105/106/107/108/109/110/111/112/113/114/115/116/117/118/119/120/121/122/123/124/125/126/127/128/129/130/131/132/133/134/135/136/137', # For each of the 137 model levels 'levtype' : 'ml', 'param' : '130/133', # Temperature (t) and specific humidity (q) 'stream' : stream, 'time' : time, 'type' : tp, 'grid' : [1.0, 1.0], # Latitude/longitude grid: east-west (longitude) and north-south resolution (latitude). Default: 0.25 x 0.25 'area' : area, #example: [60, -10, 50, 2], # North, West, South, East. Default: global }, 'tq_ml.grib') c.retrieve('reanalysis-era5-complete', { 'class' : cls, 'date' : date, 'expver' : expver, 'levelist': '1', # Geopotential (z) and Logarithm of surface pressure (lnsp) are 2D fields, archived as model level 1 'levtype' : levtype, 'param' : '129/152', # Geopotential (z) and Logarithm of surface pressure (lnsp) 'stream' : stream, 'time' : time, 'type' : tp, 'grid' : [1.0, 1.0], # Latitude/longitude grid: east-west (longitude) and north-south resolution (latitude). Default: 0.25 x 0.25 'area' : area, #example: [60, -10, 50, 2], # North, West, South, East. Default: global }, 'zlnsp_ml.grib') |
...
For users experienced in Metview, there is a built-in function called mvl_geopotential_on_ml.
Interpolation of variables from model levels to custom pressure levels
The procedure described below is to convert ERA5 model levels data to custom pressure levels data.
Input:
- variable(s) on model levels and related logarithm surface pressure in grib format
- list of custom pressure levels required for interpolation to.
Output: NetCDF file containing variable(s) at each custom pressure level
Prerequisites for interpolating variables on model levels to custom pressure levels
You will need:
- Python3
- The CDS API installed. For more details, please follow the instructions here.
- The ecCodes library to read and write data.
Step 1: Download input data
First the required ERA5 variable(s) on model levels data are downloaded. The suggested procedure is:
- Copy the script below to a text editor on your computer
- Edit the date, type, step, time and grid in the script to meet your requirements. Note the 'area' keyword can also be used. The output filename can be modified accordingly.
- Save the script (for example with the filename as 'get_data.py')
- Run the script i.e. python3 get_data.py
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
# **************************** LICENSE START ***********************************
#
# Copyright 2022 ECMWF. This software is distributed under the terms
# of the Apache License version 2.0. In applying this license, ECMWF does not
# waive the privileges and immunities granted to it by virtue of its status as
# an Intergovernmental Organization or submit itself to any jurisdiction.
#
# ***************************** LICENSE END ************************************
import cdsapi
c = cdsapi.Client()
c.retrieve('reanalysis-era5-complete', {
'class': 'ea',
'date': '2021-01-01',
'expver': '1',
'levelist':'1/to/137',
'levtype': 'ml',
'param': '130/152',
'step': '0',
'stream': 'oper',
'time': '00/to/06/by/1',
'type': 'an',
'grid': '1.0/1.0'
}, 'output_00_06_130_152_1x1.grib') |
Running the script produces a file in the current working directory called 'output_00_06_130_152_1x1.grib' (a GRIB file containing the ERA5 variables needed.).
Step 2: Interpolate variables on model levels to custom pressure levels
The suggested procedure to run the Python script to compute the conversion of the variable from model levels to the custom pressure level is:
- Copy the script below to a text editor
- Save the script as 'conversion_from_ml_to_pl.py'
- Run the script 'conversion_from_ml_to_pl.py' with the correct arguments, i.e. :
python3 conversion_from_ml_to_pl.py -p 70000 -o output.nc -i output_00_06_130_152_1x1.grib
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
# **************************** LICENSE START ***********************************
#
# Copyright 2022 ECMWF. This software is distributed under the terms
# of the Apache License version 2.0. In applying this license, ECMWF does not
# waive the privileges and immunities granted to it by virtue of its status as
# an Intergovernmental Organization or submit itself to any jurisdiction.
#
# ***************************** LICENSE END ************************************
import cfgrib
import xarray as xr
import numpy as np
from eccodes import *
import matplotlib.pyplot as plt
import argparse
import sys
import os
def parse_args():
''' Parse program arguments using ArgumentParser'''
parser = argparse.ArgumentParser(description ="Python tool to calculate the model level variable at a given pressure level and write data to a netCDF file")
parser.add_argument('-p', '--pressure', required=True, nargs='+',type=float,
help='Pressure levels (Pa) to calculate the variable')
parser.add_argument('-o', '--output', required=False, help='name of the output file (default "output.nc"')
parser.add_argument('-i', '--input', required=True, metavar='input.grib', type=str,
help=('grib file with required variable(s) on model level and surface pressure fields',
'the model levels'))
args = parser.parse_args()
if not args.output:
args.output = 'output.nc'
return args
def get_input_variable_list(fin):
f = open(fin)
var_list = []
while 1:
gid = codes_grib_new_from_file(f)
if gid is None:
break
keys = ('dataDate', 'dataTime', 'shortName')
for key in keys:
if key == 'shortName':
var_list.append(codes_get(gid, key))
codes_release(gid)
var_list_unique = list(set(var_list))
f.close()
if 'lnsp' not in var_list_unique:
print("Error - lnsp variable missing from input file -exiting")
sys.exit()
if len(var_list_unique) < 2:
print("Error - Data variable missing from input file -exiting")
sys.exit()
return var_list_unique
def check_requested_levels(plevs):
check_lev = True
if len(plevs) > 1:
error_msg = "Error - only specify 1 input pressure level to interpolate to"
else:
for lev in plevs:
if lev < 0 or lev > 110000 :
check_lev = False
error_msg = "Error - negative values and large positive values for pressure are not allowed -exiting"
if check_lev == False:
print(error_msg)
sys.exit()
return check_lev
def check_in_range(data_array,requested_levels):
amin = data_array.minimum()
amax = data_array.maximum()
print("min max ",amin,amax)
def vertical_interpolate(vcoord_data, interp_var, interp_levels):
"""A function to interpolate sounding data from each station to
every millibar. Assumes a log-linear relationship.
Input
-----
vcoord_data : A 1D array of vertical level values (e.g., pressure from a radiosonde)
interp_var : A 1D array of the variable to be interpolated to all pressure levels
vcoord_interp_levels : A 1D array containing veritcal levels to interpolate to
Return
------
interp_data : A 1D array that contains the interpolated variable on the interp_levels
"""
l_count = 0
for l in interp_levels:
if l < np.min(vcoord_data) or l > np.max(vcoord_data):
interp_data[l_count] = np.NAN
# Make vertical coordinate data and grid level log variables
lnp = np.log(vcoord_data)
lnp_intervals = [np.log(x) for x in interp_levels]
# Use numpy to interpolate from observed levels to grid levels
interp_data = np.interp(lnp_intervals, lnp, interp_var)
return interp_data[0]
def calculate_pressure_on_model_levels(ds_var,ds_lnsp):
# Get the number of model levels in the input variable
nlevs=ds_var.sizes['hybrid']
# Get the a and b coefficients from the pv array to calculate the model level pressure
pv_coeff = np.array(ds_var.GRIB_pv)
pv_coeff=pv_coeff.reshape(2,nlevs+1)
a_coeff=pv_coeff[0,:]
b_coeff=pv_coeff[1,:]
# get the surface pressure in hPa
sp = np.exp(ds_lnsp)
p_half=[]
for i in range(len(a_coeff)):
p_half.append(a_coeff[i] + b_coeff[i] * sp)
p_ml=[]
for hybrid in range(len(p_half) - 1):
p_ml.append((p_half[hybrid + 1] + p_half[hybrid]) / 2.0)
ds_p_ml = xr.concat(p_ml, 'hybrid')
return ds_p_ml
def plot_profile(var_ml,press_ml, var_int_press,var_int_plevs,tstep,lat,lon):
var_v= var_ml.sel(time = var_ml.time[tstep],longitude=lon, latitude=lat, method='nearest')
var_v_values = var_v.values
var_p= press_ml.sel(time = var_ml.time[tstep],longitude=lon, latitude=lat, method='nearest')
var_p_values = var_p.values
var_ip= var_int_press.sel(time = var_ml.time[tstep],longitude=lon, latitude=lat, method='nearest')
var_ip_values = var_ip.values
var_ip_p = var_ip.pressure
var_ip_p_values = var_ip_p.values
plt.axis([min(var_v_values), max(var_v_values), max(var_p_values), min(var_p_values)])
plt.plot(var_v_values,var_p_values, 'o', color = 'black')
plt.plot(var_ip_values,var_ip_p_values,'o', color = 'red')
plt.show()
return
def calculate_interpolated_pressure_field(data_var_on_ml, data_p_on_ml,plevs):
nlevs = len(data_var_on_ml.hybrid)
p_array = np.stack(data_p_on_ml, axis=2).flatten()
# Flatten the data array to enable faster processing
var_array = np.stack(data_var_on_ml, axis=2).flatten()
no_grid_points = int(len(var_array)/nlevs)
interpolated_var = np.empty((len(plevs), no_grid_points))
ds_shape = data_var_on_ml.shape
nlats_values = data_var_on_ml.coords['latitude']
nlons_values = data_var_on_ml.coords['longitude']
nlats = len(nlats_values)
nlons = len(nlons_values)
# Iterate over the data, selecting one vertical profile at a time
count = 0
profile_count = 0
interpolated_values=[]
for point in range(no_grid_points):
offset = count*nlevs
var_profile = var_array[offset:offset+nlevs]
p_profile = p_array[offset:offset+nlevs]
interpolated_values.append(vertical_interpolate(p_profile, var_profile, plevs))
profile_count += len(p_profile)
count = count + 1
interpolated_field=np.asarray(interpolated_values).reshape(len(plevs),nlats,nlons)
return interpolated_field
def check_data_cube(dc):
checks = True
for var_name in dc.variables:
if var_name in ['time','step','hybrid','latitude','longitude','valid_time']:
continue
if var_name == 'lnsp':
lnsp_dims = ['time','latitude','longitude']
if all(value in lnsp_dims for value in dc.variables[var_name].dims):
continue
else:
print("Not all required lnsp dimensions found -exiting ", dc.variables[var_name].dims)
checks = False
else:
var_dims = ['time','hybrid','latitude','longitude']
if all(value in var_dims for value in dc.variables[var_name].dims):
continue
else:
print("Not all required variable dimensions found -exiting ",dc.variables[var_name].dims)
checks = False
continue
return checks
def main():
'''Main function'''
print("-p <pressure level (Pa) > -o <output_file> -i <input grib file>")
print("e.g. to process a grib file containing 6 hours of lnsp and temperature data to the 500 hPa level:")
print("python3 script.py -o output_press.nc -p 50000 -i output_00_06_130_152_1x1.grib`n")
args = parse_args()
print('Arguments: %s' % ", ".join(
['%s: %s' % (k, v) for k, v in vars(args).items()]))
plevels = args.pressure
plevels.sort(reverse = True)
check_requested_levels(plevels)
input_fname = args.input
output_fname = args.output
if not os.path.isfile(input_fname):
print("Input file does not exist - exiting")
sys.exit()
variable_list = get_input_variable_list(input_fname)
# Create a data object to hold the input and derived data
data_cube = xr.merge(cfgrib.open_datasets(input_fname, backend_kwargs={'read_keys': ['pv']}), combine_attrs='override')
if not check_data_cube(data_cube):
sys.exit()
# Get the ln surface pressure
lnsp = data_cube['lnsp']
for var in variable_list:
if var == 'lnsp':
continue
else:
data_cube['pml']=data_cube[var].copy()
break
for var in variable_list:
if var == 'lnsp' :
continue
data_pressure_on_model_levels_list =[]
for time_step in range(len(data_cube[var].time)):
data_slice_var=data_cube[var][time_step,:,:,:]
data_slice_lnsp=data_cube['lnsp'][time_step,:,:]
# Get the pressure field on model levels for each timestep
data_cube['pml'][time_step,:,:,:] = calculate_pressure_on_model_levels(data_slice_var,data_slice_lnsp)
data_cube['pml'].attrs = {'units' : 'Pa','long_name':'pressure','standard_name':'air_pressure','positive':'down'}
all_interpolated_var_fields_list = []
for var in variable_list:
if var == 'lnsp' or var == 'pml':
continue
interpolated_var_field = data_cube[var].copy()
interpolated_var_field = interpolated_var_field[:,0:len(plevels),:,:]
interpolated_var_field = interpolated_var_field.rename({'hybrid':'pressure'})
interpolated_var_field['pressure'] = plevels
for time_step in range(len(data_cube[var].time)):
var_on_ml = data_cube[var][time_step,:,:,:]
p_on_ml = data_cube['pml'][time_step,:,:,:]
interpolated_var_field[time_step,:,:,:] = calculate_interpolated_pressure_field(var_on_ml,p_on_ml,plevels)
all_interpolated_var_fields_list.append(interpolated_var_field)
all_interpolated_var_fields = xr.merge(all_interpolated_var_fields_list)
all_interpolated_var_fields['pressure'].attrs = {'units' : 'Pa','long_name':'pressure','standard_name':'air_pressure','positive':'down'}
all_interpolated_var_fields.to_netcdf(output_fname)
# Write interpolated data variable to output filename
PLOT_DATA = False
if PLOT_DATA:
latitude = 45.0
longitude = 0
tstep =0
plot_profile(data_cube[var],data_cube['pml'],interpolated_var_field,plevels,tstep,latitude,longitude)
print("Finished interpolation of variables to pressure level")
if __name__ == '__main__':
main()
|
This produces a netCDF file called 'output.nc' in the current directory containing the interpolated data.
Geopotential height
In ERA5, and often in meteorology, heights (the height of the land and sea surface, or specific heights in the atmosphere) are not represented as geometric height, or altitude (in metres above the spheroid), but as geopotential height (in metres above the geoid, which is represented by the mean sea level in ERA5). Note, that ECMWF usually archive the geopotential (in m2/s2), not the geopotential height.
...