# (c) British Crown Copyright 2020, the Met Office.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#   * Redistributions of source code must retain the above copyright notice,
#     this list of conditions and the following disclaimer.
#   * Redistributions in binary form must reproduce the above copyright notice,
#     this list of conditions and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#   * Neither the name of the Met Office nor the names of its contributors may
#     be used to endorse or promote products derived from this software without
#     specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
# THE POSSIBILITY OF SUCH DAMAGE.

import netCDF4,argparse,sys
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import os

def collapse_dimensions_for_plotting(longitude, latitude, vname, vx, vd, dims):
    """
    Pre-processing of COSP variable for plotting.

    Arguments:
        longitude: longitude from COSP output file [loc].
        latitude: latitude from COSP output file [loc].
        vname: variable name.
        vx: variable from COSP output file [..., loc]
        vd: dictionary with metadata about variable.
        dims: dictionary with additional dimensions.
    Return:
        x: x axis values.
        y: y axis values.
        z: data array for plotting.
        d: dictionary with plot configuration.
    """
    yflip = False
    xticks, yticks, xticks_labels, yticks_labels, xlabel, ylabel, vmax = (None,)*7
    if vd['plot_type'] == 'map':
        d = None
        x = longitude[0]
        y = latitude[:,0]
        z = vx
        if vname == 'parasolGrid_refl': z = vx[2]
        # Roll longitude if there are values > 180
        m = (x > 180.0)
        if m is not None:
            Nroll = longitude.shape[1] // 2
            x[m] = x[m] - 360.0
            x = np.roll(x, Nroll)
            z = np.roll(z,Nroll,axis=1)
        # Calculate latitude and longitude edge points.
        # Assume they are increasing monotonically.
        # Extend length to N+2 and calculate midpoints.
        x = midpoints_to_edges(x)
        y = midpoints_to_edges(y)
        xticks = np.arange(-180,181,60)
        yticks = np.arange(-90,91,30)
        xlabel = 'Longitude (deg)'
        ylabel = 'Latitude (deg)'
    if vd['plot_type'] == '2Dhist':
        weights = np.cos(latitude * np.pi / 180.0)
        weights = weights / weights.sum()
        z = np.sum(vx * weights, axis=2)
        x = np.arange(z.shape[1]+1)
        y = np.arange(z.shape[0]+1)
    if vd['plot_type'] == 'zonal_cross_section':
        z = np.average(vx, axis=2)
        x = midpoints_to_edges(latitude[:,0])
        y = np.arange(z.shape[0] + 1)
    if vd['plot_type'] in ('2Dhist','zonal_cross_section'):
        if vd['xaxis_type'] == 'tau7':
            xticks_labels = ('0', '0.3', '1.3', '3.6', '9.4', '23', '60', '')
            xticks = x
            xlabel = 'Cloud optical depth'
        if vd['xaxis_type'] == 'cloudsat_DBZE_BINS':
            x = np.arange(-50,26,5)
            xticks = x
            xticks_labels = None
            xlabel = 'Radar reflectivity (dBZ)'
        if vd['xaxis_type'] == 'SR_BINS':
            xticks_labels = ('0', '0.01', '1.2', '3', '5', '7', '10', '15', '20',
                             '25', '30', '40', '50', '60', '80', '')
            xticks = x
            xlabel = 'Lidar scattering ratio'
        if vd['xaxis_type'] == 'latitude':
            xticks_labels = None
            xticks = np.arange(-90,91,30)
            xlabel = 'Latitude (deg)'
        if vd['yaxis_type'] == 'pres7':
            yticks_labels = ('1000', '800', '680', '560', '440', '310', '180','')
            yticks = y
            ylabel = 'Cloud Top Pressure (hPa)'
        if vd['yaxis_type'] == 'hgt16':
            yticks_labels =  ('', '0', '500', '1000', '1500', '2000', '2500',
                              '3000', '4000', '5000', '7000', '9000', '11000',
                              '13000', '15000', '17000', '')
            yticks = y
            ylabel = 'Cloud Top Height (m)'
        if vd['yaxis_type'] == 'REICE_MODIS':
            yticks_labels =  ('0', '10', '20', '30', '40', '60', '90')
            yticks = y
            ylabel = 'Ice particle size (micron)'
        if vd['yaxis_type'] == 'RELIQ_MODIS':
            yticks_labels = ('0', '8', '10', '13', '15', '20', '30')
            yticks = y
            ylabel = 'Liquid particle size (micron)'
        if vd['yaxis_type'] == 'levStat':
            y = 480*np.arange(41)
            yticks = y[0::4]
            yticks_labels = None
            ylabel = 'Altitude (m)'
            yflip = True
        if vd['yaxis_type'] == 'lev':
            yticks = y[0::4]
            yticks_labels = None
            ylabel = 'Model level'
            yflip = True
    # Extra processing for specific variables
    vmax = None
    if vname == 'cfadLidarsr355': vmax = 0.03
    if vname == 'cfadLidarsr532': vmax = 0.03
    if vname == 'cfadLidarsr532gr': vmax = 0.03
    if vname == 'cfadDbze94': vmax = 0.05
    if vname == 'iwpmodis': vmax = 2.0
    if vname == 'lwpmodis': vmax = 1.0
    if vname == 'tauisccp': vmax = 100.0
    if vname == 'tautmodis': vmax = 100.0
    if vname == 'tauwmodis': vmax = 100.0
    d = {'xticks':xticks,
         'yticks':yticks,
         'xticks_labels':xticks_labels,
         'yticks_labels':yticks_labels,
         'xlabel':xlabel,
         'ylabel':ylabel,
         'vmax':vmax}
    # Flip y axis?
    if yflip: z = np.flip(z, axis=0)
    return x, y, z, d

