Added the function to calculate spherical spreading corrections for radar images

Author: nholschuh

Diff for: SpreadingCorrection.nc

@@ -8,7 +8,7 @@
     if os.path.isfile(i): sys.path.append(i)
 ################################################################################################
 
-def find_pixelcoords(im_filename,original_width,original_height,im_pick_params=0):
+def find_pixelcoords(im_filename,original_width,original_height,im_pick_params=0,predefined_row_inds=[]):
     """
     % (C) Nick Holschuh - Amherst College -- 2022 ([email protected])
     %
@@ -20,14 +20,16 @@ def find_pixelcoords(im_filename,original_width,original_height,im_pick_params=0
     %
     %     im_filename - The filename of the image that you want to extract info from
     %     originial_width - The original matrix width (used to generate the image)
-    %     original__height - The original matrix height (used to generate the image
+    %     original__height - The original matrix height (used to generate the image)
     %     im_pick_params - This defines how the contours are treated. There should be four parameters per image:
     %                             ### ---- which bands to look for minima in (color of lines)
     %                             ### ---- minimum contour size (number of points)
     %                             ### ---- the aggregation type (0: average of marker, 1: horizontal bar)
     %                             ### ---- distance threshold for pixel combination or separation
     %                             ### ---- the distance calculation method (0:true or 1:vertical or 2:overweight horizontal)
     %                       0 defaults to [[0,5,0,5],[2,25,1,10]] which looks for blue points and red lines
+    %     predefined_row_inds - Some images have white pixels in the plot. This allows you to predefine
+    %                          which rows should be treated as "within the image"
     %
     %%%%%%%%%%%%%%%
     % The outputs are:
@@ -66,7 +68,11 @@ def find_pixelcoords(im_filename,original_width,original_height,im_pick_params=0
 
                 rinds = ndh.minmax(np.where(row_sum <= np.percentile(row_sum,80)-1))
                 cinds = ndh.minmax(np.where(col_sum <= np.percentile(col_sum,80)-1))
-                
+
+            ########################## If we want to figure out the number of rows from the full dataset
+            if len(predefined_row_inds) > 0:
+                rinds = predefined_row_inds
+            
             xrange = np.linspace(0,original_width,cinds[1]-cinds[0])
             yrange = np.linspace(0,original_height,rinds[1]-rinds[0])
 

Diff for: Untitled.ipynb

@@ -1,7 +1,7 @@
 import xarray as xr
 import numpy as np
 
-def make_nc(input_x,input_y,input_vars,input_varnames,num_dims=0,filename='temp.nc',description='No description provided',writefile0_or_dataset1_or_both2 = 1):
+def make_nc(input_x,input_y,input_vars,input_varnames,num_dims=0,filename='temp.nc',description='No description provided',coord_vars=['x','y'],writefile0_or_dataset1_or_both2 = 1):
     """
     % (C) Nick Holschuh - Amherst College -- 2022 ([email protected])
     %
