Skip to content
Extraits de code Groupes Projets
temporal_series_risk_averaged.py 5,69 ko
Newer Older
  • Learn to ignore specific revisions
  • #!/usr/bin/env python3
    #  -*- coding: utf-8 -*-
    
    """
    Generate a graph of the monthly-averaged total risk of oil spill at sensitive sites in Qatar
    
    Author: Thomas Anselain, Earth and Life Institute, UCLouvain, Belgium
    Last modified: 7 April 2022
    """
    #%% Import packages and set up
    import matplotlib.pyplot as plt
    import numpy as np
    import pandas as pd
    from matplotlib import pyplot
    from pylab import *
    import statistics
    from matplotlib import colors
    
    basedir = "/export/miro/students/tanselain/OpenOil/"
    list_mois2 = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
    
    
    #%% Define colormap
    def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100):
        from matplotlib import colors 
        new_cmap = colors.LinearSegmentedColormap.from_list(
            'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval),
            cmap(np.linspace(minval, maxval, n)))
        return new_cmap
    
    
    #%% Define main function
    def plot_series(list_year, list_mois2, vmin, vmax, outfile) :
        
        # Call pkl files and preprocess
        d = {}
        for year in list_year:
            connect_df = pd.read_pickle(basedir + f"Indicators/Proba_map/{year}/proba_map_{year}.pkl")
            connect_df = connect_df.replace(np.nan, 0)
            if year == '2016' :
                print(connect_df)
            stacked = connect_df.stack().reset_index().rename(columns={'level_0' : 'source', 'level_1' : 'month', 0 : 'total_risk'})
            stacked['total_risk'] = stacked['total_risk'].replace(0,np.nan)
            d[year] = stacked
        print(d)
        
        data = [d['2016']['source'], d['2016']['month']]
        headers = ["source", "month", "mean_total"]
        df_av = pd.concat(data, axis=1, keys=headers)
        df_av['mean_total'] = 0.0
        df_av['stdev'] = 0.0
        print(df_av)
        
        for i in range (len(df_av['mean_total'])):
            v1 = d['2016']['total_risk'][i]
            v2 = d['2017']['total_risk'][i]
            v3 = d['2018']['total_risk'][i]
            v4 = d['2019']['total_risk'][i]
            v5 = d['2020']['total_risk'][i]
            list_va = [v1, v2, v3, v4, v5]
            n_list_va = [x for x in list_va if np.isnan(x) == False]
            df_av['mean_total'][i] = statistics.mean(n_list_va)
            df_av['stdev'][i] = statistics.pstdev(n_list_va)
        print(df_av)
        
        # Take values for each sensitive areas
        indexNames1 = df_av[ df_av['source'] != 'umm_al_houl'] .index
        indexNames2 = df_av[ df_av['source'] != 'ras_laffan'] .index
        indexNames3 = df_av[ df_av['source'] != 'abu_fontas'] .index
        indexNames4 = df_av[ df_av['source'] != 'ras_laffan_port'] .index
        sta_umm = df_av.drop(indexNames1 , inplace=False)
        sta_ras = df_av.drop(indexNames2 , inplace=False)
        sta_abu = df_av.drop(indexNames3 , inplace=False)
        sta_rasp = df_av.drop(indexNames4 , inplace=False)
        
        print(sta_umm)
        
        umm_al_houl_d = sta_umm['stdev'].values.tolist()
        ras_laffan_d = sta_ras['stdev'].values.tolist()
        abu_fontas_d = sta_abu['stdev'].values.tolist()
        ras_laffan_port_d = sta_rasp['stdev'].values.tolist()
        
        umm_al_houl = sta_umm['mean_total'].values.tolist()
        ras_laffan = sta_ras['mean_total'].values.tolist()
        abu_fontas = sta_abu['mean_total'].values.tolist()
        ras_laffan_port = sta_rasp['mean_total'].values.tolist()
        
        maxi = max(ras_laffan_port)
        print(maxi)
        
        umm_al_houl[:] = [x / maxi for x in umm_al_houl]
        ras_laffan[:] = [x / maxi for x in ras_laffan]
        abu_fontas[:] = [x / maxi for x in abu_fontas]
        ras_laffan_port[:] = [x / maxi for x in ras_laffan_port]
        
        umm_al_houl_d[:] = [x / maxi for x in umm_al_houl_d]
        ras_laffan_d[:] = [x / maxi for x in ras_laffan_d]
        abu_fontas_d[:] = [x / maxi for x in abu_fontas_d]
        ras_laffan_port_d[:] = [x / maxi for x in ras_laffan_port_d]
       
            
    #%% Plot
        fig, ax = plt.subplots()
        fig.set_size_inches(9, 5)
        
    
        cmap = cm.get_cmap('Reds_r')    # PiYG
        new_cmap = truncate_colormap(cmap, 0.1, 0.9, n=5) #set colormap
        n_cmap = [new_cmap(0), new_cmap(64), new_cmap(128), new_cmap(192)]
        
        # Draw 4 lines for the 4 sensitive areas
        ax.plot(list_mois2, umm_al_houl, color=n_cmap[3], label='Umm al Houl', linewidth = 3)
        ax.plot(list_mois2, ras_laffan, color=n_cmap[1], label='Ras Laffan', linewidth = 3)
        ax.plot(list_mois2, abu_fontas, color=n_cmap[2], label='Abu Fontas', linewidth = 3)
        ax.plot(list_mois2, ras_laffan_port, color=n_cmap[0], label='Ras Laffan port', linewidth = 3)
        
        
        c = colors.to_rgba('k')
        c = c[:-1] + (0.4,)
        
        for i in range (0,3):
            n_cmap[i] = n_cmap[i][:-1] + (0.4,)
            
        #ax.errorbar(list_mois2, umm_al_houl,umm_al_houl_d, linestyle = 'None', color = c, capsize = 4)
        #ax.errorbar(list_mois2, ras_laffan,ras_laffan_d, linestyle = 'None', color = c, capsize = 4)
        #ax.errorbar(list_mois2, abu_fontas,abu_fontas_d, linestyle = 'None', color = c, capsize = 4)
        #ax.errorbar(list_mois2, ras_laffan_port,ras_laffan_port_d, linestyle = 'None', color = c, capsize = 4)
        
        
        # Set axis, legend and title
        ax.spines['top'].set_visible(False) 
        ax.spines['right'].set_visible(False)
        ax.set_ylim([vmin, vmax+0.03]) #+0.13 if errorbar  or   0.03 without colorbar
        ax.set_xlim([-0.1,11.1])
        ax.set_ylabel('Total risk') 
        ax.text(0.5,0.75,'Ras Laffan Port',color=new_cmap(0))
        ax.text(0.5,0.44,'Ras Laffan Desalination',color=new_cmap(64))
        ax.text(0.5,0.29,'Abu Fontas',color=new_cmap(128))
        ax.text(0.5,0.12,'Umm al Houl',color=new_cmap(192))
        
        # Final save
        outfile_png = outfile + '/risk_evolution_averaged.png'
        plt.savefig(outfile_png, bbox_inches='tight', dpi = 300)
    
    
    #%% Run main function
    fname = basedir + 'Indicators/Proba_map/average'
    list_year= ['2016', '2017', '2018', '2019', '2020']
    plot_series(list_year, list_mois2, 0, 1, fname)