def midpoints_to_edges(x):
    """
    Calculate edge points. Midpoints must increase monotonically.

    Arguments:
        x: vector with mid points. Dimension N.
    Return:
        y: numpy vector with edges. Dimension N+1.
    """
    y = np.append(np.append(2 * x[0] - x[1], x), 2 * x[-1] - x[-2])
    return 0.5 * (y[1:] + y[:-1])

def plot_pcolormesh(x, y, v, d, fig_name, title=None, coastlines=False):
    """
    Plot pcolormesh and write the output to a png file.

    Arguments:
        x: x axis values.
        y: y axis values.
        v: data array. Dimensions [Nx,Ny,Np]
        d: dictionary with plot configuration.
        fig_name: output file name.
    Keywords:
        title: plot title.
        coastlines: plot coast lines.
    """
    fig = plt.figure(figsize=(10,5))
    cmap = plt.get_cmap('YlOrRd', 20)
    cmap.set_bad('grey', 1)
    if coastlines:
        ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=0))
        ax.coastlines()
    h = plt.pcolormesh(x, y, v, cmap=cmap, vmax=d['vmax'])
    if d['xticks_labels']:
        plt.xticks(d['xticks'],d['xticks_labels'])
    else:
        plt.xticks(d['xticks'])
    if d['yticks_labels']:
        plt.yticks(d['yticks'],d['yticks_labels'])
    else:
        plt.yticks(d['yticks'])
    plt.xlabel(d['xlabel'])
    plt.ylabel(d['ylabel'])
    plt.colorbar(h,orientation='vertical')
    if title is not None: plt.title(title)
    plt.savefig(fig_name, dpi=200)
    plt.close()

def read_dimensions(fname):
    """
    Read useful dimensions from COSP output file.

    Arguments:
        fname: path to NetCDF file.
    Return:
        d: dictionary with the following dimensions:
           'cloudsat_DBZE_BINS', 'hgt16', 'REICE_MODIS',
           'RELIQ_MODIS', 'levStat', 'SR_BINS', 'lev'
    """
    dim_names = ['cloudsat_DBZE_BINS', 'hgt16', 'REICE_MODIS', 'RELIQ_MODIS',
                 'levStat', 'SR_BINS', 'lev']
    d = {}
    f_id = netCDF4.Dataset(fname, 'r')
    for dim in dim_names:
        d[dim] = f_id.variables[dim][:]
    f_id.close()
    return d

def read_var_to_masked_array(fname, vname, fill_value, Nlat_lon = None):
    """
    Reads a variable from a NetCDF file, and produces a masked array.

    Arguments:
        fname: path to NetCDF file.
        vname: variable name.
        fill_value: missing data value.
    Keywords:
        Nlat_lon: tuple (Nrows, Ncols). If defined, variable is
                  reshaped to a lat-lon grid.
    Return:
        x: variable data array.
        lon: longitude array.
        lat: latitude array
        units: units attribute.
        long_name: long name attribute.
    """
    f_id = netCDF4.Dataset(fname, 'r')
    x = np.ma.masked_equal(f_id.variables[vname][:], fill_value)
    lon = np.ma.masked_equal(f_id.variables['longitude'][:], fill_value)
    lat = np.ma.masked_equal(f_id.variables['latitude'][:], fill_value)
    units = f_id.variables[vname].getncattr('units')
    long_name = f_id.variables[vname].getncattr('long_name')
    f_id.close()

    if Nlat_lon is not None:
        x = np.reshape(x, x.shape[:-1]+Nlat_lon)
        lon = np.reshape(lon, lon.shape[:-1] + Nlat_lon)
        lat = np.reshape(lat, lat.shape[:-1] + Nlat_lon)
    return x, lon, lat, units, long_name

