<jupyter_start><jupyter_text>Walkthrough: Simulating Greenland glaciers with SERMeQ Welcome!  This notebook will walk you through the necessary steps to simulate one or more Greenland outlet glaciers using SERMeQ.  We will use geographic data and identifiers from [MEaSUREs](https://nsidc.org/data/NSIDC-0642/versions/1), bed topography from [BedMachine v3](https://nsidc.org/data/idbmg4), and surface mass balance from [HIRHAM/ERA-Interim](http://prudence.dmi.dk/data/temp/RUM/HIRHAM/GREENLAND/).  Importing necessary Python packages<jupyter_code>from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
import pyproj as pyproj
import csv
import collections
import datetime
from matplotlib import cm
from scipy import interpolate
from scipy.ndimage import gaussian_filter
## Special import for SERMeQ modules
import sys
sys.path.insert(0, '/Users/lizz/Documents/GitHub/plastic-networks')
from SERMeQ.plastic_utilities_v2 import *
from SERMeQ.GL_model_tools import *
from SERMeQ.flowline_class_hierarchy import *<jupyter_output><empty_output><jupyter_text>Importing input data<jupyter_code>print('Reading in surface topography')
gl_bed_path ='/Users/lizz/Documents/GitHub/Data_unsynced/BedMachine-Greenland/BedMachineGreenland-2017-09-20.nc'
fh = Dataset(gl_bed_path, mode='r')
xx = fh.variables['x'][:].copy() #x-coord (polar stereo (70, 45))
yy = fh.variables['y'][:].copy() #y-coord
s_raw = fh.variables['surface'][:].copy() #surface elevation
h_raw=fh.variables['thickness'][:].copy() # Gridded thickness
b_raw = fh.variables['bed'][:].copy() # bed topo
thick_mask = fh.variables['mask'][:].copy()
#mask values: 0=ocean, 1=ice-free land, 2=grounded ice, 3=floating ice, 4=non-Greenland land
ss = np.ma.masked_where(thick_mask !=2, s_raw) #select only grounded ice cells
hh = np.ma.masked_where(thick_mask !=2, h_raw) 
bb = np.ma.masked_where(thick_mask !=2, b_raw)
## Down-sampling
X = xx[::2]
Y = yy[::2]
S = ss[::2, ::2]
H = hh[::2, ::2] 
B = bb[::2, ::2]
M = thick_mask[::2,::2]
fh.close()

#Smoothing bed
unsmoothB = B
smoothB = gaussian_filter(B, 2)
smoothS = gaussian_filter(S, 2) #17 Jan 19 - smoothing S as well for consistency with auto-selected networks
#B_processed = np.ma.masked_where(thick_mask !=2, smoothB)

S_interp = interpolate.RectBivariateSpline(X, Y[::-1], smoothS.T[::, ::-1])
H_interp = interpolate.RectBivariateSpline(X, Y[::-1], H.T[::, ::-1])
B_interp = interpolate.RectBivariateSpline(X, Y[::-1], smoothB.T[::, ::-1])
def NearestMaskVal(x0,y0):
    x_idx = np.abs(X-x0).argmin()
    y_idx = np.abs(Y-y0).argmin()
    return M[y_idx, x_idx]
print('Reading in surface mass balance from 1981-2010 climatology')
gl_smb_path ='/Users/lizz/Documents/GitHub/Data_unsynced/HIRHAM5-SMB/DMI-HIRHAM5_GL2_ERAI_1980_2014_SMB_YM.nc'
fh2 = Dataset(gl_smb_path, mode='r')
x_lon = fh2.variables['lon'][:].copy() #x-coord (latlon)
y_lat = fh2.variables['lat'][:].copy() #y-coord (latlon)
#zs = fh2.variables['height'][:].copy() #height in m - is this surface elevation or SMB?
ts = fh2.variables['time'][:].copy()
smb_raw = fh2.variables['smb'][:].copy()
fh2.close()<jupyter_output><empty_output><jupyter_text>The bed topography from BedMachine is on a Polar Stereographic North grid (EPSG code 3413).  The surface mass balance from HIRHAM is on a lat-lon grid with WGS 84 datum (EPSG code 4326).  We use `pyproj` to align the grids, then interpolate the yearly mass balance so that we can extract at each glacier site.<jupyter_code>print('Now transforming coordinate system of SMB')
wgs84 = pyproj.Proj("+init=EPSG:4326") # LatLon with WGS84 datum used by GPS units and Google Earth
psn_gl = pyproj.Proj("+init=epsg:3413") # Polar Stereographic North used by BedMachine (as stated in NetCDF header)
xs, ys = pyproj.transform(wgs84, psn_gl, x_lon, y_lat)
Xmat, Ymat = np.meshgrid(X, Y)
## Hindcasting SMB: year-specific 2006-2014
SMB_dict = {} #set up a dictionary of surface mass balance fields indexed by year
for year in range(2006, 2015):
    index = year - 2015 #so that 2014 will be smb_raw[-1], etc.
    smb_year = smb_raw[index][0]
    regridded_smb_year = interpolate.griddata((xs.ravel(), ys.ravel()), smb_year.ravel(), (Xmat, Ymat), method='nearest')
    SMB_dict[year] = interpolate.interp2d(X, Y, regridded_smb_year, kind='linear')<jupyter_output><empty_output><jupyter_text>Loading in saved glaciers We identify flowline coordinates and optimise the yield strength for each glacier in a separate pre-processing routine.  We load in the saved results, tagged by their MEaSUREs glacier IDs.<jupyter_code>print('Reading in optimal yield strength dictionary')
optimal_taus_fpath = '/Users/lizz/Documents/GitHub/Data_unsynced/Auto_selected-networks/Optimization_analysis/bestfit_taus-B_S_smoothing-fromdate_2019-01-17.csv'
f_ot = open(optimal_taus_fpath, 'r')
header = f_ot.readline()
hdr = header.strip('\r\n')
optimal_taus = {}
lines = f_ot.readlines()
for i, l in enumerate(lines):
    linstrip = l.strip('\r\n')
    parts = linstrip.split(',')
    gid = int(parts[0]) #MEaSUREs glacier ID is first
    tau_y = float(parts[1]) #yield strength in Pa
    yieldtype = parts[2] #'constant' or 'variable' string
    optimal_taus[gid] = (tau_y, yieldtype)<jupyter_output><empty_output><jupyter_text>We identify the base filepath where glacier networks are stored, and which glacier IDs should be available to simulate.  Choose which glacier(s) to simulate by changing `gids_totest` in the cell below.<jupyter_code>base_fpath = '/Users/lizz/Documents/GitHub/Data_unsynced/Auto_selected-networks/Gld-autonetwork-GID'

glacier_ids = range(1,195) #range of MEaSUREs IDs
not_present = (93, 94, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 169) #glacier IDs missing from set
added_jan19 = (139, 140, 141, 142, 143, 159, 161, 172, 173, 177)
for n in not_present:
    try:
        glacier_ids.remove(n)
    except ValueError:
        pass
gids_totest = (3,) #test specific glaciers; set ==glacier_ids to test all<jupyter_output><empty_output><jupyter_text>Simulating glacier networks Choose the settings for the simulation.  If you had correct input data, you could change `testyears` to simulate a longer period (or change the last argument to simulate denser/coarser in time).  You can also switch debug mode on or off by changing `db`, change Glen's A parameter by adjusting `test_A`, and, if applicable, change the climate scenario by changing `persistence`.<jupyter_code>## Simulation settings
testyears = arange(0, 9, 0.25) #test only 2006-2015, for comparison
start_year=2006 #determined by which MEaSUREs termini we used to initialize a given set
db = True #debug mode, True for frequent print statements
test_A, icetemp = 3.5E-25, 'min10C' # -10 C, good guess for Greenland
scenario = 'persistence'
output_heavy = False #pref for output--should be False for efficient running, True for saving full surface profiles<jupyter_output><empty_output><jupyter_text>Now, run the simulation for each glacier selected and save the output as a pickle.<jupyter_code>bad_gids = []

for gid in gids_totest:
    print 'Reading in glacier ID: '+str(gid)
    if gid in added_jan19:
        filename = base_fpath+str(gid)+'-date_2019-01-10.csv'
    elif gid<160:
        filename = base_fpath+str(gid)+'-date_2018-10-03.csv'
    else:
        filename = base_fpath+str(gid)+'-date_2018-10-04.csv' #workaround because I ran these in batches and saved them with the date
    
    coords_list = Flowline_CSV(filename, has_width=True, flip_order=False)
    nlines = len(coords_list)
    branch_0 = Branch(coords=coords_list[0], index=0, order=0) #saving central branch as main
    branch_list = [branch_0]
    if nlines>0:
        for l in range(1, nlines):
            branch_l = Branch(coords = coords_list[l], index=l, order=1, flows_to=0)
            branch_list.append(branch_l)
    nw = PlasticNetwork(name='GID'+str(gid), init_type='Branch', branches=branch_list, main_terminus=branch_0.coords[0])
    nw.make_full_lines()

    print 'Now processing glacier ID: '+str(gid)
    nw.process_full_lines(B_interp, S_interp, H_interp)
    nw.remove_floating()
    nw.make_full_lines()
    nw.process_full_lines(B_interp, S_interp, H_interp)
    nw.network_tau = optimal_taus[gid][0]
    nw.network_yield_type = optimal_taus[gid][1]
    nw.network_ref_profiles()
    
    ## Simulations forced by SMB
    time_varying_smb = [[0.001*(1000/nw.rho_ice)*SMB_dict[int(floor(yr))](nw.flowlines[0].coords[i,0], nw.flowlines[0].coords[i,1]) for i in range(len(nw.flowlines[0].coords))] for yr in start_year+testyears]
    catchment_smb_vals = [np.mean(time_varying_smb[i]) for i in range(len(time_varying_smb))] #because grid spacing is very fine, mean of all points should approximate integral/length
    print 'shape(catchment_smb_vals) = {}'.format(shape(catchment_smb_vals))
    nw.smb_alphadot = catchment_smb_vals[0] #initial SMB
    print 'a_dot from SMB: {}'.format(nw.smb_alphadot)
    nw.terminus_adot = time_varying_smb[0][0]
    print 'Terminus a_dot: {}'.format(nw.terminus_adot)
    try:
        nw.terminus_time_evolve(testyears=testyears, alpha_dot_variable=catchment_smb_vals, 
                                dL=1/L0, separation_buffer=10000/L0, rate_factor=test_A,
                                has_smb=True, terminus_balance=nw.terminus_adot, submarine_melt = 0, 
                                debug_mode=db, output_heavy=output_heavy, sl_contribution=False
                                )
    except:
        bad_gids.append(gid)   
        continue 
    
    print 'Saving output for {}'.format(nw.name)
    fn = str(nw.name)
    fn1 = fn.replace(" ", "")
    fn2 = fn1.replace("[", "-")
    fn3 = fn2.replace("/", "_")
    fn4 = fn3.replace("]", "")
    fn5 = '/Users/lizz/Desktop/Hindcasted_networks/'+fn4+'-{}-{}-{}ice-{}a_dt025a.pickle'.format(datetime.date.today(), scenario, icetemp, int(max(testyears)))
    nw.save_network(filename=fn5)<jupyter_output><empty_output><jupyter_text>Finally, report which glaciers failed to run in a csv output of bad glacier IDs.<jupyter_code>##Output errors to csv file
error_fn = '/Users/lizz/Desktop/Hindcasted_networks/error_gids-{}.csv'.format(datetime.date.today())
np.savetxt(error_fn, np.asarray(bad_gids), delimiter=',')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1254392815306521, 'fraction_numerical': 0.03162827020695041, 'mean_word_length': 4.797962648556876, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 34, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '537449', 'license': 'MIT', 'n_tokens_mistral': 4019, 'n_tokens_neox': 3616, 'n_words': 915, 'path': 'Hindcast-MS/Gld_sim-walkthrough.ipynb', 'repo_name': 'srcaddle/plastic-networks', 'size': '10244'}