saving most recent PyGEM version before reverting to make a release for the HMA papers

Author: drounce

Farinotti_icethickness_process.ipynb

@@ -51,7 +51,7 @@
 option_watersheds_colored = 0                       # still good 11/6/2019
 option_runoff_monthlychange_and_components = 0      # updated - 11/20/2019
 runoff_erainterim_bywatershed = 0                   # updated - better to export to table
-option_excess_meltwater_diagram = 1
+option_excess_meltwater_diagram = 0
 
 option_startdate = 0
 
@@ -61,7 +61,7 @@
 option_glaciermip_table = 0                         # updated - 11/12/2019
 option_zemp_compare = 0                             # updated - 11/6/2019
 option_gardelle_compare = 0                         # updated - 11/6/2019
-option_wgms_compare = 0                             # updated - 11/6/2019
+option_wgms_compare = 1                             # updated - 11/6/2019
 option_dehecq_compare = 0
 option_uncertainty_fig = 0                          # updated - 11/12/2019
 option_nick_snowline = 0
@@ -5776,8 +5776,8 @@ def __call__(self, value, clip=None):
     # Loop through conditions:
     condition_dict = OrderedDict()
     condition_dict['All']= cal_data.index.values
-    condition_dict['All (> 1 yr)'] = (cal_data['t2'] - cal_data['t1']) >= 0.75
-    condition_dict['All (< 1 yr)'] = (cal_data['t2'] - cal_data['t1']) < 0.75
+    condition_dict['All (annual)'] = (cal_data['t2'] - cal_data['t1']) >= 0.75
+    condition_dict['All (seasonal)'] = (cal_data['t2'] - cal_data['t1']) < 0.75
     condition_dict['Geodetic'] = cal_data['obs_type'] == 'mb_geo'
     condition_dict['Glaciological'] = cal_data['obs_type'] == 'mb_glac'
     condition_dict['Glaciological (annual)'] = ((cal_data['seasonal/annual'] == 'annual') & 

analyze_simulation.py

@@ -333,24 +333,24 @@ def importGCMvarnearestneighbor_xarray(self, filename, vn, main_glac_rgi, dates_
 
 #%% Testing
 if __name__ == '__main__':
-#    gcm = GCM(name='CanESM2', rcp_scenario='rcp85')
-    gcm = GCM(name='ERA5')
-#    gcm = GCM(name='ERA-Interim')
-    
-    main_glac_rgi = modelsetup.selectglaciersrgitable(rgi_regionsO1=input.rgi_regionsO1, rgi_regionsO2 = 'all',
-                                                      rgi_glac_number=input.rgi_glac_number)
-    dates_table = modelsetup.datesmodelrun(startyear=1980, endyear=2017, spinupyears=0, 
-                                           option_wateryear=input.gcm_wateryear)
-
-    # Air temperature [degC], Precipitation [m], Elevation [masl], Lapse rate [K m-1]
-    gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table)
-    gcm_prec, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table)
-    gcm_elev = gcm.importGCMfxnearestneighbor_xarray(gcm.elev_fn, gcm.elev_vn, main_glac_rgi)
-    if gcm.name == 'ERA-Interim' or gcm.name == 'ERA5':
-        gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table)
-    if gcm.name == 'ERA5':
-        gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, 
-                                                                        dates_table)
+##    gcm = GCM(name='CanESM2', rcp_scenario='rcp85')
+#    gcm = GCM(name='ERA5')
+##    gcm = GCM(name='ERA-Interim')
+#    
+#    main_glac_rgi = modelsetup.selectglaciersrgitable(rgi_regionsO1=input.rgi_regionsO1, rgi_regionsO2 = 'all',
+#                                                      rgi_glac_number=input.rgi_glac_number)
+#    dates_table = modelsetup.datesmodelrun(startyear=1980, endyear=2017, spinupyears=0, 
+#                                           option_wateryear=input.gcm_wateryear)
+#
+#    # Air temperature [degC], Precipitation [m], Elevation [masl], Lapse rate [K m-1]
+#    gcm_temp, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.temp_fn, gcm.temp_vn, main_glac_rgi, dates_table)
+#    gcm_prec, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.prec_fn, gcm.prec_vn, main_glac_rgi, dates_table)
+#    gcm_elev = gcm.importGCMfxnearestneighbor_xarray(gcm.elev_fn, gcm.elev_vn, main_glac_rgi)
+#    if gcm.name == 'ERA-Interim' or gcm.name == 'ERA5':
+#        gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table)
+#    if gcm.name == 'ERA5':
+#        gcm_tempstd, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.tempstd_fn, gcm.tempstd_vn, main_glac_rgi, 
+#                                                                        dates_table)
 #    else:
 #        gcm_lr = np.tile(ref_lr_monthly_avg, int(gcm_temp.shape[1]/12))
 #    # COAWST data has two domains, so need to merge the two domains
