Merge branch 'feat-step3-intg' of https://github.com/brown-ccv/icesat…

…2-tracks into feat-step3-intg

analysis/SB04_2d_wavefield_emulator.py

@@ -1,6 +1,5 @@
 
 import os, sys
-#execfile(os.environ['PYTHONSTARTUP'])
 
 """
 This file open a ICEsat2 track applied filters and corections and returns smoothed photon heights on a regular grid in an .nc file.
@@ -10,8 +9,6 @@
     exec(open(os.environ['PYTHONSTARTUP']).read())
     exec(open(STARTUP_2021_IceSAT2).read())
 
-    #%matplotlib inline
-
     import ICEsat2_SI_tools.convert_GPS_time as cGPS
     import h5py
     import ICEsat2_SI_tools.io as io
@@ -22,36 +19,31 @@
     import spicke_remover
     import datetime
     import concurrent.futures as futures
-    #import s3fs
-    # %%
+
     track_name, batch_key, test_flag = io.init_from_input(sys.argv) # loads standard experiment
-    #track_name, batch_key, test_flag = '20190605061807_10380310_004_01', 'SH_batch01', False
-    #track_name, batch_key, test_flag = '20190601094826_09790312_004_01', 'SH_batch01', False
-    #track_name, batch_key, test_flag = '20190207111114_06260210_004_01', 'SH_batch02', False
+
     track_name, batch_key, test_flag = '20190219073735_08070210_004_01', 'SH_batch02', False
 
 
 
 
     #print(track_name, batch_key, test_flag)
     hemis, batch = batch_key.split('_')
-    #track_name= '20190605061807_10380310_004_01'
     ATlevel= 'ATL03'
 
     save_path   = mconfig['paths']['work'] + '/B03_spectra_'+hemis+'/'
     save_name   = 'B03_'+track_name
 
-    #plot_path   = mconfig['paths']['plot'] + '/'+hemis+'/'+batch_key+'/' + track_name + '/B_spectra/'
     plot_path   = mconfig['paths']['plot'] + '/phase_fitting_fake/2D_fake/'
     MT.mkdirs_r(plot_path)
     MT.mkdirs_r(save_path)
     bad_track_path =mconfig['paths']['work'] +'bad_tracks/'+ batch_key+'/'
-# %%
+
 if __name__ == '__main__':
     all_beams   = mconfig['beams']['all_beams']
     high_beams  = mconfig['beams']['high_beams']
     low_beams   = mconfig['beams']['low_beams']
-    #Gfilt   = io.load_pandas_table_dict(track_name + '_B01_regridded', load_path) # rhis is the rar photon data
+
 
     load_path   = mconfig['paths']['work'] +'/B01_regrid_'+hemis+'/'
     Gd      = io.load_pandas_table_dict(track_name + '_B01_binned' , load_path)  #
@@ -70,30 +62,21 @@
     Gspec.coords['T'] = 1/Gspec.coords['f']
     Gspec=Gspec.swap_dims({'k': 'f'})
 
-    #Gspec.spectral_power_optm.sel(beam='weighted_mean').plot()
-    # %%
-    #k_lim= 0.02
+
     A, B = Gspec.sel(beam= 'gt2r').Y_model_hat , Gspec.sel(beam= 'gt2l').Y_model_hat
 
     r_ave_kargs={'x':2, 'f':10, 'center':True, 'min_periods':2}
     r_ave_kargs2={'f':10, 'center':True, 'min_periods':2}
-    #(abs(B) - abs(A)).plot()
+
     S_aa = (A*A.conj()).real
     S_bb = (B*B.conj()).real
-    # co_spec = (A.conj() *B) /S_aa/S_bb
-    # abs(co_spec).plot()
+
     co_spec = (A.conj() *B).rolling(**r_ave_kargs).mean()
     np.log(abs(co_spec)).plot(levels=np.arange(-2, 3, 0.1))
 
     (abs(co_spec)).plot(levels=np.exp(np.arange(-3, 2, 0.1)))
 
     abs(co_spec).mean('x').plot()
-
-#(abs(co_spec)/(S_aa *S_bb).rolling(**r_ave_kargs).mean()).plot()
-#
-# (abs(A.conj() *B)/(S_aa *S_bb)).rolling(**r_ave_kargs).mean()[:,:].plot()
-
-# %%
 if __name__ == '__main__':
     L1 = 50
     k1 = 2* np.pi /L1
@@ -111,11 +94,9 @@
     alpha = 35
     kk, ll = np.cos(alpha * np.pi/180) * np.array([0.9*k1, k1, 1.1* k1]), np.sin(alpha * np.pi/180) * np.array([0.9* k1, 1*k1, 1.1* k1])
     M_k, M_l = kk.size, ll.size
-    #y =np.sin(k1* x) + np.sin(k2* x)
     kk_mesh, ll_mesh = np.meshgrid(kk, ll)
     kk_mesh, ll_mesh = kk_mesh.reshape(kk_mesh.size), ll_mesh.reshape(ll_mesh.size)
-    G = np.cos(np.outer(XX, kk_mesh) + np.outer(YY, ll_mesh)).T# + np.sin(np.outer(XX, kk_mesh) + np.outer(YY, ll_mesh)).T
-    #G = np.vstack([ np.cos(np.outer(x, k) + np.outer(y, l)).T ,  np.sin(np.outer(x, k)  + np.outer(y, l) ).T ] ).T
+    G = np.cos(np.outer(XX, kk_mesh) + np.outer(YY, ll_mesh)).T
     G.shape
 
     plt.contourf(x, y, G.sum(0).reshape(Ny, Nx) )
@@ -124,7 +105,6 @@
  # %% radial coordincates
 
 def gaus_2d(x, y, pos_tuple, sigma_g ):
-    #grid = ( (XX - pos_tuple[0])  * (YY - pos_tuple[1]) )
     gx = np.exp(-0.5 * (x - pos_tuple[0])**2 /sigma_g**2  )
     gy = np.exp(-0.5 * (y - pos_tuple[1])**2 /sigma_g**2  )
     return  np.outer(gx ,  gy).T
@@ -140,8 +120,6 @@ def gaus_2d(x, y, pos_tuple, sigma_g ):
     plt.contourf(k_range, l_range, gaus_lk  )
     plt.axis('equal')
 
-# %%
-
 k_0 = 0.03
 l_0 = 0
 dk = 0.01
@@ -203,7 +181,6 @@ def gaus_2d_mesh(XX,YY, pos_tuple, sigma_g ):
     plt.axis('equal')
 
 
-# %% radial coodinates
 def get_stancils_polar( amp, angle_rad,  size=1, dk = 0.01, mesh = True, plot_flag = True, amp_std= None, random=True):
     """
     inputs:
@@ -263,10 +240,6 @@ def get_stancils_polar( amp, angle_rad,  size=1, dk = 0.01, mesh = True, plot_fl
 
 if __name__ == '__main__':
     font_for_pres()
-    #k_mesh.shape
-    #for angle in np.arange(-80, 80+20, 40):
-    #for phase in np.arange(0, 2*np.pi, np.pi/3):
-    #for k_abs in np.arange(0.01, 0.09, 0.01):
 
     angle =45
     amp =1
@@ -276,24 +249,16 @@ def get_stancils_polar( amp, angle_rad,  size=1, dk = 0.01, mesh = True, plot_fl
     for dk in np.arange(0.005, 0.03, 0.002):
 
         for size in [2, 5, 10, 30, 50, 100, 200]:
-            #for angle in np.arange(-80, 80+20, 40):
-            #for k_abs in np.arange(0.01, 0.09, 0.01):
-            #for phase in np.arange(0, 2*np.pi, np.pi/3):
-
-
             F = M.figure_axis_xy(8, 3, view_scale = 0.5, container =False)
             plt.suptitle('k_abs=' + str(k_abs) +' angle=' + str(angle) + ' size=' + str(size) +' dk=' + str(dk) )
             ax = plt.subplot(1, 2, 1)
             k_list, l_list, amp_weights, stancil_shape = get_stancils_polar(k_abs, angle * np.pi/180, size=size, dk = dk, mesh = True , plot_flag= True, random = True)
-
             circle1 = plt.Circle((k_list[0], l_list[0]), dk, color='b', fill=False)
             ax.add_patch(circle1)
 
-
             k_noise, l_noise, amp_noise, stancil_shape = get_stancils_polar(0.8, 0 * np.pi/180, size=20, dk = 0.3, mesh = True , plot_flag= True, random = True)
             amp_noise = (amp_noise *0+1) * 0
 
-
             plt.xlim(0, 0.1)
             #plt.axis('equal')
             plt.ylim(-0.1, 0.1)
@@ -308,10 +273,6 @@ def get_stancils_polar( amp, angle_rad,  size=1, dk = 0.01, mesh = True, plot_fl
             amp_all = np.concatenate([amp_weights, amp_noise])
             amp_all.shape
             G = np.vstack([ np.cos(np.outer(XX, k_all) + np.outer(YY, l_all)).T ,  np.sin(np.outer(XX, k_all) + np.outer(YY, l_all)).T ] ).T
-            #G = np.vstack([ np.cos(np.outer(x, k) + np.outer(y, l)).T             ,  np.sin(np.outer(x, k)  + np.outer(y, l) ).T ] ).T
-
-            #phase1 = np.random.rand(1, amp_list.size) *  np.pi*2
-            #phase = np.arange(0, amp_list.size) *  np.pi/2
 
             b = np.hstack([ np.cos(phase)*amp_all, np.sin(phase) *amp_all]).squeeze() * amp
             z_model = (G @ b).reshape(Ny, Nx)
@@ -323,9 +284,3 @@ def get_stancils_polar( amp, angle_rad,  size=1, dk = 0.01, mesh = True, plot_fl
             F.save_light(path = plot_path, name = 'fake_2d_dk' +str(dk) +'_s' + str(int(size)))
             plt.show()
 
-
-    # class twoD_wave_packets(object):
-    #     def __init__(self):
-
-
-    # %%

analysis_db/B05_define_angle.py

@@ -280,7 +280,7 @@ def weighted_means(data, weights, x_angle, color='k'):
     exit()
 
 # %%
-class plot_polarspectra(object):
+class plot_polarspectra:
         def __init__(self,k, thetas, data, data_type='fraction' ,lims=None,  verbose=False):
 
             """

src/icesat2_tracks/ICEsat2_SI_tools/angle_optimizer.py

@@ -122,7 +122,7 @@ def simple_log_panelty(x, x0, sigma):
         return  - 0.5 * (cost_sqrt/tot_var + np.log(tot_var) ).sum() + prior_weight * penalties
 
 
-class sample_with_mcmc(object):
+class sample_with_mcmc:
     """
     sample a 2nd surface using mcmc and other methods. its made for getting a quick estimate!
 

src/icesat2_tracks/ICEsat2_SI_tools/io.py

@@ -70,13 +70,10 @@ def ID_to_str(ID_name):
     date
     return IDs[0] +' ' +date +' granule: ' + IDs[2]
 
-class case_ID(object):
+class case_ID:
     """docstring for case_ID"""
     def __init__(self, track_name):
         import re
-        super(case_ID, self).__init__()
-
-        #track_name_pattern = r'(\d{4})(\d{2})(\d{2})(\d{2})(\d{2})(\d{2})_(\d{4})(\d{2})(\d{2})_(\d{3})_(\d{2})'
         track_name_pattern = r'(\D{2}|\d{2})_?(\d{4})(\d{2})(\d{2})(\d{2})?(\d{2})?(\d{2})?_(\d{4})(\d{2})(\d{2})_?(\d{3})?_?(\d{2})?'
         case_ID_pattern = r'(\d{4})(\d{2})(\d{2})_(\d{4})(\d{2})(\d{2})'
 

src/icesat2_tracks/ICEsat2_SI_tools/spectral_estimates.py

@@ -376,7 +376,7 @@ def get_lon_lat_coords(stancil):
     return coord_positions
 
 
-class wavenumber_spectrogram(object):
+class wavenumber_spectrogram:
     def __init__(self, x_grid, data, Lpoints, ov=None, window=None):
         """
         returns a wavenumber spectrogram with the resolution L-ov
@@ -512,7 +512,7 @@ def get_stancil_var_apply(stancil):
             self.G.attrs['mean_variance_detrended_chunks']  = np.array(stancil_vars).mean()
             self.G.attrs['mean_variance_pwelch_spectrum']   = self.calc_var()
 
-class wavenumber_spectrogram_LS_even(object):
+class wavenumber_spectrogram_LS_even:
     def __init__(self, x, data, L, waven_method = 'fftX2' , dy=None ,  ov=None, window=None, kjumps=1):
         """
         returns a wavenumber spectrogram with the resolution L-ov
@@ -618,7 +618,7 @@ def parceval(self, add_attrs=True ):
     def mean_spectral_error(self, confidence = 0.95 ):
         return wavenumber_spectrogram.mean_spectral_error(self, confidence= confidence )
 
-class wavenumber_spectrogram_LS(object):
+class wavenumber_spectrogram_LS:
     def __init__(self, x, data, L, dx, dy = None, waven_method = 'fftX2', ov=None, window=None):
         """
         returns a wavenumber spectrogram with the resolution L-ov
@@ -925,7 +925,7 @@ def mean_spectral_error(self, mask=None, confidence = 0.95 ):
 
 
 # class for getting standard Pwelch spectrum. old version, deprechiate
-class wavenumber_pwelch(object):
+class wavenumber_pwelch:
     def __init__(self,data, x, L, ov=None, window=None, save_chunks=False, plot_chunks=False):
         """
         returns a wavenumber spectrum using the pwelch method
@@ -1098,7 +1098,7 @@ def calc_var(self):
 
 # %% optimze spectral variance
 
-class conserve_variance(object):
+class conserve_variance:
     def __init__(self,Z, freq, data, nan_mask= None):
 
         """

src/icesat2_tracks/analysis_fake_data/X03_MCMC_surface_smapling.py

@@ -41,7 +41,7 @@ def cost(x, y):
     z = 4+ np.sin(4*  2 * np.pi *x/Lx) +  np.sin( 3 * np.pi *x/Lx - np.pi/5) + np.cos(1*  2 * np.pi *y/Ly) +  np.sin( 3 * np.pi *y/Ly - np.pi/3)
     return z**2
 
-class sample_with_mcmc(object):
+class sample_with_mcmc:
     """
     sample a 2nd surface using mcmc. its make for getting a quick estimate!
     """

src/icesat2_tracks/local_modules/m_colormanager_ph3.py

@@ -50,7 +50,7 @@ def json_load(name, path, verbose=False):
 
 
 
-class color(object):
+class color:
         def __init__(self, path=None, name=None):
             self.white=(1,1,1)
             if (path is not None) & (name is not None):

src/icesat2_tracks/local_modules/m_general_ph3.py

@@ -27,7 +27,7 @@
 import matplotlib.pyplot as plt
 
 
-class color(object):
+class color:
         def __init__(self):
 
             self.red=(203/255, 32/255, 39/255)
@@ -61,7 +61,7 @@ def show(self):
 
 # funny massage
 
-class figure_axis_xy(object):
+class figure_axis_xy:
         """define standart  XY Plot with reduced grafics"""
 
         def __init__(self,x_size=None,y_size=None,view_scale=None, size_tuple=None , fig_scale=None, container=False, dpi=180):
@@ -171,7 +171,7 @@ class subplot_routines(figure_axis_xy):
         def __init__(self, ax):
             self.ax=ax
 
-class plot_sprecta(object):
+class plot_sprecta:
             def __init__(self,fs, Xdata,sample_unit=None,data_unit=None):
 
                 self.fs=fs
@@ -261,7 +261,7 @@ def power(self, Color='b', fig_scale=2, fax='f'):
                 self.F.make_clear()
                 plt.grid()
 
-class plot_periodogram(object):
+class plot_periodogram:
             def __init__(self,time,fs, data,clevs=None, sample_unit=None, data_unit=None,
                         ylim=None, time_unit=None, cmap=None):
                 self.fs=fs[1:]
@@ -696,7 +696,7 @@ def set_xaxis_to_days(self, **kwargs):
                 set_timeaxis_days(self.ax, **kwargs)
 
 
-class plot_polarspectra(object):
+class plot_polarspectra:
         def __init__(self,f, thetas, data,unit=None, data_type='fraction' ,lims=None,  verbose=False):
 
             self.f=f
@@ -1256,7 +1256,7 @@ def spickes_to_mean(ts, nloop=None, spreed=1, gaussian=True):
 
 ## Composites
 
-class composite_data(object):
+class composite_data:
         def __init__(self,var, index_weight=None):
             #print(var.shape)
             self.composites=var
@@ -1285,7 +1285,7 @@ def bootstrap(self, ci=[2.5, 50, 97.5], reps=1000,):
             yb = 1/np.arange(1, n+1)[:, None] * np.cumsum(xb, axis=0)
             upper, lower = np.percentile(yb, [2.5, 97.5], axis=1)
 
-class comp_iter(object):
+class comp_iter:
         def __init__(self, span, dt=None, unit=None):
             self.span=list(span)
             for s in self.span:
@@ -1306,7 +1306,7 @@ def __init__(self, span, dt=None, unit=None):
 
                 self.time_iter_string=time_str
 
-class composite(object):
+class composite:
         def __init__(self,index, time=None, weigthing=False, span=None):
             """ Initial Class for bulding composite based on:
                 index       position in the time vector 'time'

src/icesat2_tracks/local_modules/m_spectrum_ph3.py

@@ -91,7 +91,7 @@ def spicke_remover(data, nstd=20.0, spreed=500.0, max_loops=10.0 , verbose=False
 
     return data2 , act_flag
 
-class Spectrum(object):
+class Spectrum:
     """ A class that represents a single realization of
             the one-dimensional spectrum  of a given field phi """
 
@@ -149,7 +149,7 @@ def parceval(self):
 
 
 
-class moments(object):
+class moments:
     def __init__(self,data_org,dt, L=None, ov=None,window=None, save_chunks=False, plot_chunks=False, prewhite=None):
         """
         This function calculates the spectral moments from a station (buoy, GPS, or seismic station) that measures
@@ -396,7 +396,7 @@ def cal_MEM(self, theta=None, flim=(0.01, .5)):
 
 
 
-class pwelch(object):
+class pwelch:
     def __init__(self,data,dt,L=None, ov=None,window=None, save_chunks=False, plot_chunks=False, periodogram=False, prewhite=None):
         """
         prewhite    None(default)
@@ -800,7 +800,7 @@ def save_data(self, path=None, S=None):
         save_file(P, path)
 
 
-class save_data_periodogram(object):
+class save_data_periodogram:
     def __init__(self,P, S=None):
         self.meta=S.meta if S is not None else ''
         self.data_unit=S.unit if S is not None else ''