def produce_cosp_summary_plots(fname, variables, output_dir, Nlat_lon = None):
    """
    Wrapper function that iterates over a list of COSP variables and produces
    a PNG figure for each of them.

    Arguments:
        fname: COSP output filename.
        variables: list of variable names.
        output_dir: output directory.
    Keywords:
        Nlat_lon: tuple with (lat,lon) dimensions model's grid.
    """
    fill_value = -1.0e30
    dimensions = read_dimensions(fname)
    for vname, vd in variables.items():
        new_shape = None
        if vd['reshape']: new_shape = Nlat_lon
        vx, longitude, latitude, units, long_name = read_var_to_masked_array(fname, vname, fill_value, Nlat_lon = new_shape)
        x, y, z, pkw = collapse_dimensions_for_plotting(longitude, latitude, vname, vx, vd, dimensions)
        title = long_name + ' (' + units + ')'
        fig_name = os.path.join(output_dir, ".".join([os.path.basename(fname), vname, 'png']))
        coastlines = False
        if vd['plot_type'] == 'map': coastlines = True
        plot_pcolormesh(x, y, z, pkw, fig_name, title=title, coastlines=coastlines)

def variable2D_metadata(var_list, fname):
    """
    Return dictionary with metadata for each variable.

    Arguments:
        var_list: list of variable names.
        fname: COSP output filename.
    Return:
        d: dictionary of dictionaries with relevant metadata.
    """
    map_dims = (('loc',),('PARASOL_NREFL','loc'))
    hist2D_dims = (('pres7', 'tau7', 'loc'),
                  ('levStat', 'cloudsat_DBZE_BINS', 'loc'),
                  ('hgt16', 'tau7', 'loc'),
                  ('REICE_MODIS', 'tau7', 'loc'),
                  ('RELIQ_MODIS', 'tau7', 'loc'),
                  ('levStat', 'SR_BINS', 'loc'))
    zcs_dims = (('levStat','loc'), ('lev','loc'))
    f_id = netCDF4.Dataset(fname, 'r')
    vmeta = {}
    for vname in var_list:
        x = f_id.variables[vname]
        # Standard map
        if x.dimensions in map_dims:
            vmeta[vname] = {'plot_type':'map', 'reshape':True}
        # 2D histograms
        if x.dimensions in hist2D_dims:
            vmeta[vname] = {'plot_type':'2Dhist', 'reshape':False,
                            'xaxis_type': x.dimensions[1],
                            'yaxis_type': x.dimensions[0]}
        # Zonal cross section
        if x.dimensions in zcs_dims:
            vmeta[vname] = {'plot_type':'zonal_cross_section', 'reshape':True,
                            'xaxis_type': 'latitude',
                            'yaxis_type': x.dimensions[0]}
    f_id.close()
    return vmeta



