Skip to content
Extraits de code Groupes Projets
plot_exposure_all_annual_ships_qatar.py 8,78 ko
Newer Older
  • Learn to ignore specific revisions
  • #!/usr/bin/env python3
    #  -*- coding: utf-8 -*-
    
    """
    Combine risk of oil spill with shipping density to build maps of exposure to shipping pollution
    
    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 geopandas as gpd
    from cartopy import crs as ccrs
    import pathlib
    import numpy as np
    import matplotlib
    from shapely.geometry import Polygon
    from mpl_toolkits.axes_grid1 import make_axes_locatable
    from matplotlib import colors
    from mpl_toolkits.axes_grid1.inset_locator import inset_axes
    
    basedir = "/export/miro/students/tanselain/OpenOil/"
    
    
    #%% 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_intersect_ships(list_mois, list_source, vmin, vmax, outfile) :
        
        # Call files
        path = basedir + f"Indicators/Proba_map/all_qatar/proba_risk_map_all_qatar"
        df1 = gpd.read_file(path)
        df1 = df1[~np.isnan(df1['proba'])]
    
        path = basedir + f'Indicators/Oil_shipping_exposure_map/Density_shipping_2019/shp_data/Persian_shipping_density_2019'
        df2 = gpd.read_file(path)
        df2 = df2.to_crs("EPSG:4326")
        
        # Intersect risk and shipping information
        df1_df2 = gpd.sjoin(df1,df2,how='left', op='intersects')
        df1_df2i = df1_df2.reset_index()
        df1_df2I = df1_df2i.groupby('index').agg('mean')
        df1_df2I = gpd.GeoDataFrame(df1_df2I, geometry=df1['geometry'], crs="EPSG:4326")
        
        # Creation of dn_total, a multiplication of ships probability and risk probability
        df1_df2I['DN_total']= ((df1_df2I['proba']) * (df1_df2I['density']))
        print(df1_df2I)
        
        # Creation of total, the sum of dn_total
        Total2 = df1_df2I['DN_total'].sum()
        total = round(Total2, 1)
        print ("Sum risk:",total)
        total = str(total)
            
        
    #%% Draw figure
        # Main axis
        fig = plt.figure(figsize = (12,12))     
        ax = fig.add_subplot(111)
        
        # Second axis to draw cbbox
        fc = colors.to_rgba('white')
        ec = colors.to_rgba('gray')
        fc = fc[:-1] + (0.7,)
        #cbbox = inset_axes(ax, '20%', '35%',loc=1, bbox_to_anchor=(0, 0, 0.98, 0.98), bbox_transform=ax.transAxes, borderpad=0)
        #cbbox.spines['bottom'].set_color(ec)
        #cbbox.spines['top'].set_color(ec) 
        #cbbox.spines['right'].set_color(ec)
        #cbbox.spines['left'].set_color(ec)
        #cbbox.tick_params(axis='both', left='off', top='off', right='off', bottom='off', labelleft='off', labeltop='off', labelright='off', labelbottom='off')
        #cbbox.set_facecolor(fc)
        
        # Subaxis to add colorbar to cbbox
        #subax = inset_axes(cbbox,
                        #width="35%", # width = 30% of parent_bbox
                        #height="82%", # height : 1 inch
                        #loc = 6,
                        #bbox_to_anchor=(0.1,0,1,0.9),
                        #bbox_transform=cbbox.transAxes,
                        #borderpad=0,
                        #)
        
        # Set and identify sensitive areas + additionnal texts
        list_source2 = []
        for j in range(len(list_source)) :
            src = list_source[j].split("_")
            for i in range(len(src)):
                if src[i] == "al": continue
                if src[i] == "port": continue
                src[i] = src[i][0].upper() + src[i][1:]
            source2 = " ".join(src)
            list_source2.append(source2)
        #mois = mois[0].upper() + mois[1:]
        
        #ax.text(50.05, 24.1, f'{mois}', fontsize=18, horizontalalignment = 'left', verticalalignment = 'center', zorder=6)
        #ax.text(50.05, 24.25, f'Total exposure: {total}', fontsize=18, horizontalalignment = 'left', verticalalignment = 'center', zorder=6)
    
        ax.scatter(51.55, 25.9464, s=50, c='#1f77b4' ,zorder=5)
        ax.text(51.355, 25.89, list_source2[0], color = 'k', fontsize=12, fontweight='bold',horizontalalignment = 'center', verticalalignment = 'center',zorder=5)
    
        ax.scatter(51.6269, 25.2061, s=50, c='#1f77b4' ,zorder=5)
        ax.text(51.385, 25.2, list_source2[1], color = 'k', fontsize=12, fontweight='bold' , horizontalalignment = 'center', verticalalignment = 'center',zorder=5)
    
        ax.scatter(51.633, 25.1135, s=50, c='#1f77b4' ,zorder=5)
        ax.text(51.355, 25.1135, list_source2[2], color = 'k', fontsize=12, fontweight='bold', horizontalalignment = 'center', verticalalignment = 'center',zorder=5)
        
        ax.scatter(51.6682, 25.9204, s=100, color ='limegreen', marker = '*' ,zorder=5)   
        
        # Put land, border and eez in background
        bgdir = basedir + "Indicators/Background_map/"
        gpd_landp  = gpd.read_file(bgdir+"PG_layer_poly_new/PG_layer_poly_new.shp")
        gpd_bordl  = gpd.read_file(bgdir+"natural_earth/ne_10m_admin_0_boundary_lines_land.shp")
        gpd_eezl  = gpd.read_file(bgdir+"qatar_eez/qatar_eez_line.shp")
        gpd_landp.plot(ax=ax, linewidth=0.5, edgecolors="k", color='lightgray',zorder=4)
        gpd_bordl.plot(ax=ax, linewidth=0.5, linestyle='-', color="k", zorder=4)
        gpd_eezl.plot(ax=ax, linewidth=2, color="darkgreen",zorder=4)
        
        # Plot dn_total
        column = df1_df2I['DN_total'].values.astype(float)
        vmax2 = max(column)
        print(vmax2)
        cmap = plt.get_cmap('YlOrRd')
        new_cmap = truncate_colormap(cmap, 0.2, 1, n=100) #set colormap
        
        bounds = np.array([0, vmax2/10000, vmax2/1000, vmax2/100, vmax2]) #0, 0.1 , 1, 4, 20
        stretched_bounds = np.interp(np.linspace(0, 1, 257), np.linspace(0, 1, len(bounds)), bounds)
        norm = matplotlib.colors.BoundaryNorm(stretched_bounds, ncolors=256)
        
        #pcm = df1_df2I.plot(column='DN_total', ax=ax, cmap= new_cmap, norm=norm, vmax = vmax2)
        
        df1_df2I1 = df1_df2I.loc[df1_df2I['DN_total'] < vmax2/100000]
        df1_df2I2 = df1_df2I.loc[(df1_df2I['DN_total'] >= vmax2/100000) & (df1_df2I['DN_total'] < vmax2/10000)]
        df1_df2I3 = df1_df2I.loc[(df1_df2I['DN_total'] >= vmax2/10000) & (df1_df2I['DN_total'] < vmax2/1000)]
        df1_df2I4 = df1_df2I.loc[(df1_df2I['DN_total'] >= vmax2/1000) & (df1_df2I['DN_total'] < vmax2/100)]
        df1_df2I5 = df1_df2I.loc[(df1_df2I['DN_total'] >= vmax2/100) & (df1_df2I['DN_total'] <= vmax2)]
        
        df1_df2I1.plot(column='DN_total',ax=ax, color=new_cmap(0), zorder=0.8)
        df1_df2I2.plot(column='DN_total',ax=ax, color=new_cmap(64), zorder=0.8)
        df1_df2I3.plot(column='DN_total',ax=ax, color=new_cmap(128), zorder=0.8)
        df1_df2I4.plot(column='DN_total',ax=ax, color=new_cmap(192), zorder=0.8)
        df1_df2I5.plot(column='DN_total',ax=ax, color=new_cmap(256), zorder=0.8)
        
        #Creation of separated legend
        Patch1 = matplotlib.patches.Rectangle((0, 0), 0, 0, color = new_cmap(0))
        Patch2 = matplotlib.patches.Rectangle((0, 0), 0, 0, color = new_cmap(64))
        Patch3 = matplotlib.patches.Rectangle((0, 0), 0, 0, color = new_cmap(128))
        Patch4 = matplotlib.patches.Rectangle((0, 0), 0, 0, color = new_cmap(192))
        Patch5 = matplotlib.patches.Rectangle((0, 0), 0, 0, color = new_cmap(256))
        legend = plt.legend([Patch5, Patch4, Patch3, Patch2, Patch1], ['1% - 100%', '0.1% - 1%', '0.01% - 0.1%', '0.001% - 0.01%','0% - 0.001%'], markerscale = 100, title = 'Exposure', loc= 'upper right', frameon = True, fontsize = 18)
        legend.get_title().set_fontsize('18')
    
        ax.set_aspect('equal')
        
        # Set axis
        ax.set_xlim(50, 53.2)
        ax.set_ylim(24, 27.5)
        yticks = np.arange(24, 28, .5)
        ax.set_yticks(yticks)
        ax.set_yticklabels([f"{y}°N" for y in yticks])
        ax.tick_params(labelleft=True, labelbottom = True) 
        xticks = np.arange(50,53.5,0.5)
        ax.set_xticks(xticks)
        ax.set_xticklabels([f"{x}°E" for x in xticks])
        
        # Set colorbar
        #cticks = bounds 
        #cbar = matplotlib.colorbar.ColorbarBase(subax, cmap=new_cmap, norm=norm , ticks=[cticks])
        #ctickslab = cticks
        #ctickslab= np.round(ctickslab,2)
        #cbar.ax.set_yticklabels(['0%', '0.01%', '0.1%', '1%','100%']) #'Zero', 'Low', 'Medium', 'High','Very high'
        #cbar.set_label(label = 'Exposure', labelpad=-45, y=1.15, rotation=0, fontsize = 14)
        ax.set_aspect("equal")
        
        
    #%% Save in png and shp
        save = outfile + f"exposure_annual_shipping_map_all_qatar2.png"
        plt.savefig(save, dpi=300, bbox_inches="tight")
        saveshp = outfile + f"exposure_annual_shipping_map_all_qatar2"
        df1_df2I.to_file(saveshp)       
        plt.close(fig)
         
         
    #%% Beginning of loops
    # Give the list of month assessed and name of sensitive areas    
    list_source = ['ras_laffan', 'abu_fontas', 'umm_al_houl','ras_laffan_port']  
    list_mois = ['january', 'february', 'march', 'april', 'may', 'june', 'july', 'august', 'september', 'october', 'november', 'december']
    #list_mois = ['march']
    
    # Start loops for main function
    fname = basedir + f"Indicators/Oil_shipping_exposure_map/all_annual_qatar/"
    plot_intersect_ships(list_mois, list_source, 0, 5, fname)