Try this notebook in .

Computing wind shear from forecast fields and plotting it on a vertical cross section (GRIB)

In this example we will compute the vertical wind shear of horizontal wind speed from forecast fields stored in GRIB format and generate a vertical cross section out of it.

import metview as mv

File forecast_ml.grib contains ECMWF forecast fields (step=+24h) on model levels retrieved from the MARS archive. We read this data into a Fieldset which is Metview’s own class to handle GRIB data.

filename = "forecast_ml.grib"
if mv.exist(filename):
    fs = mv.read(filename)
else:
    fs = mv.gallery.load_dataset(filename)

Our dataset contains the following set of parameters:

set(fs.grib_get_string("shortName"))

{'lnsp', 'q', 't', 'u', 'v', 'z'}

Since geopotential is not archived on model levels in MARS our field z only contains the surface geopotential, which is - for technical reasons - stored on model level 1.

mv.read(data=fs,param="z").grib_get_string("levelist")

'1'

In this example we define the wind shear by dividing the wind speed difference with the distance of two adjacent model levels. Since we have neither wind speed nor geopotential in our dataset we need to compute both.

First, we compute the wind speed.

u = mv.read(data = fs, param = "u")
v = mv.read(data = fs, param = "v")
sp = mv.sqrt(u*u + v*v)

Next, we use mvl_geopotential_on_ml() to compute the geopotential on model levels. This will perform a hydrostatic vertical integration in each gridpoint.

t = mv.read(data = fs, param = "t")
q = mv.read(data = fs, param = "q")
zs = mv.read(data = fs, param = "z")
lnsp = mv.read(data = fs, param = "lnsp")

z = mv.mvl_geopotential_on_ml(t, q, lnsp, zs)

Finally, the actual windshear computation can be done in a single line using Metview’s high level fieldset operations. Here we made use of the fact that our fieldsets are sorted by level in ascending order (and we have 137 model levels in total).

shear = (sp[0:135] - sp[1:136]) / (z[0:135] - z[1:136])

We convert the result into m/s/km units for easier interpretation and plotting.

shear = shear * 9.81 * 1000

To plot our data, we need to tell Metview to send the graphical output to Jupyter.

mv.setoutput('jupyter')

Plotting is performed through Metview’s interface to the Magics library developed at ECMWF.

We will first define the cross section line and plot it onto a map together with the wind speed on model level 80 (~250 hPa) to make sure that the cross section will go through an interesting weather feature, the jet stream, in our case.

area_view = mv.geoview(
    map_area_definition = 'corners',
    area = [35,-20,80,40],
    coastlines = mv.mcoast(
        map_coastline_land_shade        = "on",
        map_coastline_land_shade_colour = "grey"
    )
)

line = [41,-2,78,32]

line_graph = mv.mgraph(
    graph_type         = "curve",
    graph_line_colour  = "coral",
    graph_line_thickness = 5
)

speed_cont = mv.mcont(legend= "on",
                contour_automatics_settings = "style_name",
                contour_style_name = "sh_grn_f30t100i10")

mv.plot(area_view,
        mv.read(data = sp, levelist = 80),speed_cont,
        mv.mvl_geoline(*line,1),line_graph)

In the cross section we will overlay the computed wind shear with the wind speed and also want to plot the orography along the section line. First, we define a cross section view with logarithmic vertical coordinates in pressure between 1015 and 80 hPa.

xs_view = mv.mxsectview(
  line = line,
  top_level = 80,
  vertical_scaling = "log"
)

Next, we define the contouring for the wind shear and wind speed, respectively. For the wind shear, being the main feature, we use isoline shading while the wind speed is simply plotted with dashed isolines. We also specify area filling to represent the orography.

sh_cont= mv.mcont(
    legend                     = "on",
    contour_line_colour        = "charcoal",
    contour_highlight           = "off",
    contour_max_level           = 40,
    contour_min_level           = -40,
    contour_level_count         = 21,
    contour_label               = "off",
    contour_shade               = "on",
    contour_shade_method        = "area_fill",
    contour_shade_colour_method = "palette",
    contour_shade_palette_name  = "eccharts_black_red_21"
)

sp_cont = mv.mcont(
    contour_line_style     = "dash",
    contour_line_thickness = 2,
    contour_line_colour    = "charcoal",
    contour_highlight      = "off",
    contour_min_level      = 20
)

orog_graph = mv.mgraph(
    graph_type         = "area",
    graph_shade_colour  = "charcoal"
)

The final step is to define the input datasets and pass them on to the cross section. Since our data is available on model levels and we want to use a pressure-based cross section we need to supply the lnsp filed (logarithm of surface pressure) so that Metview could perform the required vertical interpolation. The lnsp field is also required for plotting the orography.

xs_sh_data = mv.merge(shear,lnsp)
xs_sp_data = mv.merge(sp,lnsp)

mv.plot(xs_view,
        xs_sh_data,sh_cont,
        xs_sp_data,sp_cont,
        orog_graph)