diff --git a/code/alpgm_interface.py b/code/alpgm_interface.py index b4def1b..2381ec5 100644 --- a/code/alpgm_interface.py +++ b/code/alpgm_interface.py @@ -32,7 +32,7 @@ ######### SETTINGS ############################################################################## # projection == True -> Projections with ADAMONT for the 21st century # projection == False -> Historical simulations for the 1984 - 2015 period with SAFRAN -historical_forcing, projection_forcing, simulation_type = settings.simulation_settings(projection = False) +historical_forcing, projection_forcing, simulation_type = settings.simulation_settings(projection = True) ### Global variables ### # Set the glacier index to start the simulations @@ -55,7 +55,7 @@ # SMB machine learning models generation settings.train_smb_model(historical_forcing, - compute_forcing = True, # Compute historical climate forcings + compute_forcing = False, # Compute historical climate forcings train_model = False) # Re-train SMB machine learning models ########## DELTA H FUNCTIONS GENERATION ####################### @@ -65,9 +65,9 @@ ########## SMB PROJECTION + GLACIER GEOMETRY EVOLUTION ####### settings.glacier_simulation(simulation_type, counter_threshold, - validate_SMB = True, # SMB model(s) validation or reconstruction + validate_SMB = False, # SMB model(s) validation or reconstruction compute_projection_forcings = False, # Compute projection climate forcings - compute_evolution = False, # Compute glacier evolution + compute_evolution = True, # Compute glacier evolution reconstruct = True, # Reconstruct glacier-wide SMB timeseries overwrite = True) diff --git a/code/glacier_evolution.py b/code/glacier_evolution.py index f1a1e62..b06904b 100644 --- a/code/glacier_evolution.py +++ b/code/glacier_evolution.py @@ -253,16 +253,16 @@ def preload_ensemble_SMB_models(): print("|", end="", flush=True) member_idx = member_idx+1 - member_idx = 0 - print("\n\nPreloading ensemble full SMB models...") - for path_member in path_ensemble_members: - # We retrieve the ensemble member ANN model - ann_member_model = load_model(os.path.join(path_ensemble, path_member, 'ann_glacier_model.h5'), custom_objects={"r2_keras": r2_keras, "root_mean_squared_error": root_mean_squared_error}, compile=False) -# ensemble_members.append(ann_member_model) - ensemble_members[member_idx] = ann_member_model - print("|", end="", flush=True) - member_idx = member_idx+1 - +# member_idx = 0 +# print("\n\nPreloading ensemble full SMB models...") +# for path_member in path_ensemble_members: +# # We retrieve the ensemble member ANN model +# ann_member_model = load_model(os.path.join(path_ensemble, path_member, 'ann_glacier_model.h5'), custom_objects={"r2_keras": r2_keras, "root_mean_squared_error": root_mean_squared_error}, compile=False) +## ensemble_members.append(ann_member_model) +# ensemble_members[member_idx] = ann_member_model +# print("|", end="", flush=True) +# member_idx = member_idx+1 +# CV_ensemble_members = np.asarray(CV_ensemble_members) ensemble_members = np.asarray(ensemble_members) print("\n") @@ -307,16 +307,16 @@ def make_ensemble_simulation(ensemble_SMB_models, x_ann, batch_size, glimsID, gl member_idx = 0 # We compute the quadratic slope difference - slope_dist = (training_slopes - ref_slope)**2 +# slope_dist = (training_slopes - ref_slope)**2 # Depending if glacier is present in training dataset we use the CV or full model - if(np.any(glims_rabatel['GLIMS_ID'] == glimsID.encode('utf-8'))): + if(not evolution and np.any(glims_rabatel['GLIMS_ID'] == glimsID.encode('utf-8'))): SMB_ensemble_members = ensemble_SMB_models['full'] print("\nFull ensemble models") else: - SMB_ensemble_members = ensemble_SMB_models['CV'] - CV_ensemble = True - print("\nCross-validation ensemble models") + SMB_ensemble_members = ensemble_SMB_models['CV'] + CV_ensemble = True + print("\nCross-validation ensemble models") # We iterate the previously loaded ensemble models for ensemble_model in SMB_ensemble_members: @@ -677,51 +677,6 @@ def get_SAFRAN_glacier_coordinates(glims_dataset): return all_glacier_coordinates -@jit -# Mean daily temperature -def get_mean_temps(datetimes, hourly_data): - ref_day = -9 - daily_data = [] - idx, day_idx = 0, 0 - first = True - for time_hour in datetimes: - current_day = time_hour.astype(object).timetuple().tm_yday - if(current_day == ref_day): - daily_data[day_idx] = daily_data[day_idx] + hourly_data[idx]/24.0 - else: - ref_day = current_day - daily_data.append(hourly_data[idx]/24.0) - if(not first): - day_idx = day_idx + 1 - else: - first = False - - idx = idx + 1 - - return np.asarray(daily_data) - -@jit -# Daily precipitations -def get_precips(datetimes, hourly_data): - ref_day = -9 - daily_data = [] - idx, day_idx = 0, 0 - isFirst = True - for time_hour in datetimes: - current_day = time_hour.astype(object).timetuple().tm_yday - if(current_day == ref_day): - daily_data[day_idx] = daily_data[day_idx] + hourly_data[idx] - else: - ref_day = current_day - daily_data.append(hourly_data[idx]) - if(not isFirst): - day_idx = day_idx + 1 - else: - isFirst = False - idx = idx + 1 - return np.asarray(daily_data) - - def compute_local_anomalies(glacier_CPDD, glacier_winter_snow, glacier_summer_snow, meteo_refs): local_CPDD_anomaly = glacier_CPDD - meteo_refs['CPDD'] @@ -838,7 +793,7 @@ def get_ADAMONT_idx(massif, glacier_altitude, massif_number, zs): return ADAMONT_idx -def get_default_ADAMONT_forcings(year_start, year_end, midfolder): +def get_default_ADAMONT_forcings(year_start, year_end, midfolder, overwrite): path_temps = os.path.join(path_smb_function_adamont, midfolder, 'daily_temps_years_' + str(year_start) + '-' + str(year_end) + '.txt') path_snow = os.path.join(path_smb_function_adamont, midfolder, 'daily_snow_years_' + str(year_start) + '-' + str(year_end) + '.txt') @@ -848,7 +803,7 @@ def get_default_ADAMONT_forcings(year_start, year_end, midfolder): path_massif = os.path.join(path_smb_function_adamont, midfolder, 'massif_number.txt') path_aspect = os.path.join(path_smb_function_adamont, midfolder, 'aspects.txt') - if(os.path.exists(path_temps) & os.path.exists(path_snow) & os.path.exists(path_rain) & os.path.exists(path_dates) & os.path.exists(path_zs) & os.path.exists(path_massif) & os.path.exists(path_aspect)): + if((os.path.exists(path_temps) & os.path.exists(path_snow) & os.path.exists(path_rain) & os.path.exists(path_dates) & os.path.exists(path_zs) & os.path.exists(path_massif) & os.path.exists(path_aspect))): print("Fetching ADAMONT forcings...") with open(path_temps, 'rb') as temps_f: daily_temps_years = np.load(temps_f, allow_pickle=True) @@ -872,27 +827,27 @@ def get_default_ADAMONT_forcings(year_start, year_end, midfolder): forcing_daymean = settings.current_ADAMONT_model_daymean forcing_daysum = settings.current_ADAMONT_model_daysum - print("\nCurrent ADAMONT combination: " + str(forcing_daymean) + "\\n") + print("\nCurrent ADAMONT combination: " + str(forcing_daymean) + "\n") start = time.time() # We load the two ADAMONT files - adamont_mean_climate = xr.open_dataset(forcing_daymean, parallel=True) - adamont_sum_climate = xr.open_dataset(forcing_daysum, parallel=True) + adamont_mean_climate = xr.open_dataset(forcing_daymean) + adamont_sum_climate = xr.open_dataset(forcing_daysum) - # Rename the time coordinates to match the SAFRAN format + # Rename the time coordinates to match the ADAMONT format adamont_mean_climate = adamont_mean_climate.rename({"TIME": "time"}) adamont_sum_climate = adamont_sum_climate.rename({"TIME": "time"}) end = time.time() - print("\n-> open SAFRAN dataset processing time: " + str(end - start) + " s") + print("\n-> open ADAMONT dataset processing time: " + str(end - start) + " s") zs = adamont_mean_climate['ZS'].compute() massif_number = adamont_mean_climate['MASSIF_NUMBER'].compute() aspects = np.repeat(-1, len(zs)) - daily_temps_years, daily_snow_years, daily_rain_years, daily_datetimes = np.array([]),np.array([]),np.array([]),np.array([]) + daily_temps_years, daily_snow_years, daily_rain_years, daily_datetimes = [],[],[],[] - for year in year_period: + for year in range(year_start, year_end+1): print("Hydrological year: " + str(year-1) + "-" + str(year)) start = time.time() @@ -903,16 +858,16 @@ def get_default_ADAMONT_forcings(year_start, year_end, midfolder): safran_rain_d = (adamont_sum_climate.sel(time = slice(str(year-1)+'-10-01', str(year)+'-09-30'))['RAIN'].resample(time="1D").sum()).compute() # Store daily raw data for future re-processing - daily_temps_years = np.append(daily_temps_years, safran_tmean_d) - daily_snow_years = np.append(daily_snow_years, safran_snow_d) - daily_rain_years = np.append(daily_rain_years, safran_rain_d) - daily_datetimes = np.append(daily_datetimes, safran_tmean_d.time) + daily_temps_years.append(safran_tmean_d.data) + daily_snow_years.append(safran_snow_d.data) + daily_rain_years.append(safran_rain_d.data) + daily_datetimes.append(safran_tmean_d.time.data) daily_meteo_data = {'temps':daily_temps_years, 'snow': daily_snow_years, 'rain': daily_rain_years, 'dates': daily_datetimes, 'zs': zs} # We create the folder if it's not there - if(not os.path.exists(path_smb_function_adamont+midfolder)): - os.makedirs(path_smb_function_adamont+midfolder) + if(not os.path.exists(os.path.join(path_smb_function_adamont, midfolder))): + os.makedirs(os.path.join(path_smb_function_adamont, midfolder)) with open(os.path.join(path_smb_function_adamont, midfolder, 'daily_temps_years_' + str(year_start) + '-' + str(year_end) + '.txt'), 'wb') as dtemp_f: np.save(dtemp_f, daily_temps_years) @@ -936,47 +891,61 @@ def get_adjusted_glacier_ADAMONT_forcings(year, year_start, glacier_mean_altitud idx = year - year_start # print("ADAMONT_idx: " + str(ADAMONT_idx)) glacier_idx = int(ADAMONT_idx) + t_lim = 0.0 - safran_tmean_d = daily_meteo_data['temps'][idx] - safran_snow_d = daily_meteo_data['snow'][idx] - safran_rain_d = daily_meteo_data['rain'][idx] - zs = daily_meteo_data['zs'][idx] - daily_datetimes_years = daily_meteo_data['dates'][idx] + # We create a dataset for the current year + adamont_ds = xr.Dataset({'temperature': (['time','space'], daily_meteo_data['temps'][idx]), + 'snow': (['time','space'], daily_meteo_data['snow'][idx]), + 'rain': (['time','space'], daily_meteo_data['rain'][idx])}, + coords={'zs': daily_meteo_data['zs'], + 'time': daily_meteo_data['dates'][idx]}) + # Re-scale temperature at glacier's actual altitude - safran_tmean_d_g = copy.deepcopy(safran_tmean_d[:, glacier_idx] + ((zs[0,glacier_idx].data - glacier_alt_y)/1000.0)*6.0) + adamont_tmean_d_g = copy.deepcopy(adamont_ds['temperature'][:, glacier_idx] + ((adamont_ds['zs'][glacier_idx].data - glacier_mean_altitude)/1000.0)*6.0) # We adjust the snowfall rate at the glacier's altitude - safran_snow_d_g = copy.deepcopy(safran_snow_d[:, glacier_idx]) - safran_rain_d_g = copy.deepcopy(safran_rain_d[:, glacier_idx]) - safran_snow_d_g.data = np.where(safran_tmean_d_g.data > 0.0, 0.0, safran_snow_d_g.data) - safran_snow_d_g.data = np.where(((safran_tmean_d_g.data < 0.0) & (safran_snow_d_g.data == 0.0)), safran_rain_d_g.data, safran_snow_d_g.data) + adamont_snow_d_g = copy.deepcopy(adamont_ds['snow'][:, glacier_idx]) + adamont_rain_d_g = copy.deepcopy(adamont_ds['rain'][:, glacier_idx]) + + adamont_snow_d_g.data = np.where(adamont_tmean_d_g.data > t_lim, 0.0, + adamont_snow_d_g.data) + adamont_snow_d_g.data = np.where((adamont_tmean_d_g.data < t_lim), adamont_snow_d_g.data + adamont_rain_d_g.data, + adamont_snow_d_g.data) # Monthly data during the current hydrological year - # Compute dask arrays prior to storage - safran_tmean_m_g = safran_tmean_d_g.resample(time="1MS").mean().data.compute() - safran_snow_m_g = safran_snow_d_g.resample(time="1MS").sum().data.compute() + adamont_tmean_m_g = adamont_tmean_d_g.resample(time="1MS").mean().data + adamont_snow_m_g = adamont_snow_d_g.resample(time="1MS").sum().data + + #### /!\ If climate data is not complete for year 2099, skip this year ####### + if(adamont_tmean_m_g.size != 12): + return {'CPDD': []}, {'temps': []} # Compute CPDD - # Compute dask arrays prior to storage - glacier_CPDD = np.sum(np.where(safran_tmean_d_g.data < 0, 0, safran_tmean_d_g.data)).compute() - glacier_CPDDs_all[j]['CPDD'].append(year) - glacier_CPDDs_all[j]['CPDD'].append(glacier_CPDD) + glacier_CPDD = np.sum(np.where(adamont_tmean_d_g.data < 0, 0, + adamont_tmean_d_g.data)) # Compute snowfall - # Compute dask arrays prior to storage - glacier_winter_snow = np.sum(safran_snow_d_g.sel(time = slice(str(year-1)+'-10-01', str(year)+'-03-31')).data).compute() - glacier_summer_snow = np.sum(safran_snow_d_g.sel(time = slice(str(year)+'-04-01', str(year)+'-07-31')).data).compute() + glacier_winter_snow = np.sum(adamont_snow_d_g.sel(time = slice(str(year-1)+'-10-01', str(year)+'-03-31')).data) + glacier_summer_snow = np.sum(adamont_snow_d_g.sel(time = slice(str(year)+'-04-01', str(year)+'-07-31')).data) # Seasonal anomalies - CPDD_LocalAnomaly, winter_snow_LocalAnomaly, summer_snow_LocalAnomaly = compute_local_anomalies(glacier_CPDD, glacier_winter_snow, glacier_summer_snow, + CPDD_LocalAnomaly, winter_snow_LocalAnomaly, summer_snow_LocalAnomaly = compute_local_anomalies(glacier_CPDD, + glacier_winter_snow, + glacier_summer_snow, meteo_anomalies) # Monthly anomalies - mon_temp_anomaly, mon_snow_anomaly = compute_monthly_anomalies(safran_tmean_m_g, safran_snow_m_g, meteo_anomalies['mon_temp'], meteo_anomalies['mon_snow']) + mon_temp_anomaly, mon_snow_anomaly = compute_monthly_anomalies(adamont_tmean_m_g, + adamont_snow_m_g, + meteo_anomalies['mon_temp'], + meteo_anomalies['mon_snow']) - season_anomalies_y = {'CPDD': CPDD_LocalAnomaly, 'winter_snow':winter_snow_LocalAnomaly, 'summer_snow': summer_snow_LocalAnomaly} - monthly_anomalies_y = {'temps': mon_temp_anomaly, 'snow': mon_snow_anomaly} + season_anomalies_y = {'CPDD': CPDD_LocalAnomaly, + 'winter_snow':winter_snow_LocalAnomaly, + 'summer_snow': summer_snow_LocalAnomaly} + monthly_anomalies_y = {'temps': mon_temp_anomaly, + 'snow': mon_snow_anomaly} # print("\nCPDD_LocalAnomaly: " + str(CPDD_LocalAnomaly)) # print("winter_snow_LocalAnomaly: " + str(winter_snow_LocalAnomaly)) @@ -1083,7 +1052,7 @@ def store_rasters(masked_DEM_current_glacier_u, masked_ID_current_glacier_u, mid if not os.path.exists(os.path.join(path_glacier_evolution_DEM_rasters, midfolder)): os.makedirs(os.path.join(path_glacier_evolution_DEM_rasters, midfolder)) path_ID_raster_year = os.path.join(path_glacier_evolution_ID_rasters, midfolder, "IceDepth_Glacier_0" + str(glacierID) + "_" + str(year) + ".tif") - if not os.path.exists(path_glacier_evolution_ID_rasters + midfolder): + if not os.path.exists(os.path.join(path_glacier_evolution_ID_rasters, midfolder)): os.makedirs(os.path.join(path_glacier_evolution_ID_rasters, midfolder)) array2raster(path_DEM_raster_year, r_origin, r_pixelwidth, r_pixelheight, masked_DEM_current_glacier_u) array2raster(path_ID_raster_year, r_origin, r_pixelwidth, r_pixelheight, masked_ID_current_glacier_u) @@ -1167,6 +1136,10 @@ def glacier_evolution(masked_DEM_current_glacier, masked_ID_current_glacier, season_anomalies_y, monthly_anomalies_y = get_adjusted_glacier_ADAMONT_forcings(year, year_start, masked_DEM_current_glacier_u.compressed().mean(), SAFRAN_idx, daily_meteo_data, meteo_anomalies) + ### Check if data is not available for year 2099 ### + if(len(monthly_anomalies_y['temps']) == 0): + # Skip year + break #### CREATION OF THE MODEL TEST DATASET #### x_lasso, x_ann = create_input_array(season_anomalies_y, monthly_anomalies_y, mean_glacier_alt, max_glacier_alt, slope20, current_glacierArea, lat, lon, aspect) @@ -1317,10 +1290,9 @@ def main(compute, ensemble_SMB_models, overwrite_flag, use_cluster, counter_thre path_smb = os.path.join(workspace, 'glacier_data', 'smb') path_smb_function = os.path.join(path_smb, 'smb_function') path_glacier_outlines_shapefile = os.path.join(path_glacier_2003_shapefiles, 'GLIMS_glaciers_2003_ID_massif' + '.shp') + # SAFRAN climate forcings path_ann = settings.path_ann path_safran_forcings = os.path.join(path_smb_function, 'SAFRAN') - # SAFRAN climate forcings - path_safran_forcings = settings.path_safran # Path to be updated with ADAMONT forcings local path path_adamont_forcings = settings.path_adamont @@ -1340,7 +1312,7 @@ def main(compute, ensemble_SMB_models, overwrite_flag, use_cluster, counter_thre # glims_2015 = genfromtxt(path_glims + 'GLIMS_2015_massif.csv', delimiter=';', skip_header=1, dtype=[('Area', '= 1 and fold < 5): -# test_matrix[-1, :] = 1 -# # Fill train matrix -# train_matrix[-1, :] = 0 -# if(fold >= 1 and fold < 11): -## test_matrix[2, :] = 1 -# test_matrix[:, :27] = 1 -# # Fill train matrix -## train_matrix[2, :] = 0 -# train_matrix[:, :27] = 0 -## else: - # Fill test matrix - test_matrix[glacier_idx, :] = 1 - # Fill train matrix - train_matrix[glacier_idx, :] = 0 -# - - # Fill test matrix -# test_matrix[random_glaciers[glacier_idx], :] = 1 -# test_matrix[random_glaciers[glacier_idx+1], :] = 1 -# if(not(fold >= 1 and fold < 11)): - test_matrix[:, year_idx] = 1 -# test_matrix[:, random_years[year_idx+1]] = 1 - -# # 1 -# test_matrix[random_glaciers[glacier_idx], random_years[year_idx]] = 1 -# test_matrix[random_glaciers[glacier_idx], random_years[year_idx+1]] = 1 -# test_matrix[random_glaciers[glacier_idx], random_years[year_idx+2]] = 1 -# test_matrix[random_glaciers[glacier_idx], random_years[year_idx+3]] = 1 -# -# # 2 -# test_matrix[random_glaciers[glacier_idx+1], random_years[year_idx]] = 1 -# test_matrix[random_glaciers[glacier_idx+1], random_years[year_idx+1]] = 1 -# test_matrix[random_glaciers[glacier_idx+1], random_years[year_idx+2]] = 1 -# test_matrix[random_glaciers[glacier_idx+1], random_years[year_idx+3]] = 1 -# -# # 3 -# test_matrix[random_glaciers[glacier_idx+2], random_years[year_idx]] = 1 -# test_matrix[random_glaciers[glacier_idx+2], random_years[year_idx+1]] = 1 -# test_matrix[random_glaciers[glacier_idx+2], random_years[year_idx+2]] = 1 -# test_matrix[random_glaciers[glacier_idx+2], random_years[year_idx+3]] = 1 -# -# # 4 -# test_matrix[random_glaciers[glacier_idx+3], random_years[year_idx]] = 1 -# test_matrix[random_glaciers[glacier_idx+3], random_years[year_idx+1]] = 1 -# test_matrix[random_glaciers[glacier_idx+3], random_years[year_idx+2]] = 1 -# test_matrix[random_glaciers[glacier_idx+3], random_years[year_idx+3]] = 1 - - lsygo_test_matrixes.append(test_matrix) + if(lsygo_soft): + ################## Soft LSYGO ################################### + choice = weighted_choice(np.asarray(range(0,1048)), p_weights) + + # Weighted bagging + idx, i, j = 0, 0, 0 + for i in range(0, y_o.shape[0]): + for j in range(0, y_o.shape[1]): + if(np.isfinite(y_o[i,j])): + if(choice == idx): + glacier_idx = i + year_idx = j + # print("counter = " + str(idx)) + idx=idx+1 + + print("\nChosen glacier: " + str(glacier_idx)) + print("Chosen year: " + str(year_idx)) + print("SMB: " + str(y_o[glacier_idx,year_idx])) + avg_sampled_smb = np.concatenate((avg_sampled_smb, np.concatenate((y_o[glacier_idx,:], y_o[:,year_idx])))) + + # Fill test matrix + test_matrix[glacier_idx, :] = 1 + # Fill train matrix + train_matrix[glacier_idx, :] = 0 + + # Fill test matrix + test_matrix[:, year_idx] = 1 + + # Fill train matrix + train_matrix[:, year_idx] = 0 - # Fill train matrix -# train_matrix[random_glaciers[glacier_idx], :] = 0 -# train_matrix[random_glaciers[glacier_idx+1], :] = 0 -# train_matrix[random_glaciers[glacier_idx+2], :] = 0 -# train_matrix[random_glaciers[glacier_idx+3], :] = 0 -# if(not(fold >= 1 and fold < 11)): - train_matrix[:, year_idx] = 0 -# train_matrix[:, random_years[year_idx+1]] = 0 -# train_matrix[:, random_years[year_idx+2]] = 0 -# train_matrix[:, random_years[year_idx+3]] = 0 + else: + ############ Hard LSYGO ########################################## + +# print("\nglacier_idx: " + str(glacier_idx)) +# print("year_idx: " + str(year_idx)) +# print("random_glaciers[glacier_idx]: " + str(random_glaciers[glacier_idx])) +# print("random_years[year_idx]: " + str(random_years[year_idx])) + + # Fill test matrix + # 1 + test_matrix[random_glaciers[glacier_idx], random_years[year_idx]] = 1 + test_matrix[random_glaciers[glacier_idx], random_years[year_idx+1]] = 1 + test_matrix[random_glaciers[glacier_idx], random_years[year_idx+2]] = 1 + test_matrix[random_glaciers[glacier_idx], random_years[year_idx+3]] = 1 + + # 2 + test_matrix[random_glaciers[glacier_idx+1], random_years[year_idx]] = 1 + test_matrix[random_glaciers[glacier_idx+1], random_years[year_idx+1]] = 1 + test_matrix[random_glaciers[glacier_idx+1], random_years[year_idx+2]] = 1 + test_matrix[random_glaciers[glacier_idx+1], random_years[year_idx+3]] = 1 + + # 3 + test_matrix[random_glaciers[glacier_idx+2], random_years[year_idx]] = 1 + test_matrix[random_glaciers[glacier_idx+2], random_years[year_idx+1]] = 1 + test_matrix[random_glaciers[glacier_idx+2], random_years[year_idx+2]] = 1 + test_matrix[random_glaciers[glacier_idx+2], random_years[year_idx+3]] = 1 + + # 4 + test_matrix[random_glaciers[glacier_idx+3], random_years[year_idx]] = 1 + test_matrix[random_glaciers[glacier_idx+3], random_years[year_idx+1]] = 1 + test_matrix[random_glaciers[glacier_idx+3], random_years[year_idx+2]] = 1 + test_matrix[random_glaciers[glacier_idx+3], random_years[year_idx+3]] = 1 + + # Fill train matrix + train_matrix[random_glaciers[glacier_idx], :] = 0 + train_matrix[random_glaciers[glacier_idx+1], :] = 0 + train_matrix[random_glaciers[glacier_idx+2], :] = 0 + train_matrix[random_glaciers[glacier_idx+3], :] = 0 + + train_matrix[:, random_years[year_idx]] = 0 + train_matrix[:, random_years[year_idx+1]] = 0 + train_matrix[:, random_years[year_idx+2]] = 0 + train_matrix[:, random_years[year_idx+3]] = 0 + + ########################################################################### + # Add matrixes to folds + lsygo_test_matrixes.append(test_matrix) lsygo_train_matrixes.append(train_matrix) year_idx = year_idx+4 glacier_idx = glacier_idx+4 +#import pdb; pdb.set_trace() + # We capture the mask from the SMB data to remove the nan gaps finite_mask = np.isfinite(y) -print("\nAverage sampled SMB LSYGO: " + str(np.nanmean(avg_sampled_smb))) +if(lsygo_soft): + print("\nAverage sampled SMB LSYGO: " + str(np.nanmean(avg_sampled_smb))) X = X[finite_mask,:] y = y[finite_mask] @@ -530,8 +524,10 @@ def create_lsygo_model(n_features, final): lsygo_test_int_folds, lsygo_train_int_folds = [],[] # Filter LSYGO folds for test_fold, train_fold in zip(lsygo_test_matrixes, lsygo_train_matrixes): - lsygo_test_int_folds.append(test_fold.flatten()[finite_mask]) - lsygo_train_int_folds.append(train_fold.flatten()[finite_mask]) +# print("test_fold.flatten()[finite_mask]: " + str(np.sum(test_fold.flatten()[finite_mask]))) + if(np.sum(test_fold.flatten()[finite_mask]) > 0): + lsygo_test_int_folds.append(test_fold.flatten()[finite_mask]) + lsygo_train_int_folds.append(train_fold.flatten()[finite_mask]) # From int to boolean lsygo_test_folds = np.array(lsygo_test_int_folds, dtype=bool) @@ -730,6 +726,7 @@ def create_lsygo_model(n_features, final): SMB_nn_all = [] RMSE_nn_all, RMSE_nn_all_w = [],[] + bias_nn_all = [] r2_nn_all, r2_nn_all_w = [],[] SMB_ref_all = [] fold_idx = 1 @@ -884,6 +881,7 @@ def create_lsygo_model(n_features, final): r2_nn_all = np.concatenate((r2_nn_all, r2_score(y_test, SMB_nn)), axis=None) RMSE_nn_all = np.concatenate((RMSE_nn_all, math.sqrt(mean_squared_error(y_test, SMB_nn))), axis=None) + bias_nn_all = np.concatenate((bias_nn_all, np.average(y_test - SMB_nn)), axis=None) SMB_ref_all = np.concatenate((SMB_ref_all, y_test), axis=None) if(w_weights): @@ -900,6 +898,7 @@ def create_lsygo_model(n_features, final): r2_nn_all = np.asarray(r2_nn_all) RMSE_nn_all = np.asarray(RMSE_nn_all) + bias_nn_all = np.asarray(bias_nn_all) if(w_weights): r2_nn_all_w = np.asarray(r2_nn_all_w) RMSE_nn_all_w = np.asarray(RMSE_nn_all_w) @@ -924,11 +923,14 @@ def create_lsygo_model(n_features, final): print("\nMean overall RMSE w/ weights: " + str(math.sqrt(mean_squared_error(SMB_ref_all, SMB_nn_all, weights_validation)))) print("\nRMSE per fold: " + str(RMSE_nn_all)) + print("\nBias per fold: " + str(bias_nn_all)) print("\nAverage bias: " + str(np.average(SMB_nn_all - SMB_ref_all))) with open(os.path.join(path_ann, 'RMSE_per_fold.txt'), 'wb') as rmse_f: np.save(rmse_f, RMSE_nn_all) + with open(os.path.join(path_ann, 'bias_per_fold.txt'), 'wb') as bias_f: + np.save(bias_f, bias_nn_all) # import pdb; pdb.set_trace() diff --git a/scripts/plot_projections.py b/scripts/plot_projections.py index 3f16c70..cfe77f5 100644 --- a/scripts/plot_projections.py +++ b/scripts/plot_projections.py @@ -34,36 +34,36 @@ ###### FILE PATHS ####### # Folders -workspace = str(Path(os.getcwd()).parent) + '\\' -path_glims = workspace + 'glacier_data\\GLIMS\\' +workspace = str(Path(os.getcwd()).parent) +path_glims = os.path.join(workspace, 'glacier_data', 'GLIMS') #path_obs = 'C:\\Jordi\\PhD\\Data\\Obs\\' -path_smb = workspace + 'glacier_data\\smb\\' -path_glacier_evolution = workspace + 'glacier_data\\glacier_evolution\\' +path_smb = os.path.join(workspace, 'glacier_data', 'smb') +path_glacier_evolution = os.path.join(workspace, 'glacier_data', 'glacier_evolution') #path_glacier_evolution = 'C:\\Jordi\\PhD\\Simulations\\AGU 2018\\' -path_smb_simulations = path_smb + 'smb_simulations\\' -path_glacier_evolution_plots = path_glacier_evolution + 'plots\\' -path_glacier_area = path_glacier_evolution + 'glacier_area\\' -path_glacier_volume = path_glacier_evolution + 'glacier_volume\\' -path_glacier_zmean = path_glacier_evolution + 'glacier_zmean\\' -path_glacier_slope20 = path_glacier_evolution + 'glacier_slope20\\' -path_glacier_melt_years = path_glacier_evolution + 'glacier_melt_years\\' -path_glacier_w_errors = path_glacier_evolution + 'glacier_w_errors\\' -path_glacier_CPDDs = path_glacier_evolution + 'glacier_CPDDs\\' -path_glacier_snowfall = path_glacier_evolution + 'glacier_snowfall\\' -path_glacier_discharge = path_glacier_evolution + 'glacier_meltwater_discharge\\' +path_smb_simulations = os.path.join(path_smb, 'smb_simulations') +path_glacier_evolution_plots = os.path.join(path_glacier_evolution, 'plots') +path_glacier_area = os.path.join(path_glacier_evolution, 'glacier_area') +path_glacier_volume = os.path.join(path_glacier_evolution, 'glacier_volume') +path_glacier_zmean = os.path.join(path_glacier_evolution, 'glacier_zmean') +path_glacier_slope20 = os.path.join(path_glacier_evolution, 'glacier_slope20') +path_glacier_melt_years = os.path.join(path_glacier_evolution, 'glacier_melt_years') +path_glacier_w_errors = os.path.join(path_glacier_evolution, 'glacier_w_errors') +path_glacier_CPDDs = os.path.join(path_glacier_evolution, 'glacier_CPDDs') +path_glacier_snowfall = os.path.join(path_glacier_evolution, 'glacier_snowfall') +path_glacier_discharge = os.path.join(path_glacier_evolution, 'glacier_meltwater_discharge') #### We fetch all the data from the simulations -path_area_root = np.asarray(os.listdir(path_glacier_area + "projections\\")) -path_melt_years_root = np.asarray(os.listdir(path_glacier_melt_years + "projections\\")) -path_slope20_root = np.asarray(os.listdir(path_glacier_slope20 + "projections\\")) -path_volume_root = np.asarray(os.listdir(path_glacier_volume + "projections\\")) -path_errors_root = np.asarray(os.listdir(path_glacier_w_errors + "projections\\")) -path_zmean_root = np.asarray(os.listdir(path_glacier_zmean + "projections\\")) -path_CPDDs_root = np.asarray(os.listdir(path_glacier_CPDDs + "projections\\")) -path_snowfall_root = np.asarray(os.listdir(path_glacier_snowfall + "projections\\")) -path_SMB_root = np.asarray(os.listdir(path_smb_simulations + "projections\\")) - -glims_2003 = genfromtxt(path_glims + 'GLIMS_2003.csv', delimiter=';', skip_header=1, dtype=[('Area', ' 1): for year in range(year_start, 2100): diff --git a/scripts/plot_smb_results.py b/scripts/plot_smb_results.py index b4e6c6d..78f2600 100644 --- a/scripts/plot_smb_results.py +++ b/scripts/plot_smb_results.py @@ -8,6 +8,7 @@ ## Dependencies: ## import matplotlib.pyplot as plt import matplotlib as mpl +import scipy.stats as st import numpy as np from numpy import genfromtxt import os @@ -86,7 +87,7 @@ glims_rabatel = genfromtxt(path_glims + 'GLIMS_Rabatel_30_2003.csv', delimiter=';', skip_header=1, dtype=[('Area', ' -60): + if((np.isfinite(smb_glacier[-1]) and np.sum(smb_glacier[:-15]) < 0) and np.sum(smb_glacier[:13]) > -15): # if(True): line1, = ax11.plot(range(1967, 2016), smb_glacier, linewidth=linewidth, alpha=alpha) line2, = ax21.plot(range(1967, 2016), np.cumsum(smb_glacier), linewidth=linewidth, alpha=alpha) + if(smb_glacier.size > 40): + all_cum_smb.append(np.sum(smb_glacier)) big_glacier, small_glacier, very_small_glacier, extremely_small_glacier = False, False, False, False s_0_20, s_20_30, s_30_40, s_40_50 = False, False, False, False @@ -495,6 +499,7 @@ # All glaciers all_glacier_smb = np.asarray(all_glacier_smb) +all_cum_smb = np.asarray(all_cum_smb) flat_all_glacier_smb = np.asarray(flat_all_glacier_smb) a_avg_smb, a_median_smb, a_avg_smb_marzeion, a_median_smb_marzeion = [],[],[],[] annual_avg_smb_area = {'big':[], 'small':[], 'v_small':[], 'x_small':[]} @@ -625,8 +630,9 @@ ax21.axhline(y=0, color='black', linewidth=0.9, linestyle='-') line111, = ax11.plot(range(1967, 2016), a_avg_smb, linewidth=3, c='midnightblue', label='Area weighted mean') line112, = ax21.plot(range(1967, 2016), np.cumsum(a_avg_smb), linewidth=3, c='midnightblue', label='Area weighted mean') -ax11.legend(fontsize='large') -ax21.legend(fontsize='large') +ylim = ax21.get_ylim() +ax11.legend(fontsize=16) +ax21.legend(fontsize=16) fig11.tight_layout() for member in a_avg_ensemble_smb: @@ -643,6 +649,32 @@ #ax1.legend(fontsize='x-large') ax2.legend(fontsize='x-large') +fig1.tight_layout() + + +## Histogram for cumulative SMB #### +#ax_hist = plt.subplot(1,1,1) +#plt.hist(all_cum_smb, density=True, bins=100, alpha=0.7) +#plt.xlim(ylim) +#kde_xs = np.linspace(all_cum_smb.min(), all_cum_smb.max(), 301) +#kde = st.gaussian_kde(all_cum_smb) +#plt.plot(kde_xs, kde.pdf(kde_xs), color='darkgoldenrod', linewidth='3', label="PDF") +#plt.axvline(x=np.sum(a_avg_smb), color='midnightblue', linestyle='--', linewidth='3', label="Weighted average") +#plt.gca().invert_xaxis() +##plt.tick_params( +## axis='x', # changes apply to the x-axis +## which='both', # both major and minor ticks are affected +## bottom=False, # ticks along the bottom edge are off +## top=False, # ticks along the top edge are off +## labelbottom=False) # labels along the bottom edge are off +#ax_hist.yaxis.tick_right() +#plt.yticks(rotation='vertical') +#plt.legend() +#plt.ylabel('Probability') +##plt.xlabel('Cumulative SMB (m w.e.)') +##plt.title("Histogram"); + +#import pdb; pdb.set_trace() print("\nNumber of glaciers disappeared between 2003 and 2015: " + str(glaciers_not_2015)) @@ -669,20 +701,21 @@ ax47.set_ylabel('Cumulative area influence on glacier-wide SMB signal (m.w.e.)', fontsize=14) ax47.set_xlabel('Year', fontsize=14) ax47.tick_params(labelsize=14) -line25, = ax47.plot(range(1967, 2016)[-32:], a_median_smb[-32:], linestyle=':', linewidth=1, label='B: Reconstructed SMB', c='darkred') -line25, = ax47.plot(range(1967, 2016)[-32:], a_median_smb_marzeion[-32:], linestyle=':', linewidth=1, label='M: Reconstructed SMB', c='midnightblue') line25, = ax47.plot(range(1967, 2016)[-32:], np.cumsum(annual_avg_smb_obs[-32:] - a_median_smb[-32:]), linestyle='--', linewidth=2, label='Observations', c='olivedrab') +line25, = ax47.plot(range(1967, 2016)[-32:], a_median_smb[-32:], linestyle=':', linewidth=1, label='B: Reconstructed SMB', c='darkred') line14, = ax47.plot(range(1967, 2016)[-32:], np.cumsum(annual_avg_smb_area['x_small'][-32:] - a_median_smb[-32:]), linewidth=2, label='B: Glaciers < 0.1 km$^2$', c='darkred') line24, = ax47.plot(range(1967, 2016)[-32:], np.cumsum(annual_avg_smb_area['v_small'][-32:] - a_median_smb[-32:]), linewidth=2, label='B: Glaciers 0.1 - 0.5 km$^2$', c='crimson') line23, = ax47.plot(range(1967, 2016)[-32:], np.cumsum(annual_avg_smb_area['small'][-32:] - a_median_smb[-32:]), linewidth=2, label='B: Glaciers 0.5 - 2 km$^2$', c='darkorange') line22, = ax47.plot(range(1967, 2016)[-32:], np.cumsum(annual_avg_smb_area['big'][-32:] - a_median_smb[-32:]), linewidth=2, label='B: Glaciers > 2 km$^2$', c='tan') +line25, = ax47.plot(range(1967, 2016)[-32:], a_median_smb_marzeion[-32:], linestyle=':', linewidth=1, label='M: Reconstructed SMB', c='midnightblue') line14, = ax47.plot(range(1967, 2016)[-32:], np.cumsum(marzeion_annual_avg_smb_area['x_small'][-32:] - a_median_smb_marzeion[-32:]), linewidth=2, label='M: Glaciers < 0.1 km$^2$', c='midnightblue') line14, = ax47.plot(range(1967, 2016)[-32:], np.cumsum(marzeion_annual_avg_smb_area['v_small'][-32:] - a_median_smb_marzeion[-32:]), linewidth=2, label='M: Glaciers 0.1 - 0.5 km$^2$', c='slateblue') line13, = ax47.plot(range(1967, 2016)[-32:], np.cumsum(marzeion_annual_avg_smb_area['small'][-32:] - a_median_smb_marzeion[-32:]), linewidth=2, label='M: Glaciers 0.5 - 2 km$^2$', c='mediumorchid') line12, = ax47.plot(range(1967, 2016)[-32:], np.cumsum(marzeion_annual_avg_smb_area['big'][-32:] - a_median_smb_marzeion[-32:]), linewidth=2, label='M: Glaciers > 2 km$^2$', c='thistle') + #ax47.legend(loc='upper center', bbox_to_anchor=[-0.1, 1.25], ncol=4) -ax47.legend(loc='upper center', bbox_to_anchor=[0.5, 1.2], ncol=3, fontsize='large') +ax47.legend(ncol=2, fontsize='large') plt.subplots_adjust(top=0.80) @@ -1057,6 +1090,8 @@ avg_decadal_smb_g = np.array([avg_smb_70s_g, avg_smb_80s_g, avg_smb_90s_g, avg_smb_00s_g, avg_smb_10s_g]) xmin = np.array([1970, 1980, 1990, 2000, 2010]) xmax = np.array([1980, 1990, 2000, 2010, 2015]) +x_center = np.array([1975, 1985, 1995, 2005, 2012.5]) +decadal_uncertainties = np.array([0.18, 0.16, 0.14, 0.16, 0.11]) total_avg_smb_g = np.average(avg_glacier_smb, weights=area_glaciers) fig9, ax9 = plt.subplots(figsize=(6, 4)) @@ -1071,11 +1106,14 @@ ax9.hlines(total_avg_smb_g, 1967, 2015, color='darkblue', linewidth=7, label='Total average SMB') #ax9.hlines(total_avg_smb_marzeion, 1967, 2015, color='darkred', linewidth=4, label='Total average SMB (update of Marzeion et al., 2015)') #ax9.hlines(avg_decadal_smb_g, xmin, xmax, color='steelblue', linewidth=6, label='Decadal average SMB (this study)') +ax9.errorbar(x_center, avg_decadal_smb_g, yerr=decadal_uncertainties, fmt='none', elinewidth=1, ecolor='midnightblue') ax9.hlines(avg_decadal_smb_g, xmin, xmax, color='steelblue', linewidth=6, label='Decadal average SMB') +#import pdb; pdb.set_trace() + #ax9.hlines(avg_decadal_smb_marzeion, xmin, xmax, color='darkgoldenrod', linewidth=6, label='Decadal average SMB (update of Marzeion et al., 2015)') ax9.set_xticks([1970, 1980, 1990, 2000, 2010, 2015]) ax9.set_xticklabels(('1970', '1980', '1990', '2000', '2010', '2015')) -ax9.xaxis.grid(True) +ax9.xaxis.grid(True, linewidth=2) #ax9.legend(loc='lower left', mode= 'expand', bbox_to_anchor=(0,1.02,1,0.2)) #plt.subplots_adjust(top=0.80, left=0.15) ax9.legend() @@ -1113,8 +1151,8 @@ fig10, (ax10, ax11) = plt.subplots(2,1, figsize=(6,12)) #fig10.suptitle("Annual glacier-wide SMB by massif") ax10.axhline(y=0, color='black', linewidth=0.7, linestyle='-') -ax10.set_ylabel('Glacier-wide SMB (m.w.e. $a^{-1}$)', fontsize=13) -ax10.tick_params(labelsize=11) +ax10.set_ylabel('Glacier-wide SMB (m.w.e. $a^{-1}$)', fontsize=15) +ax10.tick_params(labelsize=14) #ax10.set_xlabel('Year') avg_smb_per_massif = copy.deepcopy(smb_massif_template) @@ -1132,9 +1170,9 @@ line101, = ax10.plot(range(1967, 2016), avg_smb_massif[massif], color=color, marker=marker, linewidth=1, label=massif) ax11.axhline(y=0, color='black', linewidth=0.7, linestyle='-') -ax11.set_ylabel('Cumulative glacier-wide SMB (m.w.e.)', fontsize=13) -ax11.set_xlabel('Year', fontsize=13) -ax11.tick_params(labelsize=11) +ax11.set_ylabel('Cumulative glacier-wide SMB (m.w.e.)', fontsize=15) +ax11.set_xlabel('Year', fontsize=14) +ax11.tick_params(labelsize=14) for massif, color, marker in zip(avg_smb_massif, colors, markers): format_str = "{color}{marker}-".format(color=color, marker=marker) line111, = ax11.plot(range(1967, 2016), np.cumsum(avg_smb_massif[massif]), color=color, marker=marker, linewidth=1, label=massif) @@ -1142,7 +1180,7 @@ handles, labels = ax11.get_legend_handles_labels() #ax11.legend(handles, labels, loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, shadow=True, ncol=5) #ax11.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=7, mode="expand", borderaxespad=0.) -ax11.legend(prop={'size': 12}) +ax11.legend(prop={'size': 14}) plt.tight_layout() plt.subplots_adjust(hspace = 0.08) diff --git a/scripts/smb_reconstructions_validation.py b/scripts/smb_reconstructions_validation.py index 6cfb924..f288d9c 100644 --- a/scripts/smb_reconstructions_validation.py +++ b/scripts/smb_reconstructions_validation.py @@ -151,8 +151,6 @@ def plot_and_save(glacier_info, glacier_annual_obs, glacier_annual_reconstructio print("\nDavaze et al. (2020) SMB 2003-2012: " + str(glacier_annual_obs[-13:-4])) print("Avg SMB: " + str(np.mean(glacier_annual_obs[-13:-4]))) - print("\nBolibar et al. (2020) avg SMB 2003-2012: " + str(glacier_annual_reconstructions[-13:-4].mean())) - # We plot the comparison and store the individual plots plot_and_save(glacier_info, glacier_annual_obs, glacier_annual_reconstructions, plot_count) @@ -160,6 +158,8 @@ def plot_and_save(glacier_info, glacier_annual_obs, glacier_annual_reconstructio else: glacier_metrics['RMSE'].append([]) + + print("\nBolibar et al. (2020) avg SMB 2003-2012: " + str(glacier_annual_reconstructions[-13:-4].mean())) print("\n----------------------------------------------") @@ -179,7 +179,7 @@ def plot_and_save(glacier_info, glacier_annual_obs, glacier_annual_reconstructio plt.show() - #import pdb; pdb.set_trace() +#import pdb; pdb.set_trace() ##### FINAL WHISKERS PLOTS ############### @@ -189,8 +189,8 @@ def plot_and_save(glacier_info, glacier_annual_obs, glacier_annual_reconstructio SMB_2000_2015 = genfromtxt(os.path.join(path_smb_validation, 'SMB_2000_2015.csv'), delimiter=';', dtype=None) SMB_2003_2012 = genfromtxt(os.path.join(path_smb_validation, 'SMB_2003_2012.csv'), delimiter=';', dtype=None) -uncertainties_2000_2015 = np.array([0.27, 0.4, 0.55]) -uncertainties_2003_2012 = np.array([0.27, 0.4, 0.55, 0.2]) +uncertainties_2000_2015 = np.array([0.27, 0.16, 0.55]) +uncertainties_2003_2012 = np.array([0.27, 0.55, 0.55, 0.2]) #import pdb; pdb.set_trace() @@ -213,8 +213,7 @@ def plot_and_save(glacier_info, glacier_annual_obs, glacier_annual_reconstructio m_colors = np.array(['golden brown', 'sienna', 'slate blue']) #e_colors = np.array(['palegreen', 'steelblue', 'wheat']) - -f1, axs1 = plot.subplots(ncols=4, nrows=3, share=2, tight=True, axwidth=1) +f1, axs1 = plot.subplots(ncols=4, nrows=4, share=2, tight=True, axwidth=1) for idx in range(0, raw_SMB_2000_2015.shape[0]): if(raw_SMB_2000_2015[idx,:].size > 0): axs1[idx].axhline(y=0, color='black', linewidth=0.5, linestyle='-', zorder=0) @@ -227,9 +226,9 @@ def plot_and_save(glacier_info, glacier_annual_obs, glacier_annual_reconstructio #import pdb; pdb.set_trace() axs1.format( - suptitle='2000-2015 annual glacier-wide SMB bias', + suptitle='2000-2015 annual glacier-wide SMB', abc=True, abcstyle='a', abcloc='ul', abcborder=True, - ylabel='Bias (m w.e. $a^{-1}$)', + ylabel='Glacier-wide SMB (m w.e. $a^{-1}$)', # xlim=(1, 10), xticks=1, ylim=(-2.5, 0.25), yticks=plot.arange(-2.5, 0.25), # xtickdir='inout', xtickminor=False, @@ -271,9 +270,9 @@ def plot_and_save(glacier_info, glacier_annual_obs, glacier_annual_reconstructio #import pdb; pdb.set_trace() axs2.format( - suptitle='2003-2012 annual glacier-wide SMB bias', + suptitle='2003-2012 glacier-wide SMB', abc=True, abcstyle='a', abcloc='ul', abcborder=True, - ylabel='Bias (m w.e. $a^{-1}$)', + ylabel='Glacier-wide SMB (m w.e. $a^{-1}$)', # xlim=(1, 10), xticks=1, ylim=(-2.75, 0.25), yticks=plot.arange(-2.75, 0.25), # xtickdir='inout', xtickminor=False, @@ -292,7 +291,7 @@ def plot_and_save(glacier_info, glacier_annual_obs, glacier_annual_reconstructio #plt.show() -import pdb; pdb.set_trace() +#import pdb; pdb.set_trace()