Model-Obs Difference Example

# Title: Model-Obs difference

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




use_mars = 0
if use_mars then
    # retrieve model data for 2m temperature (GRIB)
    t2m_fc48 = retrieve(
        type    : "fc",
        levtype : "sfc",
        param   : "2t",
        date    : -5,
        step    : 48,
        grid    : "o1280")

    # retrieve obs data (BUFR)
    synop = retrieve(
        type   : "ob",
        repres : "bu",
        date   : -3)
else
    t2m_fc48 = read('t2m_fc48.grib')
    synop = read('t2m_obs.bufr')
end if

# filter just the 2m temperature from the obs data (Geopoints)
synop_t2m = obsfilter(
    output    : "geopoints",
    parameter : "012004",
    data      : synop)

# compute the difference
diff = t2m_fc48 - synop_t2m

# define the graphical symbol plotting style
max_diff = maxvalue(abs(diff))

diff_symb = msymb(
    legend                               : "on",
    symbol_type                          : "marker",
    symbol_table_mode                    : "advanced",
    symbol_outline                       : "on",
    symbol_outline_colour                : "charcoal",
    symbol_advanced_table_selection_type : "list",
    symbol_advanced_table_level_list     : [-max_diff,-10,-5,-1,0,1,5,10,max_diff],
    symbol_advanced_table_colour_method  : "list",
    symbol_advanced_table_colour_list    : ["blue","sky","rgb(0.82,0.85,1)","white","white","rgb(0.9,0.8,0.8)","rgb(0.9,0.45,0.45)","red"],
    symbol_advanced_table_height_list    : [0.6,0.5,0.4,0.3,0.3,0.4,0.5,0.6]
    )

# shaded land and sea to make the points stand out more
grey_land_sea_shading = mcoast(
    map_coastline_land_shade        : "on",
    map_coastline_land_shade_colour : "RGB(0.89,0.89,0.89)",
    map_coastline_sea_shade         : "on",
    map_coastline_sea_shade_colour  : "grey",
    map_grid_latitude_increment     : 20,
    map_grid_longitude_increment    : 40,
    map_grid_colour                 : "charcoal"
    )

# use numpy to give us some statistics to put in the title
# - we could also have done this with Metview functions
vals = values(diff)
title = mtext(
    text_line_1 : 'Min = '  & minvalue(vals) &
                 ' Mean = ' & mean(vals) &
                 ' Max = '  & maxvalue(vals),
    text_font_size : 0.4)

# adjust the legend
legend = mlegend(legend_text_font_size : 0.35)

# set the view area
view = geoview(
    map_area_definition : "corners",
    area                : [30,-28,75,48],
    coastlines          : grey_land_sea_shading)

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

# plot the data with the style
plot(view, diff, diff_symb, title, legend)

# Title: Model-Obs difference

#  **************************** LICENSE START ***********************************
# 
#  Copyright 2018 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
import numpy as np

use_mars = 0
if use_mars:
    # retrieve model data for 2m temperature (GRIB)
    t2m_fc48 = mv.retrieve(
        type    = "fc",
        levtype = "sfc",
        param   = "2t",
        date    = -5,
        step    = 48,
        grid    = "o1280")

    # retrieve obs data (BUFR)
    synop = mv.retrieve(
        type   = "ob",
        repres = "bu",
        date   = -3)
else:
    t2m_fc48 = mv.read("t2m_fc48.grib")
    synop = mv.read("t2m_obs.bufr")


# filter just the 2m temperature from the obs data (Geopoints)
synop_t2m = mv.obsfilter(
    output    = "geopoints",
    parameter = "012004",
    data      = synop)

# compute the difference
diff = t2m_fc48 - synop_t2m

# define the graphical symbol plotting style
max_diff = mv.maxvalue(mv.abs(diff))

diff_symb = mv.msymb(
    legend                               = "on",
    symbol_type                          = "marker",
    symbol_table_mode                    = "advanced",
    symbol_outline                       = "on",
    symbol_outline_colour                = "charcoal",
    symbol_advanced_table_selection_type = "list",
    symbol_advanced_table_level_list     = [-max_diff,-10,-5,-1,0,1,5,10,max_diff],
    symbol_advanced_table_colour_method  = "list",
    symbol_advanced_table_colour_list    = ["blue","sky","rgb(0.82,0.85,1)","white","white","rgb(0.9,0.8,0.8)","rgb(0.9,0.45,0.45)","red"],
    symbol_advanced_table_height_list    = [0.6,0.5,0.4,0.3,0.3,0.4,0.5,0.6]
    )

# shaded land and sea to make the points stand out more
grey_land_sea_shading = mv.mcoast(
    map_coastline_land_shade        = "on",
    map_coastline_land_shade_colour = "RGB(0.89,0.89,0.89)",
    map_coastline_sea_shade         = "on",
    map_coastline_sea_shade_colour  = "grey",
    map_grid_latitude_increment     = 20,
    map_grid_longitude_increment    = 40,
    map_grid_colour                 = "charcoal"
    )

# use numpy to give us some statistics to put in the title
# - we could also have done this with Metview functions
vals = mv.values(diff)
title = mv.mtext(
    text_line_1 = 'Min = '  + format(np.amin(vals), '.2f') +
                 ' Mean = ' + format(np.mean(vals), '.2f') +
                 ' Max = '  + format(np.amax(vals), '.2f'),
    text_font_size = 0.4)

# adjust the legend
legend = mv.mlegend(legend_text_font_size = 0.35)

# set the view area
view = mv.geoview(
    map_area_definition = "corners",
    area                = [30,-28,75,48],
    coastlines          = grey_land_sea_shading)

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

# plot the data with the style
mv.plot(view, diff, diff_symb, title, legend)