@@ -368,4 +368,18 @@ def importGCMvarnearestneighbor_xarray(self, filename, vn, main_glac_rgi, dates_
 #                ~(input.coawst_d02_lon_min <= glac_lon <= input.coawst_d02_lon_max)):
 #                gcm_prec[glac,:] = gcm_prec_d01[glac,:]
 #                gcm_temp[glac,:] = gcm_temp_d01[glac,:]
-#                gcm_elev[glac] = gcm_elev_d01[glac]
+#                gcm_elev[glac] = gcm_elev_d01[glac]
+        
+    #%%
+#    # Get range of dates
+#    rcp_scenario = 'rcp85'
+#    gcm_names = ['bcc-csm1-1', 'CanESM2', 'CESM1-CAM5', 'CCSM4', 'CNRM-CM5', 'CSIRO-Mk3-6-0', 'FGOALS-g2', 'GFDL-CM3', 
+#             'GFDL-ESM2G', 'GFDL-ESM2M', 'GISS-E2-R', 'HadGEM2-ES', 'IPSL-CM5A-LR', 'IPSL-CM5A-MR', 'MIROC-ESM', 
+#             'MIROC-ESM-CHEM', 'MIROC5', 'MPI-ESM-LR', 'MPI-ESM-MR', 'MRI-CGCM3', 'NorESM1-M', 'NorESM1-ME']
+#    for gcm_name in gcm_names:
+#        print(gcm_name)
+#        ds = xr.open_dataset(input.cmip5_fp_var_prefix + rcp_scenario + input.cmip5_fp_var_ending + 
+#                             'tas' + '_mon_' + gcm_name + '_' + rcp_scenario + '_r1i1p1_native.nc')
+#        
+#        print('  ', ds.time[0].values,
+#              '\n  ', ds.time[-1].values)

class_climate.py

@@ -128,16 +128,16 @@ def glac_fromcsv(csv_fullfn, cn='RGIId'):
 #rgi_glac_number = glac_num_fromrange(1,5)
 #rgi_glac_number = get_same_glaciers(output_filepath + 'cal_opt1/reg1/')
 #rgi_glac_number = get_shean_glacier_nos(rgi_regionsO1[0], 1, option_random=1)
-glac_no = None
+#glac_no = None
 #glac_no = glac_fromcsv(main_directory + '/../qgis_himat/trishuli_and_naltar_RGIIds.csv')
-#glac_no = ['15.01152']
+glac_no = ['15.03473']
 if glac_no is not None:
     rgi_regionsO1 = sorted(list(set([int(x.split('.')[0]) for x in glac_no])))
 
 # ===== CLIMATE DATA ===== 
 # Reference period runs
-ref_gcm_name = 'ERA-Interim'    # reference climate dataset
-#ref_gcm_name = 'ERA5'           # reference climate dataset
+#ref_gcm_name = 'ERA-Interim'    # reference climate dataset
+ref_gcm_name = 'ERA5'           # reference climate dataset
 
 #startyear = 1980                # first year of model run (reference dataset)
 #endyear = 2018                  # last year of model run (reference dataset)
@@ -189,7 +189,7 @@ def glac_fromcsv(csv_fullfn, cn='RGIId'):
 
 #%% ===== CALIBRATION OPTIONS =====
 # Calibration option (1 = minimization, 2 = MCMC, 3=HH2015, 4=modified HH2015)
-option_calibration = 2
+option_calibration = 4
 # Calibration datasets ('shean', 'larsen', 'mcnabb', 'wgms_d', 'wgms_ee', 'group')
 cal_datasets = ['shean']
 #cal_datasets = ['shean']
@@ -286,8 +286,8 @@ def glac_fromcsv(csv_fullfn, cn='RGIId'):
 # Model parameter column names and filepaths
 modelparams_colnames = ['lrgcm', 'lrglac', 'precfactor', 'precgrad', 'ddfsnow', 'ddfice', 'tempsnow', 'tempchange']
 # Model parameter filepath
-#modelparams_fp = output_filepath + 'cal_opt' + str(option_calibration) + '/'
-modelparams_fp = output_filepath + 'cal_opt2_spc_20190806/'
+modelparams_fp = output_filepath + 'cal_opt' + str(option_calibration) + '/'
+#modelparams_fp = output_filepath + 'cal_opt2_spc_20190806/'
 
 #%% CLIMATE DATA
 # ERA-INTERIM (Reference data)
@@ -378,15 +378,24 @@ def glac_fromcsv(csv_fullfn, cn='RGIId'):
     option_shift_elevbins_20m = 0   # option to shift bins by 20 m (needed since off by 20 m, seem email 5/24/2018)
     # Dictionary of hypsometry filenames
     hyps_filepath = main_directory + '/../IceThickness_Farinotti/output/'
-    hyps_filedict = {1:  'area_km2_01_Farinotti2019_10m.csv',}
+    hyps_filedict = {1:  'area_km2_01_Farinotti2019_10m.csv',
+                     13: 'area_km2_13_Farinotti2019_10m.csv',
+                     14: 'area_km2_14_Farinotti2019_10m.csv',
+                     15: 'area_km2_15_Farinotti2019_10m.csv'}
     hyps_colsdrop = ['RGIId']
     # Thickness data
     thickness_filepath = main_directory + '/../IceThickness_Farinotti/output/'
-    thickness_filedict = {1:  'thickness_m_01_Farinotti2019_10m.csv'}
+    thickness_filedict = {1:  'thickness_m_01_Farinotti2019_10m.csv',
+                          13: 'thickness_m_13_Farinotti2019_10m.csv',
+                          14: 'thickness_m_14_Farinotti2019_10m.csv',
+                          15: 'thickness_m_15_Farinotti2019_10m.csv'}
     thickness_colsdrop = ['RGIId']
     # Width data
     width_filepath = main_directory + '/../IceThickness_Farinotti/output/'
-    width_filedict = {1:  'width_km_01_Farinotti2019_10m.csv'}
+    width_filedict = {1:  'width_km_01_Farinotti2019_10m.csv',
+                      13: 'width_km_13_Farinotti2019_10m.csv',
+                      14: 'width_km_14_Farinotti2019_10m.csv',
+                      15: 'width_km_15_Farinotti2019_10m.csv'}
     width_colsdrop = ['RGIId']
 
 elif hyps_data == 'Huss':
@@ -440,6 +449,12 @@ def glac_fromcsv(csv_fullfn, cn='RGIId'):
                     16: 'width_16_Huss_LowLatitudes_10m.csv',
                     17: 'width_17_Huss_SouthernAndes_10m.csv'}
     width_colsdrop = ['RGI-ID','Cont_range']