@@ -16,6 +16,7 @@ def make_nc(input_x,input_y,input_vars,input_varnames,num_dims=0,filename='temp.
     %    input_varnames,
     %    filename='temp.nc',
     %    description='No description provided',
+    %    coord_vars=List of names for x coord and y coord, defaults to ['x','y']
     %    writefile0_or_dataset1 = 1
     %
     %%%%%%%%%%%%%%%
@@ -41,35 +42,35 @@ def make_nc(input_x,input_y,input_vars,input_varnames,num_dims=0,filename='temp.
     if num_dims == 0:
 
         if input_vars.ndim == 1:
-            data_dict[input_varnames[0]] = (['x'],input_vars)
-            coord_dict = {'x':(['x'],input_x)}
+            data_dict[input_varnames[0]] = ([coord_vars[0]],input_vars)
+            coord_dict = {coord_vars[0]:([coord_vars[0]],input_x)}
         elif input_vars.ndim == 2:
-            data_dict[input_varnames[0]] = (['y','x'],input_vars)
-            coord_dict = {'x':(['x'],input_x),'y':(['y'],input_y)}
+            data_dict[input_varnames[0]] = ([coord_vars[1],coord_vars[0]],input_vars)
+            coord_dict = {coord_vars[0]:([coord_vars[0]],input_x),coord_vars[1]:([coord_vars[1]],input_y)}
 
         ####### If there is more than one variable provided, then it assumes that
         ####### the variable data has been provided as a 3D array, with the 
         ####### third dimension defining each variable
         else:
             for ind1,i in enumerate(input_varnames):
-                data_dict[i] = (['y','x'],input_vars[:,:,ind1])
-            coord_dict = {'x':(['x'],input_x),'y':(['y'],input_y)}
+                data_dict[i] = ([coord_vars[1],coord_vars[0]],input_vars[:,:,ind1])
+            coord_dict = {coord_vars[0]:([coord_vars[0]],input_x),coord_vars[1]:([coord_vars[1]],input_y)}
 
     ######## Here we handle the case where you specify 1D data but provide
     ######## multiple data variables
     if num_dims == 1:
         for ind1,i in enumerate(input_varnames):
-            data_dict[i] = (['x'],input_vars[ind1])
-        coord_dict = {'x':(['x'],input_x)}
+            data_dict[i] = ([coord_vars[0]],input_vars[ind1])
+        coord_dict = {coord_vars[0]:([coord_vars[0]],input_x)}
 
     if num_dims == 2:
         if input_vars.ndim == 2:
-            data_dict[input_varnames[0]] = (['y','x'],input_vars)
-            coord_dict = {'x':(['x'],input_x),'y':(['y'],input_y)}
+            data_dict[input_varnames[0]] = ([coord_vars[1],coord_vars[0]],input_vars)
+            coord_dict = {coord_vars[0]:([coord_vars[0]],input_x),coord_vars[1]:([coord_vars[1]],input_y)}
         else:
             for ind1,i in enumerate(input_varnames):
-                data_dict[i] = (['y','x'],input_vars[:,:,ind1])
-            coord_dict = {'x':(['x'],input_x),'y':(['y'],input_y)}
+                data_dict[i] = ([coord_vars[1],coord_vars[0]],input_vars[:,:,ind1])
+            coord_dict = {coord_vars[0]:([coord_vars[0]],input_x),coord_vars[1]:([coord_vars[1]],input_y)}
 
     ######## Here is where the netcdf actually gets written
     ds = xr.Dataset(

Diff for: find_pixelcoords.py

@@ -10,7 +10,7 @@
 
 
 
-def process_Standard_pickedpdf(picked_files,orig_radar_dir,layer_save, cresis_flag=1, layer_save_type=1, layer_load=''):
+def process_Standard_pickedpdf(picked_files,orig_radar_dir,layer_save, cresis_flag=1, layer_save_type=1, layer_load='', find_rows_from_fullimageset = 0):
     """
     % (C) Nick Holschuh - Amherst College -- 2022 ([email protected])
     %
@@ -26,6 +26,9 @@ def process_Standard_pickedpdf(picked_files,orig_radar_dir,layer_save, cresis_fl
     %     layer_save_type=1 - For most applications, this should be set to 1, which is, save files in your current dir. 
     %                     0 - This allows you to save layer files within the cresis file_tree
     %     layer_load='' - This is not fully implemented, but it would allow you to populate existing layer files
+    %     find_rows_from_fullimageset - Setting this to 1 will search all images to figure out which rows
+    %                                   of pixels are within the plot. For use when some of the bottom of radargrams
+    %                                   is all white.
     %
     %%%%%%%%%%%%%%%
     % The outputs are:
@@ -37,6 +40,7 @@ def process_Standard_pickedpdf(picked_files,orig_radar_dir,layer_save, cresis_fl
     import glob
     import os
     import numpy as np
+    import matplotlib.pyplot as plt
 
     from PIL import Image
     import cv2
@@ -100,6 +104,8 @@ def process_Standard_pickedpdf(picked_files,orig_radar_dir,layer_save, cresis_fl
             os.makedirs(comb_deconstruct_dir)
         if not os.path.isdir(save_dir):
             os.makedirs(save_dir)
+
+        
         ##########################################################################################################
         # Part 3 #################################################################################################
         ######## Here we define the objects that need to be populated with picks
@@ -116,6 +122,42 @@ def process_Standard_pickedpdf(picked_files,orig_radar_dir,layer_save, cresis_fl
         ##########################################################################################################
         # Part 5 ##################################################################################################
         print('Starting the information extraction.')
+
+
+        
+        ########## For images that have white at the bottom, we first look through all images to figure out the appropriate 
+        ########## row indecies to use.
+
+        if find_rows_from_fullimageset == 1:
+            for ind1,frame_fn in enumerate(frame_list):
+                im_handle = Image.open(frame_fn)
+                np_frame = np.array(im_handle)
+                np_frame_dims = np_frame.shape
+                
+               #print('Frame number %0.2d' % ind1)
+               #print('np_frame shape: ',np_frame_dims)
+    
+                ########## For PNGs with RGB+components, this works best
+                if np_frame_dims[2] == 4:
+                    if ind1 == 0:
+                        im_frame = np.array(np_frame[:,:,3] != 0).astype(float)
+                    else:
+                        im_frame = im_frame + np.array(np_frame[:,:,3] != 0).astype(float)
+                ########## For greyscale+transparent, this is what you need
+                elif np_frame_dims[2] == 2:
+                    if ind1 == 0:
+                        im_frame = np.array(np_frame[:,:,1] != 0).astype(float)
+                    else:
+                        im_frame = im_frame + np.array(np_frame[:,:,1] != 0).astype(float) 
+                    
+            if 'im_frame' in locals():   
+                #plt.imshow(im_frame)
+                selected_rows = ndh.minmax(np.where(im_frame > 0)[0])
+                #print(selected_rows)
+            else: 
+                find_rows_from_fullimageset = 0
+                print('Something went wrong with the full image-set. Defaulting to local row search.')
+
         ########## Here we actually load the images and extract pixel coordinate information
         error_frames = []
         good_frames = []
@@ -137,8 +179,11 @@ def process_Standard_pickedpdf(picked_files,orig_radar_dir,layer_save, cresis_fl
                     original_height = height_index['index'][0]+100
                 elif crop == 'nocrop':
                     original_height = len(frame_data['Time'])
-                    
-                picks = ndh.find_pixelcoords(frame_fn,original_width,original_height,im_pick_params=[[2,25,1,10,1]])  
+
+                if find_rows_from_fullimageset == 0:
+                    picks = ndh.find_pixelcoords(frame_fn,original_width,original_height,im_pick_params=[[2,25,1,10,1]])  
+                else:
+                    picks = ndh.find_pixelcoords(frame_fn,original_width,original_height,im_pick_params=[[2,25,1,10,1]], predefined_row_inds=selected_rows)
 
 
 

Diff for: make_nc.py

@@ -0,0 +1,25 @@
+
+def rotation(input_vec,rotation_angle,dims=2,axis=0):
+
+    import numpy as np
+
+    if dims == 2:
+        if axis == 0:
+            ########## A rotation around x
+            rotation = np.array([[np.cos(rotation_angle),-np.sin(rotation_angle)],[np.sin(rotation_angle),np.cos(rotation_angle)]])
+        
+    if dims == 3:
+        if axis == 0:
+            ########## A rotation around x
+            rotation = np.array([[1,0,0],[0,np.cos(ref_axial_tilt),-np.sin(ref_axial_tilt)],[0,np.sin(ref_axial_tilt),np.cos(ref_axial_tilt)]])
+        
+        if axis == 1:
+            ########## A rotation around y
+            rotation = np.array([[np.cos(ref_axial_tilt),0,np.sin(ref_axial_tilt)],[0,1,0],[-np.sin(ref_axial_tilt),0,np.cos(ref_axial_tilt)]])
+        
+        if axis == 2:
+            ########## A rotation around z
+            rotation = np.array([[np.cos(ref_axial_tilt),-np.sin(ref_axial_tilt),0],[np.sin(ref_axial_tilt),np.cos(ref_axial_tilt),0],[0,0,1]])
+
+    input_vec_rotated = np.matmul(input_vec,rotation)
+    return input_vec_rotated

Diff for: process_Standard_pickedpdf.py

@@ -39,7 +39,8 @@ def smooth_ndh(input_array,smooth_window,window_type='Boxcar'):
             x = np.linspace(-4,4,smooth_window)
             smoothing_kernel = np.exp(-np.power(x, 2.) / (2 * np.power(1, 2.)))
 
-    
+
+    input_array[input_array == np.NaN] = 0
     first_out = signal.fftconvolve(input_array, smoothing_kernel, mode = 'same')
     scalar = signal.fftconvolve(np.ones(input_array.shape), smoothing_kernel, mode = 'same')
     out_data = first_out / scalar;