Initial commit

GreenlandSG.py

@@ -0,0 +1,174 @@
+#%%
+
+import rasterio
+import netCDF4 as nc
+import numpy as np
+from scipy.interpolate import RegularGridInterpolator
+import gc
+import pickle
+from SavitzkyGolay import sgolay2d
+
+fp = r'./greenland_vel_mosaic250_vx_v1.tif'
+img = rasterio.open(fp)
+vx = img.read(1)
+
+fp = r'./greenland_vel_mosaic250_vy_v1.tif'
+img = rasterio.open(fp)
+vy = img.read(1)
+
+
+y = np.array(img.xy(np.arange(vy.shape[0]),np.zeros(vy.shape[0]))[1])
+x = np.array(img.xy(np.zeros(vy.shape[1]),np.arange(vy.shape[1]))[0])
+
+
+
+
+##Add nans
+vx[vx<-1e9]=np.nan
+vy[vy<-1e9]=np.nan
+
+#flip on axis 0 to makes postitive
+vx = np.flip(vx,0)
+vy = np.flip(vy,0)
+y = np.flip(y,0)
+
+#Trim
+# xmin = 100000
+# xmax = 295000
+# ymin = -2260000
+# ymax = -1040000
+xmin = 110000
+xmax = 310000
+ymin = -1940000
+ymax = -1500000
+
+
+xminind = np.abs(x-xmin).argmin()
+xmaxind = np.abs(x-xmax).argmin()
+yminind = np.abs(y-ymin).argmin()
+ymaxind = np.abs(y-ymax).argmin()
+
+x = x[xminind:xmaxind]
+y = y[yminind:ymaxind]
+vx = vx[yminind:ymaxind,xminind:xmaxind]
+vy = vy[yminind:ymaxind,xminind:xmaxind]
+
+
+#Savitzky-Golay filtering
+window=51
+order =5
+vx_filt = sgolay2d(vx,window,order)
+vy_filt = sgolay2d(vy,window,order)
+
+dudx,dudy = sgolay2d(vx,window,order,'both')/(x[1]-x[0])
+dvdx,dvdy = sgolay2d(vy,window,order,'both')/(y[1]-y[0])
+
+
+fn = 'BedMachineGreenland-v5.nc'
+ds = nc.Dataset(fn)
+
+
+surface=np.squeeze(ds['surface'][:])
+thickness=np.squeeze(ds['thickness'][:])
+bed=np.squeeze(ds['bed'][:])
+y2 = np.squeeze(ds['y'][:])
+x2 = np.squeeze(ds['x'][:])
+
+y2=np.flip(y2,0)
+surface = np.flip(surface,0)
+thickness = np.flip(thickness,0)
+bed = np.flip(bed,0)
+
+window=51
+order =5
+dbeddx,dbeddy = sgolay2d(bed,window,order,'both')/(x[1]-x[0])
+dsurfdx,dsurfdy = sgolay2d(surface,window,order,'both')/(y[1]-y[0])
+
+
+fn = 'Greenland1km.nc'
+ds = nc.Dataset(fn)
+acc = np.squeeze(ds['presprcp'][:])
+temp = np.squeeze(ds['presartm'][:])
+y3 = np.squeeze(ds['y'][:])
+x3 = np.squeeze(ds['x'][:])
+
+
+
+
+
+
+class gl:
+    def __init__(self):
+        self.vx_interp = RegularGridInterpolator((x,y), vx.T)
+        self.vy_interp = RegularGridInterpolator((x,y), vy.T)
+        self.dudx_interp = RegularGridInterpolator((x,y), dudx.T)
+        self.dudy_interp = RegularGridInterpolator((x,y), dudy.T)
+        self.dvdx_interp = RegularGridInterpolator((x,y), dvdx.T)
+        self.dvdy_interp = RegularGridInterpolator((x,y), dvdy.T)
+        self.surf_interp = RegularGridInterpolator((x2,y2), surface.T)
+        self.thick_interp = RegularGridInterpolator((x2,y2), thickness.T)
+        self.bed_interp = RegularGridInterpolator((x2,y2), bed.T)
+        self.accumulation_interp = RegularGridInterpolator((x3,y3), acc.T)
+        self.surftemp_interp = RegularGridInterpolator((x3,y3), temp.T)
+        self.dbeddx_interp = RegularGridInterpolator((x2,y2),dbeddx.T)
+        self.dbeddy_interp = RegularGridInterpolator((x2,y2),dbeddy.T)
+        self.dsurfdx_interp = RegularGridInterpolator((x2,y2), dsurfdx.T)
+        self.dsurfdy_interp = RegularGridInterpolator((x2,y2), dsurfdy.T)
+
+        self.xmin = np.max([x.min(),x2.min(),x3.min()])
+        self.xmax = np.min([x.max(),x2.max(),x3.max()])
+        self.ymin = np.max([y.min(),y2.min(),y3.min()])
+        self.ymax = np.min([y.max(),y2.max(),y3.max()])
+
+
+    def interps(self,xi,yi):
+        self.xi=xi
+        self.yi=yi
+
+        if xi.ndim==2:
+            pts = np.reshape([xi,yi],(2,-1)).T
+        elif xi.ndim==1:
+            pts = np.stack((xi,yi),1)
+        else:
+            pts=[xi,yi]
+
+        self.vx=self.vx_interp(pts)
+        self.vy=self.vy_interp(pts)
+        self.dudx = self.dudx_interp(pts)
+        self.dudy = self.dudy_interp(pts)
+        self.dvdx = self.dvdx_interp(pts)
+        self.dvdy = self.dvdy_interp(pts)
+        self.surface = self.surf_interp(pts)
+        self.thickness = self.thick_interp(pts)
+        self.bed = self.bed_interp(pts)
+        self.accumulation = self.accumulation_interp(pts)
+        self.surftemp = self.surftemp_interp(pts)
+
+        self.surftemp = self.surftemp_interp(pts)
+        self.dbeddx = self.dbeddx_interp(pts)
+        self.dbeddy = self.dbeddy_interp(pts)
+        self.dsurfdx = self.dsurfdx_interp(pts)
+        self.dsurfdy = self.dsurfdy_interp(pts)
+        if xi.ndim==2:
+            self.vx=self.vx.reshape(xi.shape)
+            self.vy=self.vy.reshape(xi.shape)
+            self.dudx=self.dudx.reshape(xi.shape)
+            self.dudy=self.dudy.reshape(xi.shape)
+            self.dvdx=self.dvdx.reshape(xi.shape)
+            self.dvdy=self.dvdy.reshape(xi.shape)
+            self.surface=self.surface.reshape(xi.shape)
+            self.thickness=self.thickness.reshape(xi.shape)
+            self.bed=self.bed.reshape(xi.shape)
+            self.accumulation=self.accumulation.reshape(xi.shape)
+            self.surftemp=self.surftemp.reshape(xi.shape)
+            self.dbeddx=self.dbeddx.reshape(xi.shape)
+            self.dbeddy=self.dbeddy.reshape(xi.shape)
+            self.dsurfdx=self.dsurfdx.reshape(xi.shape)
+            self.dsurfdy=self.dsurfdy.reshape(xi.shape)
+
+        self.vmag=np.sqrt(self.vx**2+self.vy**2)
+
+
+
+
+# %%

