Adding changes from HPC HolschuhLab

Author: Nicholas Holschuh

find_nearest_xy.py

@@ -34,15 +34,15 @@ def find_nearest_xy(vector_2_search,value,how_many=1):
 
         comp_vec = vector_2_search[:,0] + vector_2_search[:,1]*np.sqrt(np.array(-1,dtype=complex));
         ss = value.shape;
-        if np.max(ss) > 1 & ss[0] > ss[1]:
+        if np.all([np.min(ss) > 1,ss[0] >= ss[1]]):
             value2 = value.transpose();
         elif ss[0] == 1:
             value2 = value.transpose();
         else:
             value2 = value;
 
 
-        value_search = value2[:,0] + value2[:,1]*np.sqrt(np.array(-1,dtype=complex));
+        value_search = value2[0,:] + value2[1,:]*np.sqrt(np.array(-1,dtype=complex));
 
         complex_dists = np.tile(comp_vec,(len(value_search),1)).transpose()-np.tile(value_search,(len(comp_vec),1))
         dists = np.abs(complex_dists);

flatten_list.py

@@ -1,4 +1,7 @@
-def flatten_list(list_of_lists):
+import sys
+sys.path.append('/mnt/data01/Code/')
+
+def flatten_list(list_of_lists,recursive=1):
     """
     % (C) Nick Holschuh - Amherst College - 2022 ([email protected])
     % This will take a list of lists (or an array of arrays) and flattens
@@ -7,6 +10,7 @@ def flatten_list(list_of_lists):
     % The inputs are as follows:
     %
     % flatten_list -- The list of lists to be flattened
+    % recursive -- This will keep flattening until the list stops changing
     % 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % The outputs are as follows:
@@ -18,19 +22,39 @@ def flatten_list(list_of_lists):
     """
 
     import numpy as np
+    from NDH_Tools import compare_list
 
     flat_list = []
     array_flag = 0
-    for sublist in list_of_lists:
-        if isinstance(sublist,type([])):
-            for item in sublist:
-                flat_list.append(item)
-        elif isinstance(sublist,type(np.array([]))):
-            array_flag = 1
-            for item in sublist:
-                flat_list.append(item)            
-        else:
-            flat_list.append(sublist)
+    recursive_flag = 1
+    flat_list_check = list_of_lists.copy()
+    
+    while recursive_flag == 1:
+        for sublist in list_of_lists:
+            if isinstance(sublist,type([])):
+                for item in sublist:
+                    flat_list.append(item)
+            elif isinstance(sublist,type(np.array([]))):
+                array_flag = 1
+                for item in sublist:
+                    flat_list.append(item)            
+            else:
+                flat_list.append(sublist)
 
-    flat_list = np.array(flat_list)
-    return flat_list
+        flat_list = list(np.array(flat_list))
+        
+        ######### Here we test if the lists are truly identical
+        comb_result,r1,r2 = compare_list(flat_list,flat_list_check)
+
+        if comb_result:
+            recursive_flag = 0
+        if recursive == 0:
+            recursive_flag = 0
+        
+        if recursive_flag == 1:
+            flat_list_check = flat_list.copy()
+            list_of_lists = flat_list.copy()
+            flat_list = []
+        
+    
+    return flat_list

generate_animation.py

