Cross Section Going Through the South Pole 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 input (pressure level) data
use_mars = 0 # 0 or 1

# retrieve analysis data from MARS
if use_mars = 1 then
    f_o3 = retrieve(
        levelist : [1000,850,700,500,400,300,250,200,150,100,70,50,30,20,10,7,5,3,2,1],
        param    : "o3",
        date     : 20191015,
        time     : 12,
        area     : [-65,-180,-90,180],
        grid     : [1,1]
    )
# read data from file
else
    f_o3  = read("ozone_pl.grib")
end if

# ---------------------------
# Cross section
# ----------------------------

# define cross section line via the South Pole going 
# from [65S, 10E] to [65S, 170W].
line = [-65,10,-115,10]

# define vertical axis
vertical_axis = maxis(
    axis_tick_label_height: 0.4
    )

# define cross section view wind projected to the normal vector of the
# cross section plane
xs = mxsectview(
    bottom_level       : 100,
    top_level          : 1,
    line               : line,
    wind_parallel      : "off",
    wind_perpendicular : "on",
    wind_intenisty     : "off",
    vertical_scaling   : "log",
    vertical_axiix     : vertical_axis
    )

# define contour shading for ozone
o3_shade = mcont(
    legend                      : "on",
    contour_line_colour         : "charcoal",
    contour_highlight           : "off",
    contour_max_level           : 1.5E-5,
    contour_min_level           : 0,
    contour_shade               : "on",
    contour_shade_colour_method : "palette",
    contour_shade_method        : "area_fill",
    contour_shade_palette_name  : "m_blue_green2_10",
    grib_scaling_of_retrieved_fields: "off"
    )

# ---------------------------------------------
#  Map view with cross section line
# ---------------------------------------------

# define coastlines
coast = mcoast(
    map_coastline_thickness      : 3,
    map_grid_longitude_increment : 30,
    map_label                    : "off",
    map_coastline_land_shade      : "on",
    map_coastline_land_shade_colour: "grey"
    )
 
# define map view
map = geoview(
    map_projection         : "polar_stereographic",
    map_area_definition    : "corners",
    area                   : [-50.67,-148.88,-50.2,30.29],
    map_hemisphere         : "south",
    map_vertical_longitude : 160,
    coastlines          : coast,
    subpage_y_lenght : 75
    )

# cross section line visualiser
vis_line = input_visualiser(
    input_plot_type        : "geo_points",
    input_longitude_values : [line[2],line[4]],
    input_latitude_values  : [line[1],line[3]]
    )

graph_line = mgraph(
    graph_line_colour    : "red",
    graph_line_thickness : 4
    )

# ---------------------------------------------
#  Define layout
# --------------------------------------------

map_page = plot_page(
    bottom: 30,
    view   : map
    )

xs_page = plot_page(
    top  : 30,
    view : xs
    )

dw = plot_superpage(
    pages : [map_page, xs_page]
    )

# ---------------------------------------------
#  Generate plot
# --------------------------------------------

# xs legend
xs_legend = mlegend(legend_text_font_size: "0.35")

# xs title
xs_title = mtext(text_font_size: "0.45")


# define the output plot file
setoutput(pdf_output(output_name : 'cross_section_via_pole'))

# generate plot
plot(dw[1], vis_line, graph_line)
plot(dw[2], f_o3, o3_shade, xs_legend, xs_title)

"""
Cross Section Going Through the South Pole
============================================
"""

#  **************************** 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 input (pressure level) data
use_mars = False

# retrieve analysis data from MARS
if use_mars:
    f_o3 = mv.retrieve(
        levelist=[
            1000,
            850,
            700,
            500,
            400,
            300,
            250,
            200,
            150,
            100,
            70,
            50,
            30,
            20,
            10,
            7,
            5,
            3,
            2,
            1,
        ],
        param="o3",
        date=20191015,
        time=12,
        area=[-65, -180, -90, 180],
        grid=[1, 1],
    )
# read data from file
else:
    f_o3 = mv.read("ozone_pl.grib")

# ---------------------------
# Cross section
# ----------------------------

# define cross section line via the South Pole going
# from [65S, 10E] to [65S, 170W].
line = [-65, 10, -115, 10]

# define vertical axis
vertical_axis = mv.maxis(axis_tick_label_height=0.4)

# define cross section view wind projected to the normal vector of the
# cross section plane
xs = mv.mxsectview(
    bottom_level=100,
    top_level=1,
    line=line,
    wind_parallel="off",
    wind_perpendicular="on",
    wind_intenisty="off",
    vertical_scaling="log",
    vertical_axiix=vertical_axis,
)

# define contour shading for ozone
o3_shade = mv.mcont(
    legend="on",
    contour_line_colour="charcoal",
    contour_highlight="off",
    contour_max_level=1.5e-5,
    contour_min_level=0,
    contour_shade="on",
    contour_shade_colour_method="palette",
    contour_shade_method="area_fill",
    contour_shade_palette_name="m_blue_green2_10",
    grib_scaling_of_retrieved_fields="off",
)

# ---------------------------------------------
#  Map view with cross section line
# ---------------------------------------------

# define coastlines
coast = mv.mcoast(
    map_coastline_thickness=3,
    map_grid_longitude_increment=30,
    map_label="off",
    map_coastline_land_shade="on",
    map_coastline_land_shade_colour="grey",
)

# define map view
map = mv.geoview(
    map_projection="polar_stereographic",
    map_area_definition="corners",
    area=[-50.67, -148.88, -50.2, 30.29],
    map_hemisphere="south",
    map_vertical_longitude=160,
    coastlines=coast,
    subpage_y_lenght=75,
)

# cross section line visualiser
vis_line = mv.input_visualiser(
    input_plot_type="geo_points",
    input_longitude_values=[line[1], line[3]],
    input_latitude_values=[line[0], line[2]],
)

graph_line = mv.mgraph(graph_line_colour="red", graph_line_thickness=4)

# ---------------------------------------------
#  Define layout
# --------------------------------------------

map_page = mv.plot_page(bottom=30, view=map)

xs_page = mv.plot_page(top=30, view=xs)

dw = mv.plot_superpage(pages=[map_page, xs_page])

# ---------------------------------------------
#  Generate plot
# --------------------------------------------

# xs legend
xs_legend = mv.mlegend(legend_text_font_size="0.35")

# xs title
xs_title = mv.mtext(text_font_size="0.45")

# define the output plot file
mv.setoutput(mv.pdf_output(output_name = 'cross_section_via_pole'))

# generate plot
mv.plot(dw[0], vis_line, graph_line, 
        dw[1], f_o3, o3_shade, xs_legend, xs_title)