ecmwf / magics-python

Python interface to Magics meteorological plotting package
Apache License 2.0
49 stars 12 forks source link

Trying to generate a vertical cross section plot of a geopotential height anomaly #25

Closed winash12 closed 4 years ago

winash12 commented 4 years ago

I am trying to generate a vertical cross section of a geopotential height anomaly from NOAA reanalysis data. as shown in this link https://imgur.com/a/ASUVWFn

The data is organized as follows - the latitude is chosen from 25 to 40 degrees north. And the longitude is from 60 to 100 degrees east. So the pressure data is at 17 levels. And it is a netCDF data file that is having the variable name 'hgt'. What I want is a plot of longitude by height and by that I mean - I have a longitude range (will be x-axis) and then I choose a latitude range over which the data is averaged over. If the two latitudes are the same, then the plot is of a single latitude..

When I run the following code this is the plot I get - https://imgur.com/a/1qRMKHt

This is the link of netCDF file - https://filebin.net/jsj7lzf62zfhvj3n

Where am I going wrong ?

`# importing Magics module

from Magics.macro import *

ref = 'xsection'

def xsection(title_size=0.4, axis_size=0.6, tick_size=0.2, thumbnail = True) :

Setting the cartesian view

projection = mmap(
    subpage_map_projection='cartesian',
    subpage_x_axis_type='geoline',
    subpage_y_axis_type='logarithmic',
    subpage_x_min_latitude=25.,
    subpage_x_max_latitude=40.,
    subpage_x_min_longitude=60.,
    subpage_x_max_longitude=100.,
    subpage_y_min=1000.,
    subpage_y_max=10.,
    )

# Vertical axis

vertical = maxis(
    axis_orientation='vertical',
    axis_grid='on',
    axis_type='logarithmic',
    axis_tick_label_height=tick_size,
    axis_tick_label_colour='charcoal',
    axis_grid_colour='charcoal',
    axis_grid_thickness=1,
    axis_grid_reference_level=0.,
    axis_grid_reference_line_style='solid',
    axis_grid_reference_thickness=1,
    axis_grid_line_style='dash',
    axis_title='on',
    axis_title_text='Pressure',
    axis_title_height=axis_size,
    )

# Horizontal axis

horizontal = maxis(
    axis_orientation='horizontal',
    axis_type='geoline',
    axis_tick_label_height=tick_size,
    axis_tick_label_colour='charcoal',
    axis_grid='on',
    axis_grid_colour='charcoal',
    axis_grid_thickness=1,
    axis_grid_line_style='dash',
    )
# Definition of the Netcdf data and interpretation
data = mnetcdf(netcdf_filename = "section.nc",
  netcdf_value_variable = "hgt",
  netcdf_field_scaling_factor = 100000.,
  netcdf_y_variable = "levels",
  netcdf_x_variable = "longitude",
  netcdf_x_auxiliary_variable = "latitude"
)
# Definition of the contour 
contour = mcont(contour= "on",
            contour_highlight= "off",
            contour_hilo= "off",
            contour_label= "off",
            legend='on',
            contour_level_list= [-200., -100., -75., -50., -30., -20., 
                -15., -13., -11., -9., -7., -5., -3., -1., 1., 3., 5., 
                7., 9., 11., 13., 15., 20., 30., 50., 75., 100., 200],
            contour_level_selection_type= "level_list",
            contour_shade= "on",
            contour_shade_colour_list= ["rgb(0,0,0.3)", "rgb(0,0,0.5)", 
                "rgb(0,0,0.7)", "rgb(0,0,0.9)", "rgb(0,0.15,1)", 
                "rgb(0,0.3,1)", "rgb(0,0.45,1)", "rgb(0,0.6,1)", 
                "rgb(0,0.75,1)", "rgb(0,0.85,1)", "rgb(0.2,0.95,1)", 
                "rgb(0.45,1,1)", "rgb(0.75,1,1)", "none", "rgb(1,1,0)", 
                "rgb(1,0.9,0)", "rgb(1,0.8,0)", "rgb(1,0.7,0)", 
                "rgb(1,0.6,0)", "rgb(1,0.5,0)", "rgb(1,0.4,0)", 
                "rgb(1,0.3,0)", "rgb(1,0.15,0)", "rgb(0.9,0,0)", 
                "rgb(0.7,0,0)", "rgb(0.5,0,0)", "rgb(0.3,0,0)"],
            contour_shade_colour_method= "list",
            contour_shade_method= "area_fill")
lines = \
    ['Example of Xsection to demonstrate the use of NetCDF data...', 
    '<magics_title/>']

title = mtext(
    text_lines= lines,
    text_html= 'true',
    text_justification= 'left',
    text_font_size= title_size,
    text_colour= 'charcoal',
    )
if ( thumbnail ) :
    legend = mlegend(legend='on', 
            legend_display_type='continuous', 
            legend_text_colour='charcoal',       
            legend_text_font_size=tick_size, 
            legend_box_mode = "positional",
            legend_box_x_position = 13.00,
            legend_box_y_position = 2.5,
            legend_box_x_length = 1.00,
            legend_box_y_length = 6.00,
            legend_box_blanking = "on",
            legend_border = "on",
            legend_border_colour='charcoal'
        )
else :
    legend = mlegend(legend='on', 
            legend_display_type='continuous', 
            legend_text_colour='charcoal',       
            legend_text_font_size=0.4, 
            legend_box_mode = "positional",
            legend_box_x_position = 27.00,
            legend_box_y_position = 3.00,
            legend_box_x_length = 2.00,
            legend_box_y_length = 13.00,
            legend_box_blanking = "on",
            legend_border = "on",
            legend_border_colour='charcoal'
        )
return [projection, vertical, horizontal, data, contour, legend, title]

if ( name == "main"):

Setting of the output file name

output = output(output_formats=['png'],
                output_name_first_page_number='off',
                output_name=ref)
# To the plot
plot(output, xsection(title_size=1.0, axis_size=0.8, tick_size=0.4, thumbnail=False))

`