#######################
# Main
#######################
if __name__ == '__main__':

    # Command line arguments.
    parser = argparse.ArgumentParser()
    parser.add_argument("--tst_file", default="./data/outputs/UKMO/cosp2_output.um_global.nc",
                        help="Test output file.")
    parser.add_argument("--out_dir", default="./data/outputs/UKMO",
                        help="Output directory.")
    parser.add_argument("--Nlat", type=int, default=36, help="Number of latitude points.")
    parser.add_argument("--Nlon", type=int, default=48, help="Number of longitude points.")
    args = parser.parse_args()

    # Dictonaries with list of variables.
    v2D_maps_names = ['cllcalipsoice', 'clmcalipsoice', 'clhcalipsoice', 'cltcalipsoice', 'cllcalipsoliq', 'clmcalipsoliq',
                      'clhcalipsoliq', 'cltcalipsoliq', 'cllcalipsoun', 'clmcalipsoun', 'clhcalipsoun', 'cltcalipsoun',
                      'cllcalipso', 'clmcalipso', 'clhcalipso', 'cltcalipso', 'clopaquecalipso', 'clthincalipso',
                      'clzopaquecalipso', 'clopaquetemp', 'clthintemp', 'clzopaquetemp', 'clopaquemeanz', 'clthinmeanz',
                      'clthinemis', 'clopaquemeanzse', 'clthinmeanzse', 'clzopaquecalipsose', 'cllgrLidar532',
                      'clmgrLidar532', 'clhgrLidar532', 'cltgrLidar532', 'cllatlid', 'clmatlid', 'clhatlid',
                      'cltatlid', 'cloudsatpia', 'cltisccp', 'meantbisccp', 'meantbclrisccp', 'pctisccp', 'tauisccp',
                      'albisccp', 'misr_meanztop', 'misr_cldarea', 'cltmodis', 'clwmodis', 'climodis', 'clhmodis',
                      'clmmodis', 'cllmodis', 'tautmodis', 'tauwmodis', 'tauimodis', 'tautlogmodis', 'tauwlogmodis',
                      'tauilogmodis', 'reffclwmodis', 'reffclimodis', 'pctmodis', 'lwpmodis', 'iwpmodis',
                      'cltlidarradar', 'cloudsat_tcc', 'cloudsat_tcc2','parasolGrid_refl']
    v2D_hists_names = ['clisccp', 'clmodis', 'cfadDbze94', 'clMISR',
                       'modis_Optical_Thickness_vs_ReffICE',
                       'modis_Optical_Thickness_vs_ReffLIQ',
                       'cfadLidarsr532', 'cfadLidarsr532gr', 'cfadLidarsr355']
    v2D_zcs_names = ['clcalipsoice','clcalipsoliq','clcalipsoun','clcalipsotmp','clcalipsotmpice','clcalipsotmpliq',
                     'clcalipsotmpun','clcalipso','clcalipsoopaque','clcalipsothin','clcalipsozopaque',
                     'clcalipsoopacity','clgrLidar532','clatlid','clcalipso2',
                     'lidarBetaMol532gr','lidarBetaMol532','lidarBetaMol355']
    v2D_all_names = v2D_maps_names + v2D_hists_names + v2D_zcs_names
    # Plots for these variables are not yet developed
    # atb532_perp(lev, cosp_scol, loc);
    # atb532(lev, cosp_scol, loc);
    # calipso_tau(lev, cosp_scol, loc);
    # atb532gr(lev, cosp_scol, loc);
    # atb355(lev, cosp_scol, loc);
    # dbze94(lev, cosp_scol, loc);
    # parasolPix_refl(PARASOL_NREFL, cosp_scol, loc);
    # boxtauisccp(cosp_scol, loc);
    # boxptopisccp(cosp_scol, loc);

    # Build dictionary with metadata and produce plots.
    vars2D  = variable2D_metadata(v2D_all_names, args.tst_file)
    produce_cosp_summary_plots(args.tst_file, vars2D, args.out_dir,
                               Nlat_lon=(args.Nlat, args.Nlon))

    print("===== Summary plots produced =====")

# (c) British Crown Copyright 2020, the Met Office.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#   * Redistributions of source code must retain the above copyright notice,
#     this list of conditions and the following disclaimer.
#   * Redistributions in binary form must reproduce the above copyright notice,
#     this list of conditions and the following disclaimer in the documentation
#     and/or other materials provided with the distribution.
#   * Neither the name of the Met Office nor the names of its contributors may
#     be used to endorse or promote products derived from this software without
#     specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
# THE POSSIBILITY OF SUCH DAMAGE.

import netCDF4,argparse,sys
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import os