+    
+# Debris datasets
+debris_fp = main_directory + '/../IceThickness_Farinotti/output/'
+debris_filedict = {15:  'meltfactor_15_10m.csv'}
+debris_colsdrop = ['RGIId']
+
 
 #%% MODEL TIME FRAME DATA
 # Models require complete data for each year such that refreezing, scaling, etc. can be calculated

pygem_input.py

@@ -314,6 +314,12 @@ def runmassbalance(modelparameters, glacier_rgi_table, glacier_area_initial, ice
                     # DDF based on surface type [m w.e. degC-1 day-1]
                     for surfacetype_idx in surfacetype_ddf_dict: 
                         surfacetype_ddf[surfacetype == surfacetype_idx] = surfacetype_ddf_dict[surfacetype_idx]
+                    if input.option_surfacetype_debris == 1:
+                        print('\n\nLOAD THE MELTFACTOR DATASET over areas that are not firn\n\n')
+                        
+                        if year == 0 and month == 0:
+                            print('\nDELETE ME\n surfacetype_ddf[glac_idx]:', surfacetype_ddf[glac_idx_t0])
+                        
                     bin_meltglac[glac_idx_t0,step] = surfacetype_ddf[glac_idx_t0] * melt_energy_available[glac_idx_t0]
                     # TOTAL MELT (snow + glacier)
                     #  off-glacier need to include melt of refreeze because there are no glacier dynamics,

pygemfxns_massbalance.py

@@ -2,17 +2,18 @@
 pygemfxns_preprocessing.py is a list of the model functions that are used to preprocess the data into the proper format.
 
 """
+print('PROCESS FILE USING RASTERENV ENVIRONMENT')
 
 # Built-in libraries
 import os
 import argparse
 # External libraries
 from osgeo import gdal
-import geopandas as gpd
+#import geopandas as gpd
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
-import shapely
+#import shapely
 
 from pygeotools.lib import iolib, warplib, geolib, timelib, malib
 
@@ -45,39 +46,43 @@ def plot3panel(dem_list, clim=None, titles=None, cmap='inferno', label=None, ove
         fig.savefig(fn, bbox_inches='tight', pad_inches=0, dpi=150)
 
 #Input DEM filenames
-dem_ref_fn = '/Users/davidrounce/Documents/Dave_Rounce/HiMAT/DEMs/Alaska_albers_V3_mac/Alaska_albers_V3.tif'
+dem_ref_fn = None
+# dem_ref_fn = '/Users/davidrounce/Documents/Dave_Rounce/HiMAT/DEMs/Alaska_albers_V3_mac/Alaska_albers_V3.tif'
 thickness_fp_prefix = ('/Users/davidrounce/Documents/Dave_Rounce/HiMAT/IceThickness_Farinotti/' +
                        'composite_thickness_RGI60-all_regions/')
-dem_farinotti_fp = ('/Users/davidrounce/Documents/Dave_Rounce/HiMAT/IceThickness_Farinotti/surface_DEMs_RGI60/' +
-                    'surface_DEMs_RGI60-01/')
+# dem_farinotti_fp = ('/Users/davidrounce/Documents/Dave_Rounce/HiMAT/IceThickness_Farinotti/surface_DEMs_RGI60/' +
+#                     'surface_DEMs_RGI60-01/')
+dem_farinotti_fp_prefix = '/Users/davidrounce/Documents/Dave_Rounce/HiMAT/IceThickness_Farinotti/surface_DEMs_RGI60/'
 output_fp = '/Users/davidrounce/Documents/Dave_Rounce/HiMAT/IceThickness_Farinotti/output/'
 fig_fp = output_fp + 'figures/'
 if os.path.exists(output_fp) == False:
     os.makedirs(output_fp)
 if os.path.exists(fig_fp) == False:
     os.makedirs(fig_fp)
 
-rgi_regionsO1 = [1]                 # RGI Order 1 regions
+rgi_regionsO1 = [14]                 # RGI Order 1 regions
 binsize = 10                        # elevation bin (must be an integer greater than 1)
+dem_poorquality_switch = True       # Switch to filter poor quality DEMs if another DEM is available
 dem_poorquality_threshold = 200     # threshold used to identify problems with Farinotti DEM
 option_plot_DEMsraw = True          # Option to plot the raw DEMs
 option_plot_DEMs = False             # Option to plot the masked DEMs
 debug = False
 
-# ===== LOAD GLACIERS ======
+# ======
 glacno_wpoor_DEM = []
 for region in rgi_regionsO1:
 
     thickness_fp = thickness_fp_prefix + 'RGI60-' + str(region).zfill(2) + '/'
+    dem_farinotti_fp = dem_farinotti_fp_prefix + 'surface_DEMs_RGI60-' + str(region).zfill(2) + '/'
 
     glacno_list = []
     for i in os.listdir(thickness_fp):
         if i.endswith('_thickness.tif'):
             glacno_list.append(i.split('-')[1].split('_')[0])
     glacno_list = sorted(glacno_list)
 
-    # print('\n\nDELETE ME - SWITCH TO COMPLETE LIST\n\n')
-    # glacno_list = ['01.03622']
+#    print('\n\nDELETE ME - SWITCH TO COMPLETE LIST\n\n')
+#    glacno_list = ['15.02228']
     # glacno_list = glacno_list[10000:10010]
 
     # Load RGI glacier data
@@ -97,32 +102,39 @@ def plot3panel(dem_list, clim=None, titles=None, cmap='inferno', label=None, ove
     main_glac_width['RGIId'] = main_glac_rgi.RGIId.values
     main_glac_length['RGIId'] = main_glac_rgi.RGIId.values
     main_glac_slope['RGIId'] = main_glac_rgi.RGIId.values
-
+        
     # ===== PROCESS EACH GLACIER ======
     for nglac, glacno in enumerate(glacno_list):
         # print(nglac, glacno)
         thickness_fn = thickness_fp + 'RGI60-' + glacno + '_thickness.tif'
         dem_farinotti_fn = dem_farinotti_fp + 'surface_DEM_RGI60-' + glacno + '.tif'
-
+        
         # Reproject, resample, warp rasters to common extent, grid size, etc.
         #  note: use thickness for the reference to avoid unrealistic extrapolations, e.g., negative thicknesses
         #        also using equal area increases areas significantly compared to RGI
-        raster_fn_list = [dem_ref_fn, dem_farinotti_fn, thickness_fn]
+        raster_fn_list = [dem_farinotti_fn, thickness_fn]
+        if dem_ref_fn is not None:
+            raster_fn_list.append(dem_ref_fn)
+        
+        print(raster_fn_list)
+        
         ds_list = warplib.memwarp_multi_fn(raster_fn_list, extent='intersection', res='min', t_srs=thickness_fn)
-        # print('\n\nSWITCH BACK TO THICKNESS_FN AFTER OTHERS CORRECTED!\n\n')
-        # ds_list = warplib.memwarp_multi_fn(raster_fn_list, extent='intersection', res='min', t_srs=dem_ref_fn)
 
         # masked arrays using ice thickness estimates
-        dem_ref_raw, dem_far_raw, thickness = [iolib.ds_getma(i) for i in ds_list]
-        dem_ref = dem_ref_raw.copy()
-        dem_ref.mask = thickness.mask
+        if dem_ref_fn is not None:
+            dem_ref_raw, dem_far_raw, thickness = [iolib.ds_getma(i) for i in ds_list]
+            dem_ref = dem_ref_raw.copy()
+            dem_ref.mask = thickness.mask
+        else:
+            dem_far_raw, thickness = [iolib.ds_getma(i) for i in ds_list]
         dem_far = dem_far_raw.copy()
         dem_far.mask = thickness.mask
-
+        
         # DEM selection for binning computations
         # if exceeds threshold, then use the reference
-        if (abs(main_glac_rgi.loc[nglac,'Zmin'] - dem_far.min()) > dem_poorquality_threshold or
-            abs(main_glac_rgi.loc[nglac,'Zmax'] - dem_far.max()) > dem_poorquality_threshold):
+        if ((abs(main_glac_rgi.loc[nglac,'Zmin'] - dem_far.min()) > dem_poorquality_threshold or
+             abs(main_glac_rgi.loc[nglac,'Zmax'] - dem_far.max()) > dem_poorquality_threshold) 
+             and dem_ref_fn is not None):
             print('  Check Glacier ' + glacno + ': use Christian DEM instead of Farinotti')
             print('\n     RGI Zmin/Zmax:', main_glac_rgi.loc[nglac,'Zmin'], '/', main_glac_rgi.loc[nglac,'Zmax'])
             print('     Farinotti Zmin/Zmax:', np.round(dem_far.min(),0), '/', np.round(dem_far.max(),0))
@@ -148,7 +160,7 @@ def plot3panel(dem_list, clim=None, titles=None, cmap='inferno', label=None, ove
         else:
             dem = dem_far
             dem_raw = dem_far_raw
-
+            
         #Extract x and y pixel resolution (m) from geotransform
         gt = ds_list[0].GetGeoTransform()
         px_res = (gt[1], -gt[5])
@@ -212,7 +224,7 @@ def plot3panel(dem_list, clim=None, titles=None, cmap='inferno', label=None, ove
 
         elev_bin_edges = np.arange(elev_bin_min, elev_bin_max+binsize, binsize)
         elev_bins = (elev_bin_edges[0:-1] + binsize/2).astype(int)
-
+        
         # Hypsometry [km2]
         #  must used .compressed() in histogram to exclude masked values
         hist, elev_bin_edges = np.histogram(dem.reshape(-1).compressed(), bins=elev_bin_edges)
@@ -243,6 +255,13 @@ def plot3panel(dem_list, clim=None, titles=None, cmap='inferno', label=None, ove
 
         # Width [km] - based on length (inherently slope) and bin area
         bin_width = bin_hyps / bin_length
+        
+        # Remove negative values
+        bin_hyps[bin_hyps < 0] = 0
+        bin_thickness[bin_thickness < 0] = 0
+        bin_width[bin_width < 0] = 0
+        bin_length[bin_length < 0] = 0
+        bin_slope[bin_slope < 0] = 0
 
         # Record properties
         # Check if need to expand columns
@@ -267,12 +286,12 @@ def plot3panel(dem_list, clim=None, titles=None, cmap='inferno', label=None, ove
     main_glac_width = main_glac_width.fillna(0)
     main_glac_length = main_glac_length.fillna(0)
     main_glac_slope = main_glac_slope.fillna(0)
-    # Remove negative values
-    main_glac_hyps[main_glac_hyps < 0] = 0
-    main_glac_thickness[main_glac_thickness < 0] = 0
-    main_glac_width[main_glac_width < 0] = 0
-    main_glac_length[main_glac_length < 0] = 0
-    main_glac_slope[main_glac_slope < 0] = 0
+#     # Remove negative values
+#     main_glac_hyps[main_glac_hyps < 0] = 0
+#     main_glac_thickness[main_glac_thickness < 0] = 0
+#     main_glac_width[main_glac_width < 0] = 0
+#     main_glac_length[main_glac_length < 0] = 0
+#     main_glac_slope[main_glac_slope < 0] = 0
     # Export results
     main_glac_hyps.to_csv(output_fp + 'area_km2_' + "{:02d}".format(region) + '_Farinotti2019_' +
                           str(binsize) + 'm.csv', index=False)
@@ -284,12 +303,3 @@ def plot3panel(dem_list, clim=None, titles=None, cmap='inferno', label=None, ove
                             str(binsize) + 'm.csv', index=False)
     main_glac_slope.to_csv(output_fp + 'slope_deg_' + "{:02d}".format(region) + '_Farinotti2019_' +
                            str(binsize) + 'm.csv', index=False)
-
-##%%
-#import pandas as pd
-#import pygem_input as input
-#area = pd.read_csv(input.hyps_filepath + 'area_km2_01_Farinotti2019_10m_old.csv')
-#thickness = pd.read_csv(input.hyps_filepath + 'thickness_m_01_Farinotti2019_10m_old.csv')
-#length = pd.read_csv(input.hyps_filepath + 'length_km_01_Farinotti2019_10m_old.csv')
-#slope = pd.read_csv(input.hyps_filepath + 'slope_deg_01_Farinotti2019_10m_old.csv')
-#width = pd.read_csv(input.hyps_filepath + 'width_km_01_Farinotti2019_10m_old.csv')