Save2DpathGerber.py

@@ -0,0 +1,22 @@
+#%%
+
+import numpy as np
+import track
+import pickle
+import GreenlandSG as Greenland
+
+dt=50
+
+xc = 244691
+yc = -1544921
+
+data = Greenland.gl()
+
+p=track.path2d(100000,dt,xc,yc,data)
+
+filename = 'path2dSGdt'+str(dt)+'.pkl'
+
+with open(filename, 'wb') as f:
+   pickle.dump(p,f)
+
+# %%

SavitzkyGolay.py

@@ -0,0 +1,82 @@
+# from https://scipy.github.io/old-wiki/pages/Cookbook/SavitzkyGolay
+import numpy as np
+import scipy.signal
+
+def sgolay2d ( z, window_size, order, derivative=None):
+    """
+    """
+    # number of terms in the polynomial expression
+    n_terms = ( order + 1 ) * ( order + 2)  / 2.0
+
+    if  window_size % 2 == 0:
+        raise ValueError('window_size must be odd')
+
+    if window_size**2 < n_terms:
+        raise ValueError('order is too high for the window size')
+
+    half_size = window_size // 2
+
+    # exponents of the polynomial. 
+    # p(x,y) = a0 + a1*x + a2*y + a3*x^2 + a4*y^2 + a5*x*y + ... 
+    # this line gives a list of two item tuple. Each tuple contains 
+    # the exponents of the k-th term. First element of tuple is for x
+    # second element for y.
+    # Ex. exps = [(0,0), (1,0), (0,1), (2,0), (1,1), (0,2), ...]
+    exps = [ (k-n, n) for k in range(order+1) for n in range(k+1) ]
+
+    # coordinates of points
+    ind = np.arange(-half_size, half_size+1, dtype=np.float64)
+    dx = np.repeat( ind, window_size )
+    dy = np.tile( ind, [window_size, 1]).reshape(window_size**2, )
+
+    # build matrix of system of equation
+    A = np.empty( (window_size**2, len(exps)) )
+    for i, exp in enumerate( exps ):
+        A[:,i] = (dx**exp[0]) * (dy**exp[1])
+
+    # pad input array with appropriate values at the four borders
+    new_shape = z.shape[0] + 2*half_size, z.shape[1] + 2*half_size
+    Z = np.zeros( (new_shape) )
+    # top band
+    band = z[0, :]
+    Z[:half_size, half_size:-half_size] =  band -  np.abs( np.flipud( z[1:half_size+1, :] ) - band )
+    # bottom band
+    band = z[-1, :]
+    Z[-half_size:, half_size:-half_size] = band  + np.abs( np.flipud( z[-half_size-1:-1, :] )  -band ) 
+    # left band
+    band = np.tile( z[:,0].reshape(-1,1), [1,half_size])
+    Z[half_size:-half_size, :half_size] = band - np.abs( np.fliplr( z[:, 1:half_size+1] ) - band )
+    # right band
+    band = np.tile( z[:,-1].reshape(-1,1), [1,half_size] )
+    Z[half_size:-half_size, -half_size:] =  band + np.abs( np.fliplr( z[:, -half_size-1:-1] ) - band )
+    # central band
+    Z[half_size:-half_size, half_size:-half_size] = z
+
+    # top left corner
+    band = z[0,0]
+    Z[:half_size,:half_size] = band - np.abs( np.flipud(np.fliplr(z[1:half_size+1,1:half_size+1]) ) - band )
+    # bottom right corner
+    band = z[-1,-1]
+    Z[-half_size:,-half_size:] = band + np.abs( np.flipud(np.fliplr(z[-half_size-1:-1,-half_size-1:-1]) ) - band ) 
+
+    # top right corner
+    band = Z[half_size,-half_size:]
+    Z[:half_size,-half_size:] = band - np.abs( np.flipud(Z[half_size+1:2*half_size+1,-half_size:]) - band ) 
+    # bottom left corner
+    band = Z[-half_size:,half_size].reshape(-1,1)
+    Z[-half_size:,:half_size] = band - np.abs( np.fliplr(Z[-half_size:, half_size+1:2*half_size+1]) - band ) 
+
+    # solve system and convolve
+    if derivative == None:
+        m = np.linalg.pinv(A)[0].reshape((window_size, -1))
+        return scipy.signal.fftconvolve(Z, m, mode='valid')
+    elif derivative == 'col':
+        c = np.linalg.pinv(A)[1].reshape((window_size, -1))
+        return scipy.signal.fftconvolve(Z, -c, mode='valid')        
+    elif derivative == 'row':
+        r = np.linalg.pinv(A)[2].reshape((window_size, -1))
+        return scipy.signal.fftconvolve(Z, -r, mode='valid')        
+    elif derivative == 'both':
+        c = np.linalg.pinv(A)[1].reshape((window_size, -1))
+        r = np.linalg.pinv(A)[2].reshape((window_size, -1))
+        return scipy.signal.fftconvolve(Z, -r, mode='valid'), scipy.signal.fftconvolve(Z, -c, mode='valid')   

