Cross Section in Pressure with Orography and Boundary Layer Height Example

#Metview Macro

#  **************************** LICENSE START ***********************************
# 
#  Copyright 2020 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 ************************************
# 

# get data
use_mars = 0 # 0 or 1
if use_mars then
    # retrieve data from MARS
    ret_core = (
        date     : 20200723,
        time     : 0,
        area     : [-8,25,-25,55],
        grid     : [0.2,0.2]
    )
    
    # forecast fields on all the model levels (bottom=ML-137, top=ML-1)
    # Note: log surface pressure (lnsp) is defined on ML-1!
    fs_ml = retrieve(ret_core,
        type: "fc",
        levtype: "ml",
        levelist: [1, "TO", 137],
        step: 12,
        param: ["t", "q", "u", "v", "lnsp"])
     
    # surface geopotential is available in  
    # the analysis only! It is available on ML-1!     
    zs = retrieve(ret_core,
        type: "an",
        levtype: "ml",
        levelist: 1,
        param: "z")
        
    # boundary layer height forecast
    blh = retrieve(ret_core,
        type: "fc",
        levtype: "sfc",
        param: "blh",
        step: 12)
else
    # read data from GRIB file
    fs_in = read("xs_blh.grib")   
    fs_ml = read(data: fs_in, levtype: "ml")
    zs = read(data: fs_ml, param: "z")
    blh = read(data: fs_in, param: "blh") 
end if

# extract ml data
t = read(data: fs_ml, param: "t")
q = read(data: fs_ml, param: "q")
lnsp = read(data_ml: fs_ml, param: "lnsp")
u = read(data: fs_ml, param: "u")
v = read(data: fs_ml, param: "v")

# define cross section line
line = [-10,28,-21,52]

# -------------------------------------------
# Generate cross section data for wind speed
# -------------------------------------------

# compute wind speed and set its paramId 
sp = sqrt(u*u + v*v)
sp = grib_set_long(sp, ["paramId", 10])

# compute cross section data for sp 
# (this is a NetCDF object)
xs_sp = mcross_sect(data: sp & lnsp, line: line)

# -------------------------------------------
# Generate curve for BL height 
# -------------------------------------------

# compute geopotential on model levels
z = mvl_geopotential_on_ml(t, q, lnsp, zs)

# compute pressure on model levels
p = unipressure(lnsp)

# interpolate pressure to the height of the BL
p_blh = ml_to_hl(p, z, zs, blh, "ground", "linear")

# define a curve object (in hPa) for the pressure of BL height
p_blh_curve = xs_build_curve(xs_sp, p_blh/100, "red", "solid", 3)

# define shading for wind speed using a palette
sp_cont = mcont(
    legend                       : "on",
    contour_line_colour          : "charcoal",
    contour_highlight            : "off",
    contour_level_selection_type : "interval",
    contour_max_level            : 18,
    contour_min_level            : 0,
    contour_interval             : 2,
    contour_shade                : "on",
    contour_shade_colour_method  : "palette",
    contour_shade_method         : "area_fill",
    contour_shade_palette_name   : "m_purple_9"
    ) 
  
# define vertical axis
vertical_axis = maxis(
    axis_orientation        : "vertical",
    axis_type               : "position_list",
    axis_tick_position_list : [1000, 925, 850, 700, 600, 500],
    axis_tick_label_height: 0.4
    )
  
# define cross section in log pressure (hPa)       
xs_view = mxsectview(
  line : line,
  top_level : 500,
  bottom_level: 1030,
  vertical_scaling : "log",
  vertical_axis: vertical_axis
)

# define orography area shading
orog_graph = mgraph(
    graph_type          : "area",
    graph_shade_colour  : "charcoal"
)

# define legend
legend = mlegend(
    legend_text_font_size : 0.35
    )
    
# define title
title = mtext(
    text_font_size : 0.4
    )
        
# define the output plot file
setoutput(pdf_output(output_name : 'cross_section_orog_and_blh'))

# generate plot
plot(xs_view,
        xs_sp, sp_cont,
        orog_graph,
        p_blh_curve,
        legend, title)

