# -*- coding: utf-8 -*- """ Created on Wed Mar 13 11:24:21 2024 @author: cenv1021 This script is used to analyze future changes in ML-based flood projections """ import pandas as pd import os,cmaps,sys,glob import numpy as np import matplotlib as mpl import matplotlib.ticker as ticker import matplotlib.pyplot as plt from matplotlib.collections import LineCollection import multiprocessing as mp import cartopy.feature as cf from multiprocessing import Process from pyogrio import read_dataframe, write_dataframe import cartopy.crs as ccrs from src.plot_utils import plot_grit import seaborn as sns import scienceplots import warnings warnings.filterwarnings("ignore") plt.style.use(['science','no-latex','nature']) # require install SciencePlots #%% --------------------------------------------- Calculating climate_predictors = ['AnnSum', 'Pvol', 'Rx1day', 'Rx5day', 'hurs', 'rsds', 'sfcWind', 'tasmax', 'tasmin', 'snw', 'AnnSum_std', 'Pvol_std', 'Rx1day_std', 'Rx5day_std'] reachNames = glob.glob('../data/GRIT/segments/GRITv05_segments_*_EPSG8857_correct_grwl0p5.gpkg') predHistCSVs = glob.glob('../results/*/GRIT_historical_flood_1981-2014_addLonLat.csv') basin_list = [] for i in reachNames: a = read_dataframe(i, read_geometry = False) basin_list += a.global_id.values.tolist() # calculate ensemble difference write = False if write: df_all = [] for predHistCSV in predHistCSVs: model = os.path.basename(os.path.dirname(predHistCSV)) # read projections try: climateHist = glob.glob('../results/%s/GRIT_historical_flood_1981-2014_addLonLat.csv'%model)[0] climate245 = glob.glob('../results/%s/GRIT_ssp245_flood_2071-2100_addLonLat.csv'%model)[0] climate585 = glob.glob('../results/%s/GRIT_ssp585_flood_2071-2100_addLonLat.csv'%model)[0] df_hist = pd.read_csv(climateHist) df_245 = pd.read_csv(climate245) df_585 = pd.read_csv(climate585) except: continue # calculate difference df = pd.merge(df_hist, df_245, left_on = 'reach_id', right_on = 'reach_id', suffixes = ['_hist','_245']) df = df.merge(df_585, left_on = 'reach_id', right_on = 'reach_id') df = df.rename(columns={'pred':'pred_585'}).set_index('reach_id') df['ch245'] = (df['pred_245'] - df['pred_hist']) / df['pred_hist'] * 100 df['ch585'] = (df['pred_585'] - df['pred_hist']) / df['pred_hist'] * 100 df['model'] = model df = df.loc[df.index.isin(basin_list),:] df_all.append(df) print(model) df_all = pd.concat(df_all).reset_index() # calculate ensemble mean and model agreement df33 = [] for a in ['hist','245','585']: df1 = df_all[['reach_id','pred_'+a,'model']].pivot_table(index = 'reach_id', columns = 'model', values = 'pred_'+a) ens = df1.mean(axis=1) df1['ens_'+a] = ens df33.append(df1[['ens_'+a]]) if a == 'hist': continue df0 = df_all[['reach_id','ch'+a,'model']].pivot_table(index = 'reach_id', columns = 'model', values = 'ch'+a) ens = df0.mean(axis=1) p = np.sum(df0.mul(ens, axis = 0) > 0, axis = 1) df0['diff_'+a] = ens df0['p_'+a] = p df33.append(df0[['diff_'+a,'p_'+a]]) df33 = pd.concat(df33, axis = 1) df33['global_id'] = df33.index.values df33.to_csv('../results/GRIT_climate_diff_ens_flood.csv') # calculate ensemble difference and model agreement between SSP245 and SSP585 df_all['change'] = (df_all['pred_585'] - df_all['pred_245']) / df_all['pred_245'] * 100 df_all = df_all[['reach_id','change','model']].pivot_table(index = 'reach_id', columns = 'model', values = 'change') ens = df_all.mean(axis = 1) p = np.sum(df_all.mul(ens, axis = 0) > 0, axis = 1) df_all = pd.DataFrame({'global_id':df_all.index.values, 'diff4585':ens, 'p4585':p}) df_all.to_csv('../results/GRIT_climate_245_585_diff_ens.csv', index = False) # add prediction results to gpkg df_pred = pd.read_csv('../results/GRIT_climate_diff_ens_flood.csv') df_change = pd.read_csv('../results/GRIT_climate_245_585_diff_ens.csv') for name in reachNames: gdf = read_dataframe(name) gdf = gdf.merge(df_pred, on = 'global_id') gdf = gdf.merge(df_change, on = 'global_id') write_dataframe(gdf, name[:-5]+'_RP10_diff.gpkg') #%% --------------------------------------------- Plotting # get file names reachNames = glob.glob('../data/GRIT/segments/GRITv05_segments_*_EPSG8857_correct_grwl0p5.gpkg') reachNames00 = [i[:-5]+'_RP10_diff.gpkg' for i in reachNames] if __name__ == '__main__': fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 10), subplot_kw={'projection':ccrs.EqualEarth()}) #%% ------------------------------------------------------------------------ plot map for colName, pName, ax0 in zip(['diff_245','diff4585'], ['p_245','p4585'], [ax1, ax2]): if '245' in colName: title = 'SSP2-4.5' clabel = 'Future change in 10-year flood magnitude (%)' # colorbar ticks bounds = np.array([-50, -40, -30, -20, -10, 10, 20, 30, 40, 50]) bounds_ticklabel = [str(i) for i in bounds] else: title = 'From SSP2-4.5 to SSP5-8.5' clabel = '( SSP5-8.5 - SSP2-4.5) / SSP2-4.5 (%)' # colorbar ticks bounds = np.array([-50, -40, -30, -20, -10, 10, 20, 30, 40, 50]) bounds_ticklabel = [str(i) for i in bounds] plot_grit(ax0, reachNames00, colName, pName, bounds, bounds_ticklabel, clabel, title = title, addLatMean = True, fontsize = 8) #%% ------------------------------------------------------------------------ plot 2 boxplot change for two scenarios 2071-2100 region_dict = { 'SP':'South Pacific\nand Australia', 'AF':'Africa', 'AS' : 'Asia\n(excl. Siberia)', 'EU' : 'Europe', 'NA' : 'North\nAmerica', 'SA' : 'South\nAmerica', 'SI' : 'Siberia'} tmp = [] for reachName in reachNames00: df0 = read_dataframe(reachName, read_geometry = False) region = os.path.basename(reachName).split('_')[2] df0['region'] = region_dict[region] tmp.append(df0) tmp0 = pd.concat(tmp) tmp = tmp0[['region','ens_245','ens_585']].copy() tmp['region'] = 'Global' tmp = pd.concat([tmp, tmp0[['region','ens_245','ens_585']]]).rename(columns={'ens_245':'SSP2-4.5','ens_585':'SSP5-8.5'}).melt(id_vars='region') ax = ax1.inset_axes([1.18, 0, .45, 1]) sns.boxplot(data = tmp, x = 'region', y = 'value', hue = 'variable', showfliers = False, ax = ax, palette = 'tab10', whis = (2.5, 97.5)) ax.set_xlabel(None) ax.axhline(y = 0, color = 'k', lw = .3, linestyle = '--') ax.tick_params(labelsize=8, labelrotation = 90) sns.move_legend(ax, title = None, frameon = False, loc = 'upper right', fontsize = 8) ax.set_ylabel('Change in 10-year flood magnitude (%)', fontsize = 8) ax.set_yscale('log') ax.set_ylim(None, 1e+6) #%% ------------------------------------------------------------------------ plot 3 change for three different periods and two scenarios tmp = [] for fname in glob.glob('../results/*/GRIT_*_flood_*_addLonLat.csv'): period = os.path.basename(fname).split('_') scenario, period = period[1], period[3] if period != '1981-2014' and scenario == 'historical': continue df0 = pd.read_csv(fname) df0['period'] = scenario + '_' + period tmp.append(df0) tmp = pd.concat(tmp) tmp = tmp.pivot_table(index = 'reach_id', columns = 'period', values = 'pred') change = tmp.iloc[:,1:].sub(tmp.iloc[:,0], axis = 0).div(tmp.iloc[:,0], axis = 0) * 100 change = change.reset_index().melt(id_vars = 'reach_id') ax = ax2.inset_axes([1.18, 0, .45, 1]) dict0 = {1:'2015-2040',2:'2041-2070',3:'2071-2100'} for name in ['ssp245','ssp585']: if name == 'ssp245': color_list = ['#8EA9F3', '#567DED', '#1D52E7'] else: color_list = ['#EA978E', '#DF6456', '#D4301D'] for color,period,ls in zip(color_list, ['2015-2040','2041-2070','2071-2100'], ['dashdot','dashed','solid']): data = change.loc[change.period == name+'_'+period,'value'].values # sort data data_sorted = np.sort(data) # Calculate proportions for CDF p = 1. * np.arange(len(data_sorted)) / (len(data_sorted) - 1) # subsampling data points to accelerate visualization idx = np.arange(0, len(p), 4) # Create CDF plot sns.lineplot(x=data_sorted[idx], y = p[idx], ax = ax, label = period + ' ' + name.upper(), color = color, ls = ls) print(name, period) ax.set_xlabel('Relative change (%)') ax.set_ylabel('CDF') ax.set_xlim(-20, 150) sns.move_legend(ax, "lower right", frameon = False) # add subplot ID fig.text(0.15, 0.88, 'a', weight = 'bold', size = 12, va = 'center', ha = 'center') fig.text(0.15, 0.47, 'b', weight = 'bold', size = 12, va = 'center', ha = 'center') fig.text(.99, 0.88, 'c', weight = 'bold', size = 12, va = 'center', ha = 'center') fig.text(.99, 0.47, 'd', weight = 'bold', size = 12, va = 'center', ha = 'center') fig.savefig('../picture/fig4.png', dpi = 600)