@@ -9,9 +9,9 @@ def generate_animation(fps,title='Matplotlib Animation',comment='Matplotlib Anim
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % The inputs are as follows:
     %
-    % spacing -- the gap between sequential values
-    % fps -- the offset from 0 to round to
-    %
+    % fps -- The framerate (in frames per second)
+    % title -- This is the title of the video as stored in the file
+    % comment -- This stores some meta-text about your videoo
     % 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     % The outputs are as follows:

is2_fileselect.py

@@ -6,7 +6,7 @@
 from scipy.io import loadmat
 import inspect
 
-def is2_fileselect(xbounds, ybounds, fdir='/data/rd08/projects/IS2/ATL06', ant_or_green=0, filetype='ATL06', cycle=6, print_flag=0):
+def is2_fileselect(xbounds, ybounds, fdir='/mnt/data01/Data/IS2/ATL11/', ant_or_green=0, filetype='ATL11', cycle=6, print_flag=0):
     """
     % (C) Nick Holschuh - Amherst College -- 2022 ([email protected])
     % This function identifies the files and the segids for ICESat-2 data that falls within
@@ -17,9 +17,9 @@ def is2_fileselect(xbounds, ybounds, fdir='/data/rd08/projects/IS2/ATL06', ant_o
     %     xbounds -- two values that define the minimum x and maximum x coord (polar stereographic)
     %     ybounds -- two values that define the minumum y and maximum y coord (polar stereographic)
     %     fdir -- the absolute path to the directory where the ATL06/11 files are stored
-    %     ant_or_green -- Specifices which region range to select from
-    %     filetype -- string, currently only 'ATL06' is implemented
-    %     cycle -- the cycle number for plotting, from 1-6 [default 6]
+    %     ant_or_green -- Specifices which region range to select from [0: Ant, 1:Gre]
+    %     filetype -- string, currently only 'ATL06' and 'ATL11' are implemented
+    %     cycle -- the cycle number for plotting, from 1-6 [default 6], only required for ATL06
     %
     %%%%%%%%%%%%%%%
     % The outputs are:

loadmat.py

@@ -1,16 +1,52 @@
 import hdf5storage
-from scipy.io import loadmat
+import mat73
+import scipy
 
-def loadmat(fn):
+import sys
+sys.path.append('/mnt/data01/Code/')
+
+def loadmat(fn,varnames=None):
     """
 
 
     """
+    
+    from NDH_Tools import read_h5
 
     try:
-        data = hdf5storage.loadmat(fn)
+        if varnames == None:
+            data = mat73.loadmat(fn)
+            #print('Used the Mat73 Method')
+        else:
+            data = read_h5(fn,varnames)
+            
     except:
-        data = loadmat(fn)
-        
+        try:
+            
+            data = scipy.io.loadmat(fn,varnames,squeeze_me=True)
+            
+#            ############ This loop collapses unnecessary dimensions
+#            for i in data.keys():
+#                if type(data[i]) == type({}):
+#                    for j in data[i].keys():
+#                        if type(data[i][j]) == type({}):
+#                            pass
+#                        else:
+#                            try:
+#                                data[i][j] = np.squeeze(data[i][j])
+#                            except:
+#                                pass
+#                else:
+#                    try:
+#                        data[i] = np.squeeze(data[i])
+#                    except:
+#                        pass                       
+                        
+            
+            #print('Used the SciPy Method')
+        except:
+            print('Something is wrong with this .mat file')
+            
     return data
 
+    

regrid.py

@@ -1,5 +1,9 @@
 import numpy as np
 from scipy.interpolate import RegularGridInterpolator
+import matplotlib.pyplot as plt
+
+class CustomError(Exception):
+     pass
 
 def regrid(inx,iny,indata,newx,newy):
     """
@@ -22,6 +26,8 @@ def regrid(inx,iny,indata,newx,newy):
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     """
 
+    ############### Here we decide if we need to flip the input data to make axes increasing --
+    ############### (This is important so that we know to re-flip the output before returning)
     dx = inx[1]-inx[0]
     dy = iny[1]-iny[0]
 
@@ -39,12 +45,33 @@ def regrid(inx,iny,indata,newx,newy):
     else:
         flipy_flag = 0
 
+    ############### But we also need to know if we need to flip the comparison grid --
+    ############### (This is only important for the regridding step)        
+    dx2 = newx[1]-newx[0]
+    dy2 = newy[1]-newy[0]
+    if dx2 < 0:
+        newx = newx[::-1]
+    if dy2 < 0:
+        newy = newy[::-1]
+        
+        
     ############## Here is the actual function
-    my_interpolater = RegularGridInterpolator((iny,inx),indata)
+    my_interpolater = RegularGridInterpolator((iny,inx),indata,fill_value=np.NaN,bounds_error=False)
 
     finalx,finaly = np.meshgrid(newx,newy)
 
-    outdata = my_interpolater((finaly,finalx))
+    try:
+        outdata = my_interpolater((finaly,finalx))
+    except:
+        print('Input Grid (must have values for all output grid: [',np.min(inx),np.max(inx),'],[',np.min(iny),np.max(iny),']')
+        print('Output Grid: [',np.min(newx),np.max(newx),'],[',np.min(newy),np.max(newy),']')
+        box1x = [np.min(inx),np.max(inx),np.max(inx),np.min(inx),np.min(inx)]
+        box1y = [np.max(iny),np.max(iny),np.min(iny),np.min(iny),np.max(iny)]
+        box2x = [np.min(newx),np.max(newx),np.max(newx),np.min(newx),np.min(newx)]
+        box2y = [np.max(newy),np.max(newy),np.min(newy),np.min(newy),np.max(newy)]
+        plt.plot(box1x,box1y,'--',c='black',label='Input Grid')
+        plt.plot(box2x,box2y,'-',c='red',label='Output Grid')
+        raise CustomError("You have a bounds error for your grids")
 
     if flipy_flag == 1:
         outdata = np.flipud(outdata)

smooth_ndh.py

@@ -1,7 +1,7 @@
 from scipy import signal
 import numpy as np
 
-def smooth_ndh(input_array,smooth_window):
+def smooth_ndh(input_array,smooth_window,window_type='Boxcar'):
     """
     % (C) Nick Holschuh - Amherst College - 2022 ([email protected])
     % Convolution based smoothing function -- both 1d and 2d
@@ -28,6 +28,17 @@ def smooth_ndh(input_array,smooth_window):
         smoothing_kernel = smoothing_kernel/np.prod(smoothing_kernel.shape)
     else:
         smoothing_kernel = smooth_window
+        
+    if window_type == 'Gaussian':
+        if len(input_array.shape) == 2:
+            x = np.linspace(-4,4,smooth_window)
+            y = np.linspace(-4,4,smooth_window)
+            xmesh, ymesh = np.meshgrid(x,y);
+            smoothing_kernel = np.exp(-(np.power(xmesh, 2.) / (2 * np.power(1, 2.)) + np.power(ymesh, 2.) / (2 * np.power(1, 2.))))
+        elif len(input_array.shape) == 1:
+            x = np.linspace(-4,4,smooth_window)
+            smoothing_kernel = np.exp(-np.power(x, 2.) / (2 * np.power(1, 2.)))
+        
 
     first_out = signal.fftconvolve(input_array, smoothing_kernel, mode = 'same')
     scalar = signal.fftconvolve(np.ones(input_array.shape), smoothing_kernel, mode = 'same')

yearfrac.py

@@ -55,8 +55,9 @@ def yearfrac(date):
             years = i.astype('datetime64[Y]').astype(int) + 1970
             months = i.astype('datetime64[M]').astype(int) % 12 + 1
             days = i - i.astype('datetime64[M]') + 1
-            days = (days*1.15741e-14).astype('int')
+            days = (days).astype('int')
 
+            #print(years,months,days)
             date_output.append(float(years)+float(months)/12+float(days)/365.25)
 
     return date_output