agedepth.py

@@ -0,0 +1,37 @@
+#%%
+
+import numpy as np
+from scipy.interpolate import interp1d,griddata
+import pandas as pd
+
+
+
+data = pd.read_excel("./SupplementaryMaterial/UpstreamCorrection_data_10Jun2021.xls",
+                     sheet_name=4)
+
+
+
+
+
+depth = data['depth'].to_numpy()[2:].astype(float)
+age = data['age'].to_numpy()[2:].astype(float)
+dist = data['upstream distance'].to_numpy()[2:].astype(float)*1000
+
+
+def time2depth(t):
+    return np.interp(t,age,depth)
+
+
+
+def dist2depth(s):
+    return np.interp(s,dist,depth)
+
+def depth2time(d):
+    return np.interp(d,depth,age)
+
+def depth2dist(d):
+    return np.interp(d,depth,dist)
+
+
+def depth2time(d):
+    return np.interp(d,depth,age)

densityfromdepth.py

@@ -0,0 +1,34 @@
+import numpy as np
+from scipy.interpolate import interp1d
+
+#Digitised from Fig. 9 "Initial results from geophysical surveys and shallow coring of the
+# Northeast Greenland Ice Stream (NEGIS)" by Vallelonga et al. (2014)"
+rawdata = np.array([[0.3383916990920882, -0.35901626644076146],\
+[0.4369649805447471, 5.930440910760732],\
+[0.5783398184176395, 16.985123041666967],\
+[0.6846952010376135, 30.20859729750667],\
+[0.8040207522697794, 49.28549115454492],\
+[0.8695201037613488, 63.75470137373827],\
+[0.9019455252918288, 76.07611164841097],\
+[0.9136186770428015, 86.5542085573124],\
+[0.9149156939040207, 101.19735398564033]])
+
+raw_d = rawdata[:,1]
+raw_rho = rawdata[:,0]
+raw_rho = raw_rho/raw_rho[-1]
+
+
+
+
+def densityfromdepth(depth):
+    # depth in m
+    # normalised density where 1 = ice sheet density
+    # from Fig. 9 "Initial results from geophysical surveys and shallow coring of the
+    # Northeast Greenland Ice Stream (NEGIS)" by Vallelonga et al. (2014)
+    # Digitised using https://apps.automeris.io/wpd/
+    #
+    f = interp1d(raw_d,raw_rho,kind='cubic',fill_value=1,bounds_error=False)
+    density = f(depth)
+    density[density>1]=1
+    return density
+