def collapse_dimensions_for_plotting(longitude, latitude, vname, vx, vd, dims):
    """
    Pre-processing of COSP variable for plotting.

    Arguments:
        longitude: longitude from COSP output file [loc].
        latitude: latitude from COSP output file [loc].
        vname: variable name.
        vx: variable from COSP output file [..., loc]
        vd: dictionary with metadata about variable.
        dims: dictionary with additional dimensions.
    Return:
        x: x axis values.
        y: y axis values.
        z: data array for plotting.
        d: dictionary with plot configuration.
    """
    yflip = False
    xticks, yticks, xticks_labels, yticks_labels, xlabel, ylabel, vmax = (None,)*7
    if vd['plot_type'] == 'map':
        d = None
        x = longitude[0]
        y = latitude[:,0]
        z = vx
        if vname == 'parasolGrid_refl': z = vx[2]
        # Roll longitude if there are values > 180
        m = (x > 180.0)
        if m is not None:
            Nroll = longitude.shape[1] // 2
            x[m] = x[m] - 360.0
            x = np.roll(x, Nroll)
            z = np.roll(z,Nroll,axis=1)
        # Calculate latitude and longitude edge points.
        # Assume they are increasing monotonically.
        # Extend length to N+2 and calculate midpoints.
        x = midpoints_to_edges(x)
        y = midpoints_to_edges(y)
        xticks = np.arange(-180,181,60)
        yticks = np.arange(-90,91,30)
        xlabel = 'Longitude (deg)'
        ylabel = 'Latitude (deg)'
    if vd['plot_type'] == '2Dhist':
        weights = np.cos(latitude * np.pi / 180.0)
        weights = weights / weights.sum()
        z = np.sum(vx * weights, axis=2)
        x = np.arange(z.shape[1]+1)
        y = np.arange(z.shape[0]+1)
    if vd['plot_type'] == 'zonal_cross_section':
        z = np.average(vx, axis=2)
        x = midpoints_to_edges(latitude[:,0])
        y = np.arange(z.shape[0] + 1)
    if vd['plot_type'] in ('2Dhist','zonal_cross_section'):
        if vd['xaxis_type'] == 'tau7':
            xticks_labels = ('0', '0.3', '1.3', '3.6', '9.4', '23', '60', '')
            xticks = x
            xlabel = 'Cloud optical depth'
        if vd['xaxis_type'] == 'cloudsat_DBZE_BINS':
            x = np.arange(-50,26,5)
            xticks = x
            xticks_labels = None
            xlabel = 'Radar reflectivity (dBZ)'
        if vd['xaxis_type'] == 'SR_BINS':
            xticks_labels = ('0', '0.01', '1.2', '3', '5', '7', '10', '15', '20',
                             '25', '30', '40', '50', '60', '80', '')
            xticks = x
            xlabel = 'Lidar scattering ratio'
        if vd['xaxis_type'] == 'latitude':
            xticks_labels = None
            xticks = np.arange(-90,91,30)
            xlabel = 'Latitude (deg)'
        if vd['yaxis_type'] == 'pres7':
            yticks_labels = ('1000', '800', '680', '560', '440', '310', '180','')
            yticks = y
            ylabel = 'Cloud Top Pressure (hPa)'
        if vd['yaxis_type'] == 'hgt16':
            yticks_labels =  ('', '0', '500', '1000', '1500', '2000', '2500',
                              '3000', '4000', '5000', '7000', '9000', '11000',
                              '13000', '15000', '17000', '')
            yticks = y
            ylabel = 'Cloud Top Height (m)'
        if vd['yaxis_type'] == 'REICE_MODIS':
            yticks_labels =  ('0', '10', '20', '30', '40', '60', '90')
            yticks = y
            ylabel = 'Ice particle size (micron)'
        if vd['yaxis_type'] == 'RELIQ_MODIS':
            yticks_labels = ('0', '8', '10', '13', '15', '20', '30')
            yticks = y
            ylabel = 'Liquid particle size (micron)'
        if vd['yaxis_type'] == 'levStat':
            y = 480*np.arange(41)
            yticks = y[0::4]
            yticks_labels = None
            ylabel = 'Altitude (m)'
            yflip = True
        if vd['yaxis_type'] == 'lev':
            yticks = y[0::4]
            yticks_labels = None
            ylabel = 'Model level'
            yflip = True
    # Extra processing for specific variables
    vmax = None
    if vname == 'cfadLidarsr355': vmax = 0.03
    if vname == 'cfadLidarsr532': vmax = 0.03
    if vname == 'cfadLidarsr532gr': vmax = 0.03
    if vname == 'cfadDbze94': vmax = 0.05
    if vname == 'iwpmodis': vmax = 2.0
    if vname == 'lwpmodis': vmax = 1.0
    if vname == 'tauisccp': vmax = 100.0
    if vname == 'tautmodis': vmax = 100.0
    if vname == 'tauwmodis': vmax = 100.0
    d = {'xticks':xticks,
         'yticks':yticks,
         'xticks_labels':xticks_labels,
         'yticks_labels':yticks_labels,
         'xlabel':xlabel,
         'ylabel':ylabel,
         'vmax':vmax}
    # Flip y axis?
    if yflip: z = np.flip(z, axis=0)
    return x, y, z, d

def midpoints_to_edges(x):
    """
    Calculate edge points. Midpoints must increase monotonically.

    Arguments:
        x: vector with mid points. Dimension N.
    Return:
        y: numpy vector with edges. Dimension N+1.
    """
    y = np.append(np.append(2 * x[0] - x[1], x), 2 * x[-1] - x[-2])
    return 0.5 * (y[1:] + y[:-1])