"""
Cross Section with Orography and Boundary Layer Height
=======================================================
"""

#  **************************** LICENSE START ***********************************
#
#  Copyright 2020 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 metview as mv

# get data
use_mars = False
if use_mars:
    # get data from MARS
    ret_core = {"date": "20200723", "time": 0, "area": [-8, 25, -25, 55], "grid": [0.2, 0.2]}

    # forecast fields on all the model levels (bottom=ML-137, top=ML-1)
    # Note= log surface pressure (lnsp) is defined on ML-1!
    fs_ml = mv.retrieve(
        **ret_core,
        type="fc",
        levtype="ml",
        levelist=[1, "TO", 137],
        step=12,
        param=["t", "q", "u", "v", "lnsp"]
    )

    # surface geopotential is available in
    # the analysis only! It is available on ML-1!
    zs = mv.retrieve(**ret_core, type="an", levtype="ml", levelist=1, param="z")

    # boundary layer height forecast
    blh = mv.retrieve(**ret_core, type="fc", levtype="sfc", param="blh", step=12)
else:
    # read data from GRIB file
    fs_in = mv.read("xs_blh.grib")
    fs_ml = mv.read(data=fs_in, levtype="ml")
    zs = mv.read(data=fs_ml, param="z")
    blh = mv.read(data=fs_in, param="blh")

# extract ml data
t = mv.read(data=fs_ml, param="t")
q = mv.read(data=fs_ml, param="q")
lnsp = mv.read(data_ml=fs_ml, param="lnsp")
u = mv.read(data=fs_ml, param="u")
v = mv.read(data=fs_ml, param="v")

# define cross section line
line = [-10, 28, -21, 52]

# -------------------------------------------
# Generate cross section data for wind speed
# -------------------------------------------

# compute wind speed and set its paramId
sp = mv.sqrt(u * u + v * v)
sp = mv.grib_set_long(sp, ["paramId", 10])

# compute cross section data for sp
# (this is a NetCDF object)
sp_fs = sp
sp_fs.append(lnsp)
xs_sp = mv.mcross_sect(data=sp_fs, line=line)

# -------------------------------------------
# Generate curve for BL height
# -------------------------------------------

# compute geopotential on model levels
z = mv.mvl_geopotential_on_ml(t, q, lnsp, zs)

# compute pressure on model levels
p = mv.unipressure(lnsp)

# interpolate pressure to the height of the BL
p_blh = mv.ml_to_hl(p, z, zs, blh, "ground", "linear")

# define a curve object (in hPa) for the pressure of BL height
p_blh_curve = mv.xs_build_curve(xs_sp, p_blh/100, "red", "solid", 3)

# define shading for wind speed using a palette
sp_cont = mv.mcont(
    legend="on",
    contour_line_colour="charcoal",
    contour_highlight="off",
    contour_level_selection_type="interval",
    contour_max_level=18,
    contour_min_level=0,
    contour_interval=2,
    contour_shade="on",
    contour_shade_colour_method="palette",
    contour_shade_method="area_fill",
    contour_shade_palette_name="m_purple_9",
)

# define vertical axis
vertical_axis = mv.maxis(
    axis_orientation="vertical",
    axis_type="position_list",
    axis_tick_position_list=[1000, 925, 850, 700, 600, 500],
    axis_tick_label_height=0.4,
)

# define cross section in log pressure (hPa)
xs_view = mv.mxsectview(
    line=line,
    top_level=500,
    bottom_level=1030,
    vertical_scaling="log",
    vertical_axis=vertical_axis,
)

# define orography area shading
orog_graph = mv.mgraph(graph_type="area", graph_shade_colour="charcoal")

# define legend
legend = mv.mlegend(legend_text_font_size=0.35)

# define title
title = mv.mtext(text_font_size=0.4)

# define the output plot file
mv.setoutput(mv.pdf_output(output_name="cross_section_orog_and_blh"))

# generate plot
mv.plot(xs_view, xs_sp, sp_cont, orog_graph, p_blh_curve, legend, title)