def plot_pcolormesh(x, y, v, d, fig_name, title=None, coastlines=False):
    """
    Plot pcolormesh and write the output to a png file.

    Arguments:
        x: x axis values.
        y: y axis values.
        v: data array. Dimensions [Nx,Ny,Np]
        d: dictionary with plot configuration.
        fig_name: output file name.
    Keywords:
        title: plot title.
        coastlines: plot coast lines.
    """
    fig = plt.figure(figsize=(10,5))
    cmap = plt.get_cmap('YlOrRd', 20)
    cmap.set_bad('grey', 1)
    if coastlines:
        ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=0))
        ax.coastlines()
    h = plt.pcolormesh(x, y, v, cmap=cmap, vmax=d['vmax'])
    if d['xticks_labels']:
        plt.xticks(d['xticks'],d['xticks_labels'])
    else:
        plt.xticks(d['xticks'])
    if d['yticks_labels']:
        plt.yticks(d['yticks'],d['yticks_labels'])
    else:
        plt.yticks(d['yticks'])
    plt.xlabel(d['xlabel'])
    plt.ylabel(d['ylabel'])
    plt.colorbar(h,orientation='vertical')
    if title is not None: plt.title(title)
    plt.savefig(fig_name, dpi=200)
    plt.close()

def read_dimensions(fname):
    """
    Read useful dimensions from COSP output file.

    Arguments:
        fname: path to NetCDF file.
    Return:
        d: dictionary with the following dimensions:
           'cloudsat_DBZE_BINS', 'hgt16', 'REICE_MODIS',
           'RELIQ_MODIS', 'levStat', 'SR_BINS', 'lev'
    """
    dim_names = ['cloudsat_DBZE_BINS', 'hgt16', 'REICE_MODIS', 'RELIQ_MODIS',
                 'levStat', 'SR_BINS', 'lev']
    d = {}
    f_id = netCDF4.Dataset(fname, 'r')
    for dim in dim_names:
        d[dim] = f_id.variables[dim][:]
    f_id.close()
    return d

def read_var_to_masked_array(fname, vname, fill_value, Nlat_lon = None):
    """
    Reads a variable from a NetCDF file, and produces a masked array.

    Arguments:
        fname: path to NetCDF file.
        vname: variable name.
        fill_value: missing data value.
    Keywords:
        Nlat_lon: tuple (Nrows, Ncols). If defined, variable is
                  reshaped to a lat-lon grid.
    Return:
        x: variable data array.
        lon: longitude array.
        lat: latitude array
        units: units attribute.
        long_name: long name attribute.
    """
    f_id = netCDF4.Dataset(fname, 'r')
    x = np.ma.masked_equal(f_id.variables[vname][:], fill_value)
    lon = np.ma.masked_equal(f_id.variables['longitude'][:], fill_value)
    lat = np.ma.masked_equal(f_id.variables['latitude'][:], fill_value)
    units = f_id.variables[vname].getncattr('units')
    long_name = f_id.variables[vname].getncattr('long_name')
    f_id.close()

    if Nlat_lon is not None:
        x = np.reshape(x, x.shape[:-1]+Nlat_lon)
        lon = np.reshape(lon, lon.shape[:-1] + Nlat_lon)
        lat = np.reshape(lat, lat.shape[:-1] + Nlat_lon)
    return x, lon, lat, units, long_name

def produce_cosp_summary_plots(fname, variables, output_dir, Nlat_lon = None):
    """
    Wrapper function that iterates over a list of COSP variables and produces
    a PNG figure for each of them.

    Arguments:
        fname: COSP output filename.
        variables: list of variable names.
        output_dir: output directory.
    Keywords:
        Nlat_lon: tuple with (lat,lon) dimensions model's grid.
    """
    fill_value = -1.0e30
    dimensions = read_dimensions(fname)
    for vname, vd in variables.items():
        new_shape = None
        if vd['reshape']: new_shape = Nlat_lon
        vx, longitude, latitude, units, long_name = read_var_to_masked_array(fname, vname, fill_value, Nlat_lon = new_shape)
        x, y, z, pkw = collapse_dimensions_for_plotting(longitude, latitude, vname, vx, vd, dimensions)
        title = long_name + ' (' + units + ')'
        fig_name = os.path.join(output_dir, ".".join([os.path.basename(fname), vname, 'png']))
        coastlines = False
        if vd['plot_type'] == 'map': coastlines = True
        plot_pcolormesh(x, y, z, pkw, fig_name, title=title, coastlines=coastlines)

def variable2D_metadata(var_list, fname):
    """
    Return dictionary with metadata for each variable.

    Arguments:
        var_list: list of variable names.
        fname: COSP output filename.
    Return:
        d: dictionary of dictionaries with relevant metadata.
    """
    map_dims = (('loc',),('PARASOL_NREFL','loc'))
    hist2D_dims = (('pres7', 'tau7', 'loc'),
                  ('levStat', 'cloudsat_DBZE_BINS', 'loc'),
                  ('hgt16', 'tau7', 'loc'),
                  ('REICE_MODIS', 'tau7', 'loc'),
                  ('RELIQ_MODIS', 'tau7', 'loc'),
                  ('levStat', 'SR_BINS', 'loc'))
    zcs_dims = (('levStat','loc'), ('lev','loc'))
    f_id = netCDF4.Dataset(fname, 'r')
    vmeta = {}
    for vname in var_list:
        x = f_id.variables[vname]
        # Standard map
        if x.dimensions in map_dims:
            vmeta[vname] = {'plot_type':'map', 'reshape':True}
        # 2D histograms
        if x.dimensions in hist2D_dims:
            vmeta[vname] = {'plot_type':'2Dhist', 'reshape':False,
                            'xaxis_type': x.dimensions[1],
                            'yaxis_type': x.dimensions[0]}
        # Zonal cross section
        if x.dimensions in zcs_dims:
            vmeta[vname] = {'plot_type':'zonal_cross_section', 'reshape':True,
                            'xaxis_type': 'latitude',
                            'yaxis_type': x.dimensions[0]}
    f_id.close()
    return vmeta



#######################
# Main
#######################
if __name__ == '__main__':

    # Command line arguments.
    parser = argparse.ArgumentParser()
    parser.add_argument("--tst_file", default="./data/outputs/UKMO/cosp2_output.um_global.nc",
                        help="Test output file.")
    parser.add_argument("--out_dir", default="./data/outputs/UKMO",
                        help="Output directory.")
    parser.add_argument("--Nlat", type=int, default=36, help="Number of latitude points.")
    parser.add_argument("--Nlon", type=int, default=48, help="Number of longitude points.")
    args = parser.parse_args()

    # Dictonaries with list of variables.
    v2D_maps_names = ['cllcalipsoice', 'clmcalipsoice', 'clhcalipsoice', 'cltcalipsoice', 'cllcalipsoliq', 'clmcalipsoliq',
                      'clhcalipsoliq', 'cltcalipsoliq', 'cllcalipsoun', 'clmcalipsoun', 'clhcalipsoun', 'cltcalipsoun',
                      'cllcalipso', 'clmcalipso', 'clhcalipso', 'cltcalipso', 'clopaquecalipso', 'clthincalipso',
                      'clzopaquecalipso', 'clopaquetemp', 'clthintemp', 'clzopaquetemp', 'clopaquemeanz', 'clthinmeanz',
                      'clthinemis', 'clopaquemeanzse', 'clthinmeanzse', 'clzopaquecalipsose', 'cllgrLidar532',
                      'clmgrLidar532', 'clhgrLidar532', 'cltgrLidar532', 'cllatlid', 'clmatlid', 'clhatlid',
                      'cltatlid', 'cloudsatpia', 'cltisccp', 'meantbisccp', 'meantbclrisccp', 'pctisccp', 'tauisccp',
                      'albisccp', 'misr_meanztop', 'misr_cldarea', 'cltmodis', 'clwmodis', 'climodis', 'clhmodis',
                      'clmmodis', 'cllmodis', 'tautmodis', 'tauwmodis', 'tauimodis', 'tautlogmodis', 'tauwlogmodis',
                      'tauilogmodis', 'reffclwmodis', 'reffclimodis', 'pctmodis', 'lwpmodis', 'iwpmodis',
                      'cltlidarradar', 'cloudsat_tcc', 'cloudsat_tcc2','parasolGrid_refl']
    v2D_hists_names = ['clisccp', 'clmodis', 'cfadDbze94', 'clMISR',
                       'modis_Optical_Thickness_vs_ReffICE',
                       'modis_Optical_Thickness_vs_ReffLIQ',
                       'cfadLidarsr532', 'cfadLidarsr532gr', 'cfadLidarsr355']
    v2D_zcs_names = ['clcalipsoice','clcalipsoliq','clcalipsoun','clcalipsotmp','clcalipsotmpice','clcalipsotmpliq',
                     'clcalipsotmpun','clcalipso','clcalipsoopaque','clcalipsothin','clcalipsozopaque',
                     'clcalipsoopacity','clgrLidar532','clatlid','clcalipso2',
                     'lidarBetaMol532gr','lidarBetaMol532','lidarBetaMol355']
    v2D_all_names = v2D_maps_names + v2D_hists_names + v2D_zcs_names
    # Plots for these variables are not yet developed
    # atb532_perp(lev, cosp_scol, loc);
    # atb532(lev, cosp_scol, loc);
    # calipso_tau(lev, cosp_scol, loc);
    # atb532gr(lev, cosp_scol, loc);
    # atb355(lev, cosp_scol, loc);
    # dbze94(lev, cosp_scol, loc);
    # parasolPix_refl(PARASOL_NREFL, cosp_scol, loc);
    # boxtauisccp(cosp_scol, loc);
    # boxptopisccp(cosp_scol, loc);

    # Build dictionary with metadata and produce plots.
    vars2D  = variable2D_metadata(v2D_all_names, args.tst_file)
    produce_cosp_summary_plots(args.tst_file, vars2D, args.out_dir,
                               Nlat_lon=(args.Nlat, args.Nlon))

    print("===== Summary plots produced =====")

# (c) British Crown Copyright 2020, the Met Office. # All rights reserved. # # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: # #   * Redistributions of source code must retain the above copyright notice, #     this list of conditions and the following disclaimer. #   * Redistributions in binary form must reproduce the above copyright notice, #     this list of conditions and the following disclaimer in the documentation #     and/or other materials provided with the distribution. #   * Neither the name of the Met Office nor the names of its contributors may #     be used to endorse or promote products derived from this software without #     specific prior written permission. # # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" # AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE # IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE # ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE # LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF # SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS # INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN # CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) # ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF # THE POSSIBILITY OF SUCH DAMAGE. Pre-processing of COSP variable for plotting.      Arguments:         longitude: longitude from COSP output file [loc].         latitude: latitude from COSP output file [loc].         vname: variable name.         vx: variable from COSP output file [..., loc]         vd: dictionary with metadata about variable.         dims: dictionary with additional dimensions.     Return:         x: x axis values.         y: y axis values.         z: data array for plotting.         d: dictionary with plot configuration. # Roll longitude if there are values > 180 # Calculate latitude and longitude edge points. # Assume they are increasing monotonically. # Extend length to N+2 and calculate midpoints. # Extra processing for specific variables # Flip y axis? Calculate edge points. Midpoints must increase monotonically.      Arguments:         x: vector with mid points. Dimension N.     Return:         y: numpy vector with edges. Dimension N+1. Plot pcolormesh and write the output to a png file.      Arguments:         x: x axis values.         y: y axis values.         v: data array. Dimensions [Nx,Ny,Np]         d: dictionary with plot configuration.         fig_name: output file name.     Keywords:         title: plot title.         coastlines: plot coast lines. Read useful dimensions from COSP output file.      Arguments:         fname: path to NetCDF file.     Return:         d: dictionary with the following dimensions:            'cloudsat_DBZE_BINS', 'hgt16', 'REICE_MODIS',            'RELIQ_MODIS', 'levStat', 'SR_BINS', 'lev' Reads a variable from a NetCDF file, and produces a masked array.      Arguments:         fname: path to NetCDF file.         vname: variable name.         fill_value: missing data value.     Keywords:         Nlat_lon: tuple (Nrows, Ncols). If defined, variable is                   reshaped to a lat-lon grid.     Return:         x: variable data array.         lon: longitude array.         lat: latitude array         units: units attribute.         long_name: long name attribute. Wrapper function that iterates over a list of COSP variables and produces     a PNG figure for each of them.      Arguments:         fname: COSP output filename.         variables: list of variable names.         output_dir: output directory.     Keywords:         Nlat_lon: tuple with (lat,lon) dimensions model's grid. Return dictionary with metadata for each variable.      Arguments:         var_list: list of variable names.         fname: COSP output filename.     Return:         d: dictionary of dictionaries with relevant metadata. # Standard map # 2D histograms # Zonal cross section ####################### # Main ####################### # Command line arguments. # Dictonaries with list of variables. # Plots for these variables are not yet developed # atb532_perp(lev, cosp_scol, loc); # atb532(lev, cosp_scol, loc); # calipso_tau(lev, cosp_scol, loc); # atb532gr(lev, cosp_scol, loc); # atb355(lev, cosp_scol, loc); # dbze94(lev, cosp_scol, loc); # parasolPix_refl(PARASOL_NREFL, cosp_scol, loc); # boxtauisccp(cosp_scol, loc); # boxptopisccp(cosp_scol, loc); # Build dictionary with metadata and produce plots.