diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml index fcf209ef..bb8aebd2 100644 --- a/.github/workflows/run_tests.yml +++ b/.github/workflows/run_tests.yml @@ -24,4 +24,4 @@ jobs: - name: install pytest run: pip install pytest pytest-xdist - name: Run pytest - run: pytest -n 2 tests + run: pytest tests diff --git a/analyis_publish/PB03_plot_spectra_3x3_ov.py b/analyis_publish/PB03_plot_spectra_3x3_ov.py index 0095fd45..7cd52025 100644 --- a/analyis_publish/PB03_plot_spectra_3x3_ov.py +++ b/analyis_publish/PB03_plot_spectra_3x3_ov.py @@ -152,7 +152,7 @@ def plot_wavenumber_spectrogram(ax, Gi, clev, plot_photon_density=True , cmap=No #Gmean = G_gFT_wmean['gFT_PSD_data'].rolling(k=5, center=True).mean() -Gmean = gFT.rebin(G_gFT_wmean['gFT_PSD_data'], 10) +Gmean,_ = gFT.rebin(G_gFT_wmean['gFT_PSD_data'], 10) #Gmean = Gmean.where(~np.isnan(Gmean), 0) diff --git a/pyproject.toml b/pyproject.toml index 4aa96da3..86e7e7f4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -116,6 +116,7 @@ dependencies = [ # Optional "lmfit", "tables", "emcee==3.1.4", + "siphon", "siphon" ] diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/io.py b/src/icesat2_tracks/ICEsat2_SI_tools/io.py index 521e305c..e84f9751 100644 --- a/src/icesat2_tracks/ICEsat2_SI_tools/io.py +++ b/src/icesat2_tracks/ICEsat2_SI_tools/io.py @@ -1,41 +1,50 @@ +from sliderule import icesat2 +from icesat2_tracks.ICEsat2_SI_tools import sliderule_converter_tools as sct def init_from_input(arguments): - if (len(arguments) <= 1) | ('-f' in set(arguments) ) : + """ + Initializes the variables track_name, batch_key, and test_flag based on the input arguments. + + Parameters: + arguments (list): A list of input arguments. + + Returns: + tuple: A tuple containing the values of track_name, batch_key, and test_flag. + """ - track_name='20190605061807_10380310_004_01' - batch_key='SH_batch01' + if (len(arguments) <= 1) | ("-f" in set(arguments)): + track_name = "20190605061807_10380310_004_01" + batch_key = "SH_batch01" test_flag = True - print('use standard values') + print("use standard values") else: + track_name = arguments[1] + batch_key = arguments[2] + # $(hemisphere) $(coords) $(config) - track_name=arguments[1] - batch_key =arguments[2] - #$(hemisphere) $(coords) $(config) + print("read vars from file: " + str(arguments[1])) - print('read vars from file: ' +str(arguments[1]) ) - - if (len(arguments) >= 4): - if arguments[3] == 'True': + if len(arguments) >= 4: + if arguments[3] == "True": test_flag = True - elif arguments[3] == 'False': + elif arguments[3] == "False": test_flag = False else: - test_flag= arguments[3] + test_flag = arguments[3] - print('test_flag found, test_flag= '+str(test_flag) ) + print("test_flag found, test_flag= " + str(test_flag)) else: - test_flag=False + test_flag = False print(track_name) - - print('----- batch ='+ batch_key) - print('----- test_flag: ' + str(test_flag)) + print("----- batch =" + batch_key) + print("----- test_flag: " + str(test_flag)) return track_name, batch_key, test_flag -def init_data(ID_name, batch_key, ID_flag, ID_root, prefix ='A01b_ID'): +def init_data(ID_name, batch_key, ID_flag, ID_root, prefix="A01b_ID"): """ Takes inputs and retrieves the ID, track_names that can be loaded, hemis, batch inputs: are the outputs from init_from_input, specifically @@ -49,103 +58,150 @@ def init_data(ID_name, batch_key, ID_flag, ID_root, prefix ='A01b_ID'): """ print(ID_name, batch_key, ID_flag) - hemis, batch = batch_key.split('_') + hemis, batch = batch_key.split("_") if ID_flag: - ID_path = ID_root +'/'+batch_key+'/'+prefix+'/' - ID = json_load( prefix +'_'+ID_name, ID_path ) - track_names = ID['tracks']['ATL03'] + ID_path = ID_root + "/" + batch_key + "/" + prefix + "/" + ID = json_load(prefix + "_" + ID_name, ID_path) + track_names = ID["tracks"]["ATL03"] else: - track_names = ['ATL03_'+ID_name] - ID = ID_name + track_names = ["ATL03_" + ID_name] + ID = ID_name return ID, track_names, hemis, batch + def ID_to_str(ID_name): from datetime import datetime - IDs = ID_name.split('_') - date = datetime.strptime(IDs[1], '%Y%m%d').strftime('%Y-%m-%d')#.strftime('%m/%d/%Y') + IDs = ID_name.split("_") + date = datetime.strptime(IDs[1], "%Y%m%d").strftime( + "%Y-%m-%d" + ) # .strftime('%m/%d/%Y') date - return IDs[0] +' ' +date +' granule: ' + IDs[2] + return IDs[0] + " " + date + " granule: " + IDs[2] + + +def get_gdf(ATL03_track_name, params_yapc, maximum_height): + print("Retrieving data from sliderule ...") + gdf = icesat2.atl06p(params_yapc, resources=[ATL03_track_name]) + + if gdf.empty: + raise Exception("Empty Geodataframe. No data could be retrieved.") + + print("Initial data retrieved. Correcting and removing height ...") + gdf = sct.correct_and_remove_height(gdf, maximum_height) + return gdf + class case_ID: """docstring for case_ID""" + def __init__(self, track_name): import re - 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})' - track_name_rx = re.compile(track_name_pattern) - self.hemis,self.YY,self.MM,self.DD,self.HH,self.MN,self.SS,self.TRK,self.CYC,self.GRN,self.RL,self.VRS = track_name_rx.findall(track_name).pop() + 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})" - if self.hemis == '01': - self.hemis = 'NH' - elif self.hemis == '02': - self.hemis = 'SH' + track_name_rx = re.compile(track_name_pattern) + ( + self.hemis, + self.YY, + self.MM, + self.DD, + self.HH, + self.MN, + self.SS, + self.TRK, + self.CYC, + self.GRN, + self.RL, + self.VRS, + ) = track_name_rx.findall(track_name).pop() + + if self.hemis == "01": + self.hemis = "NH" + elif self.hemis == "02": + self.hemis = "SH" else: self.hemis = self.hemis - #self.hemis = hemis + # self.hemis = hemis self.set() self.track_name_init = track_name def set(self): - block1 = (self.YY,self.MM,self.DD) - block2 = (self.TRK,self.CYC,self.GRN) + block1 = (self.YY, self.MM, self.DD) + block2 = (self.TRK, self.CYC, self.GRN) - self.ID = self.hemis+'_'+''.join(block1) +'_'+ ''.join(block2) + self.ID = self.hemis + "_" + "".join(block1) + "_" + "".join(block2) return self.ID def get_granule(self): - return ''.join((self.TRK,self.CYC,self.GRN)) + return "".join((self.TRK, self.CYC, self.GRN)) def set_dummy(self): - block1 = (self.YY,self.MM,self.DD) - block2 = (self.TRK,self.CYC,self.GRN) + block1 = (self.YY, self.MM, self.DD) + block2 = (self.TRK, self.CYC, self.GRN) - self.ID_dummy = ''.join(block1) +'_'+ ''.join(block2) + self.ID_dummy = "".join(block1) + "_" + "".join(block2) return self.ID_dummy def set_ATL03_trackname(self): - - block1 = (self.YY,self.MM,self.DD) - block1b = (self.HH,self.MN,self.SS) - block2 = (self.TRK,self.CYC,self.GRN) - if self.RL is '': + block1 = (self.YY, self.MM, self.DD) + block1b = (self.HH, self.MN, self.SS) + block2 = (self.TRK, self.CYC, self.GRN) + if self.RL is "": raise ValueError("RL not set") - if self.VRS is '': + if self.VRS is "": raise ValueError("VRS not set") - block3 = (self.RL,self.VRS) + block3 = (self.RL, self.VRS) - self.ID_ATL03 = ''.join(block1) +''.join(block1b) +'_'+ ''.join(block2) +'_'+ '_'.join(block3) + self.ID_ATL03 = ( + "".join(block1) + + "".join(block1b) + + "_" + + "".join(block2) + + "_" + + "_".join(block3) + ) return self.ID_ATL03 def set_ATL10_trackname(self): - - block1 = (self.YY,self.MM,self.DD) - block1b = (self.HH,self.MN,self.SS) - block2 = (self.TRK,self.CYC, '01') # granule is alwasy '01' for ATL10 - if self.RL is '': + block1 = (self.YY, self.MM, self.DD) + block1b = (self.HH, self.MN, self.SS) + block2 = (self.TRK, self.CYC, "01") # granule is alwasy '01' for ATL10 + if self.RL is "": raise ValueError("RL not set") - if self.VRS is '': + if self.VRS is "": raise ValueError("VRS not set") - block3 = (self.RL,self.VRS) + block3 = (self.RL, self.VRS) - if self.hemis == 'NH': - hemis = '01' - elif self.hemis == 'SH': - hemis = '02' + if self.hemis == "NH": + hemis = "01" + elif self.hemis == "SH": + hemis = "02" else: hemis = self.hemis - self.ID_ATL10 = hemis+'_'+''.join(block1) +''.join(block1b) +'_'+ ''.join(block2) +'_'+ '_'.join(block3) + self.ID_ATL10 = ( + hemis + + "_" + + "".join(block1) + + "".join(block1b) + + "_" + + "".join(block2) + + "_" + + "_".join(block3) + ) return self.ID_ATL10 -def nsidc_icesat2_get_associated_file(file_list, product, build=True, username=None, password=None): +def nsidc_icesat2_get_associated_file( + file_list, product, build=True, username=None, password=None +): """ THis method returns assocociated files names and paths for files given in file_list for the "product" ICEsat2 product @@ -163,34 +219,39 @@ def nsidc_icesat2_get_associated_file(file_list, product, build=True, username=N import posixpath import os import icesat2_toolkit.utilities - AUXILIARY=False - #product='ATL03' - DIRECTORY= None - FLATTEN=False - TIMEOUT=120 - MODE=0o775 - #file_list = ['ATL07-01_20210301023054_10251001_005_01'] + + AUXILIARY = False + # product='ATL03' + DIRECTORY = None + FLATTEN = False + TIMEOUT = 120 + MODE = 0o775 + # file_list = ['ATL07-01_20210301023054_10251001_005_01'] if build and not (username or password): - urs = 'urs.earthdata.nasa.gov' - username,login,password = netrc.netrc().authenticators(urs) - #-- build urllib2 opener and check credentials + urs = "urs.earthdata.nasa.gov" + username, login, password = netrc.netrc().authenticators(urs) + # -- build urllib2 opener and check credentials if build: - #-- build urllib2 opener with credentials + # -- build urllib2 opener with credentials icesat2_toolkit.utilities.build_opener(username, password) - #-- check credentials + # -- check credentials icesat2_toolkit.utilities.check_credentials() parser = lxml.etree.HTMLParser() - #-- remote https server for ICESat-2 Data - HOST = 'https://n5eil01u.ecs.nsidc.org' - #-- regular expression operator for extracting information from files - rx = re.compile(r'(processed_)?(ATL\d{2})(-\d{2})?_(\d{4})(\d{2})(\d{2})' - r'(\d{2})(\d{2})(\d{2})_(\d{4})(\d{2})(\d{2})_(\d{3})_(\d{2})') - #-- regular expression pattern for finding specific files - regex_suffix = '(.*?)$' if AUXILIARY else '(h5)$' - remote_regex_pattern = (r'{0}(-\d{{2}})?_(\d{{4}})(\d{{2}})(\d{{2}})' - r'(\d{{2}})(\d{{2}})(\d{{2}})_({1})({2})({3})_({4})_(\d{{2}})(.*?).{5}') + # -- remote https server for ICESat-2 Data + HOST = "https://n5eil01u.ecs.nsidc.org" + # -- regular expression operator for extracting information from files + rx = re.compile( + r"(processed_)?(ATL\d{2})(-\d{2})?_(\d{4})(\d{2})(\d{2})" + r"(\d{2})(\d{2})(\d{2})_(\d{4})(\d{2})(\d{2})_(\d{3})_(\d{2})" + ) + # -- regular expression pattern for finding specific files + regex_suffix = "(.*?)$" if AUXILIARY else "(h5)$" + remote_regex_pattern = ( + r"{0}(-\d{{2}})?_(\d{{4}})(\d{{2}})(\d{{2}})" + r"(\d{{2}})(\d{{2}})(\d{{2}})_({1})({2})({3})_({4})_(\d{{2}})(.*?).{5}" + ) # rx = re.compile(r'(processed_)?(ATL\d{2})(-\d{2})?_(\d{4})(\d{2})(\d{2})' # r'(\d{2})(\d{2})(\d{2})_(\d{4})(\d{2})(\d{2})_(\d{3})_(\d{2})(.*?).h5$') @@ -199,69 +260,70 @@ def nsidc_icesat2_get_associated_file(file_list, product, build=True, username=N # remote_regex_pattern = (r'{0}(-\d{{2}})?_(\d{{4}})(\d{{2}})(\d{{2}})' # r'(\d{{2}})(\d{{2}})(\d{{2}})_({1})({2})({3})_({4})_(\d{{2}})(.*?).{5}') - #-- build list of remote files, remote modification times and local files + # -- build list of remote files, remote modification times and local files original_files = [] remote_files = [] remote_mtimes = [] local_files = [] - remote_names =[] + remote_names = [] for input_file in file_list: - #print(input_file) - #-- extract parameters from ICESat-2 ATLAS HDF5 file name - SUB,PRD,HEM,YY,MM,DD,HH,MN,SS,TRK,CYC,GRN,RL,VRS = \ - rx.findall(input_file).pop() - #-- get directories from remote directory - product_directory = '{0}.{1}'.format(product,RL) - sd = '{0}.{1}.{2}'.format(YY,MM,DD) - PATH = [HOST,'ATLAS',product_directory,sd] - #-- local and remote data directories - remote_dir=posixpath.join(*PATH) - temp=os.path.dirname(input_file) if (DIRECTORY is None) else DIRECTORY - local_dir=os.path.expanduser(temp) if FLATTEN else os.path.join(temp,sd) - #-- create output directory if not currently existing + # print(input_file) + # -- extract parameters from ICESat-2 ATLAS HDF5 file name + SUB, PRD, HEM, YY, MM, DD, HH, MN, SS, TRK, CYC, GRN, RL, VRS = rx.findall( + input_file + ).pop() + # -- get directories from remote directory + product_directory = "{0}.{1}".format(product, RL) + sd = "{0}.{1}.{2}".format(YY, MM, DD) + PATH = [HOST, "ATLAS", product_directory, sd] + # -- local and remote data directories + remote_dir = posixpath.join(*PATH) + temp = os.path.dirname(input_file) if (DIRECTORY is None) else DIRECTORY + local_dir = os.path.expanduser(temp) if FLATTEN else os.path.join(temp, sd) + # -- create output directory if not currently existing # if not os.access(local_dir, os.F_OK): # os.makedirs(local_dir, MODE) - #-- compile regular expression operator for file parameters - args = (product,TRK,CYC,GRN,RL,regex_suffix) + # -- compile regular expression operator for file parameters + args = (product, TRK, CYC, GRN, RL, regex_suffix) R1 = re.compile(remote_regex_pattern.format(*args), re.VERBOSE) - #-- find associated ICESat-2 data file - #-- find matching files (for granule, release, version, track) - colnames,collastmod,colerror=icesat2_toolkit.utilities.nsidc_list(PATH, - build=False, - timeout=TIMEOUT, - parser=parser, - pattern=R1, - sort=True) + # -- find associated ICESat-2 data file + # -- find matching files (for granule, release, version, track) + colnames, collastmod, colerror = icesat2_toolkit.utilities.nsidc_list( + PATH, build=False, timeout=TIMEOUT, parser=parser, pattern=R1, sort=True + ) print(colnames) - #-- print if file was not found + # -- print if file was not found if not colnames: print(colerror) continue - #-- add to lists - for colname,remote_mtime in zip(colnames,collastmod): - #-- save original file to list (expands if getting auxiliary files) + # -- add to lists + for colname, remote_mtime in zip(colnames, collastmod): + # -- save original file to list (expands if getting auxiliary files) original_files.append(input_file) - #-- remote and local versions of the file - remote_files.append(posixpath.join(remote_dir,colname)) - local_files.append(os.path.join(local_dir,colname)) + # -- remote and local versions of the file + remote_files.append(posixpath.join(remote_dir, colname)) + local_files.append(os.path.join(local_dir, colname)) remote_mtimes.append(remote_mtime) remote_names.append(colname) - return original_files, remote_files, remote_names #product_directory, sd, + return original_files, remote_files, remote_names # product_directory, sd, + def json_load(name, path, verbose=False): import json import os - full_name= (os.path.join(path,name+ '.json')) - with open(full_name, 'r') as ifile: - data=json.load(ifile) + full_name = os.path.join(path, name + ".json") + + with open(full_name, "r") as ifile: + data = json.load(ifile) if verbose: - print('loaded from: ',full_name) + print("loaded from: ", full_name) return data -def ATL03_download(username,password, dpath, product_directory, sd, file_name): + +def ATL03_download(username, password, dpath, product_directory, sd, file_name): """ inputs: username: username for https://urs.earthdata.nasa.gov @@ -273,101 +335,119 @@ def ATL03_download(username,password, dpath, product_directory, sd, file_name): """ import icesat2_toolkit.utilities from icesat2_toolkit.read_ICESat2_ATL03 import read_HDF5_ATL03 - HOST = ['https://n5eil01u.ecs.nsidc.org','ATLAS',product_directory,sd, file_name] + + HOST = ["https://n5eil01u.ecs.nsidc.org", "ATLAS", product_directory, sd, file_name] # HOST = ['https://n5eil01u.ecs.nsidc.org','ATLAS','ATL03.003','2018.10.14', # 'ATL03_20181014000347_02350101_003_01.h5'] - print('download to:', dpath+'/'+HOST[-1]) - buffer,error=icesat2_toolkit.utilities.from_nsidc(HOST,username=username, - password=password,local=dpath+'/'+HOST[-1],verbose=True) - #-- raise exception if download error + print("download to:", dpath + "/" + HOST[-1]) + buffer, error = icesat2_toolkit.utilities.from_nsidc( + HOST, + username=username, + password=password, + local=dpath + "/" + HOST[-1], + verbose=True, + ) + # -- raise exception if download error if not buffer: raise Exception(error) -def save_pandas_table(table_dict, name , save_path): + +def save_pandas_table(table_dict, name, save_path): import os + if not os.path.exists(save_path): os.makedirs(save_path) import warnings from pandas import HDFStore from pandas.io.pytables import PerformanceWarning - warnings.filterwarnings('ignore',category=PerformanceWarning) - with HDFStore(save_path+'/'+name+'.h5') as store: - for name,table in table_dict.items(): - store[name]=table + warnings.filterwarnings("ignore", category=PerformanceWarning) -def load_pandas_table_dict(name , save_path): + with HDFStore(save_path + "/" + name + ".h5") as store: + for name, table in table_dict.items(): + store[name] = table + + +def load_pandas_table_dict(name, save_path): import warnings from pandas import HDFStore from pandas.io.pytables import PerformanceWarning - warnings.filterwarnings('ignore',category=PerformanceWarning) - return_dict=dict() - with HDFStore(save_path+'/'+name+'.h5') as store: - #print(store) - #print(store.keys()) + warnings.filterwarnings("ignore", category=PerformanceWarning) + + return_dict = dict() + with HDFStore(save_path + "/" + name + ".h5") as store: + # print(store) + # print(store.keys()) for k in store.keys(): - return_dict[k[1:]]=store.get(k) + return_dict[k[1:]] = store.get(k) return return_dict -def get_beam_hdf_store(ATL03_k): +def get_beam_hdf_store(ATL03_k): import pandas as pd - DD = pd.DataFrame()#columns = ATL03.keys()) + + DD = pd.DataFrame() # columns = ATL03.keys()) for ikey in ATL03_k.keys(): DD[ikey] = ATL03_k[ikey] return DD -def get_beam_var_hdf_store(ATL03_k, ikey): +def get_beam_var_hdf_store(ATL03_k, ikey): import pandas as pd - DD = pd.DataFrame()#columns = ATL03.keys()) + + DD = pd.DataFrame() # columns = ATL03.keys()) DD[ikey] = ATL03_k[ikey] return DD -def write_track_to_HDF5(data_dict, name, path, verbose=False, mode='w'): +def write_track_to_HDF5(data_dict, name, path, verbose=False, mode="w"): import os import h5py - mode = 'w' if mode is None else mode + + mode = "w" if mode is None else mode if not os.path.exists(path): os.makedirs(path) - full_name= (os.path.join(path,name+ '.h5')) + full_name = os.path.join(path, name + ".h5") store = h5py.File(full_name, mode) for k in data_dict.keys(): store1 = store.create_group(k) for kk, I in list(data_dict[k].items()): - store1[kk]=I - #store1.close() + store1[kk] = I + # store1.close() store.close() if verbose: - print('saved at: ' +full_name) + print("saved at: " + full_name) def get_time_for_track(delta_time, atlas_epoch): "returns pandas dataframe" import pandas as pd import convert_GPS_time as cGPS + # Conversion of delta_time to a calendar date temp = cGPS.convert_GPS_time(atlas_epoch[0] + delta_time, OFFSET=0.0) - year = temp['year'][:].astype('int') - month = temp['month'][:].astype('int') - day = temp['day'][:].astype('int') - hour = temp['hour'][:].astype('int') - minute = temp['minute'][:].astype('int') - second = temp['second'][:].astype('int') + year = temp["year"][:].astype("int") + month = temp["month"][:].astype("int") + day = temp["day"][:].astype("int") + hour = temp["hour"][:].astype("int") + minute = temp["minute"][:].astype("int") + second = temp["second"][:].astype("int") + + return pd.DataFrame( + {"year": year, "month": month, "day": day, "hour": hour, "second": second} + ) - return pd.DataFrame({'year':year, 'month':month, 'day':day, 'hour':hour, 'second':second}) -def getATL03_beam(fileT, numpy=False, beam='gt1l', maxElev=1e6): +def getATL03_beam(fileT, numpy=False, beam="gt1l", maxElev=1e6): """ returns 'beam' from fileT as pandas table. fillT path of file @@ -379,38 +459,38 @@ def getATL03_beam(fileT, numpy=False, beam='gt1l', maxElev=1e6): import h5py import convert_GPS_time as cGPS import pandas as pd + # Open the file - ATL03 = h5py.File(fileT, 'r') - lons = ATL03[beam+'/heights/lon_ph'][:] - lats = ATL03[beam+'/heights/lat_ph'][:] + ATL03 = h5py.File(fileT, "r") + lons = ATL03[beam + "/heights/lon_ph"][:] + lats = ATL03[beam + "/heights/lat_ph"][:] # Along track distance from equator i think. - along_track_distance=ATL03[beam+'/heights/dist_ph_along'][:] - across_track_distance=ATL03[beam+'/heights/dist_ph_across'][:] - #dem_h = ATL03[beam+'/geophys_corr/dem_h'][:] - #delta_time_dem_h = ATL03[beam+'/geophys_corr/delta_time'][:] - segment_dist_x=ATL03[beam+'/geolocation/segment_dist_x'][:] - segment_length=ATL03[beam+'/geolocation/segment_length'][:] - segment_id = ATL03[beam+'/geolocation/segment_id'][:] - - delta_time_geolocation = ATL03[beam+'/geolocation/delta_time'][:] - reference_photon_index= ATL03[beam+'/geolocation/reference_photon_index'][:] - ph_index_beg= ATL03[beam+'/geolocation/ph_index_beg'][:] - - - ph_id_count = ATL03[beam+'/heights/ph_id_count'][:] + along_track_distance = ATL03[beam + "/heights/dist_ph_along"][:] + across_track_distance = ATL03[beam + "/heights/dist_ph_across"][:] + # dem_h = ATL03[beam+'/geophys_corr/dem_h'][:] + # delta_time_dem_h = ATL03[beam+'/geophys_corr/delta_time'][:] + segment_dist_x = ATL03[beam + "/geolocation/segment_dist_x"][:] + segment_length = ATL03[beam + "/geolocation/segment_length"][:] + segment_id = ATL03[beam + "/geolocation/segment_id"][:] + + delta_time_geolocation = ATL03[beam + "/geolocation/delta_time"][:] + reference_photon_index = ATL03[beam + "/geolocation/reference_photon_index"][:] + ph_index_beg = ATL03[beam + "/geolocation/ph_index_beg"][:] + + ph_id_count = ATL03[beam + "/heights/ph_id_count"][:] # Nathan says it's the number of seconds since the GPS epoch on midnight Jan. 6, 1980 - delta_time = ATL03[beam+'/heights/delta_time'][:] - #podppd_flag=ATL03[beam+'/geolocation/podppd_flag'][:] + delta_time = ATL03[beam + "/heights/delta_time"][:] + # podppd_flag=ATL03[beam+'/geolocation/podppd_flag'][:] # #Add this value to delta time parameters to compute full gps_seconds - atlas_epoch = ATL03['/ancillary_data/atlas_sdp_gps_epoch'][:] + atlas_epoch = ATL03["/ancillary_data/atlas_sdp_gps_epoch"][:] # Conversion of delta_time to a calendar date temp = cGPS.convert_GPS_time(atlas_epoch[0] + delta_time, OFFSET=0.0) # Express delta_time relative to start time of granule - #delta_time_granule=delta_time-delta_time[0] + # delta_time_granule=delta_time-delta_time[0] # year = temp['year'][:].astype('int') # month = temp['month'][:].astype('int') @@ -422,33 +502,41 @@ def getATL03_beam(fileT, numpy=False, beam='gt1l', maxElev=1e6): # Primary variables of interest # Photon height - heights=ATL03[beam+'/heights/h_ph'][:] - #print(heights.shape) + heights = ATL03[beam + "/heights/h_ph"][:] + # print(heights.shape) # Flag for signal confidence # column index: 0=Land; 1=Ocean; 2=SeaIce; 3=LandIce; 4=InlandWater # values: - #-- -1: Events not associated with a specific surface type - #-- 0: noise - #-- 1: buffer but algorithm classifies as background - #-- 2: low - #-- 3: medium - #-- 4: high - - mask_ocean = ATL03[beam+'/heights/signal_conf_ph'][:, 1] > 2 # ocean points medium or high quality - mask_seaice = ATL03[beam+'/heights/signal_conf_ph'][:, 2] > 2 # sea ice points medium or high quality - mask_total = (mask_seaice | mask_ocean) - - if sum(~mask_total) == (ATL03[beam+'/heights/signal_conf_ph'][:, 1]).size: - print('zero photons, lower photon quality to 2 or higher') - mask_ocean = ATL03[beam+'/heights/signal_conf_ph'][:, 1] > 1 # ocean points medium or high quality - mask_seaice = ATL03[beam+'/heights/signal_conf_ph'][:, 2] > 1 # sea ice points medium or high quality - mask_total = (mask_seaice | mask_ocean) - - signal_confidence = ATL03[beam+'/heights/signal_conf_ph'][:, 1:3].max(1) - #print(signal_confidence.shape) - - #return signal_confidence + # -- -1: Events not associated with a specific surface type + # -- 0: noise + # -- 1: buffer but algorithm classifies as background + # -- 2: low + # -- 3: medium + # -- 4: high + + mask_ocean = ( + ATL03[beam + "/heights/signal_conf_ph"][:, 1] > 2 + ) # ocean points medium or high quality + mask_seaice = ( + ATL03[beam + "/heights/signal_conf_ph"][:, 2] > 2 + ) # sea ice points medium or high quality + mask_total = mask_seaice | mask_ocean + + if sum(~mask_total) == (ATL03[beam + "/heights/signal_conf_ph"][:, 1]).size: + print("zero photons, lower photon quality to 2 or higher") + mask_ocean = ( + ATL03[beam + "/heights/signal_conf_ph"][:, 1] > 1 + ) # ocean points medium or high quality + mask_seaice = ( + ATL03[beam + "/heights/signal_conf_ph"][:, 2] > 1 + ) # sea ice points medium or high quality + mask_total = mask_seaice | mask_ocean + + signal_confidence = ATL03[beam + "/heights/signal_conf_ph"][:, 1:3].max(1) + # print(signal_confidence.shape) + + # return signal_confidence # Add photon rate and background rate to the reader here ATL03.close() @@ -458,26 +546,45 @@ def getATL03_beam(fileT, numpy=False, beam='gt1l', maxElev=1e6): return along_track_dist, elev else: - dF = pd.DataFrame({'heights':heights, 'lons':lons, 'lats':lats, 'signal_confidence':signal_confidence, 'mask_seaice':mask_seaice, - 'delta_time':delta_time, 'along_track_distance':along_track_distance, #'delta_time_granule':delta_time_granule, - 'across_track_distance':across_track_distance,'ph_id_count':ph_id_count})#, - #'year':year, 'month':month, 'day':day, 'hour':hour,'minute':minute , 'second':second}) - - dF_seg = pd.DataFrame({'delta_time':delta_time_geolocation, 'segment_dist_x':segment_dist_x, 'segment_length':segment_length, 'segment_id':segment_id, - 'reference_photon_index':reference_photon_index, 'ph_index_beg':ph_index_beg}) + dF = pd.DataFrame( + { + "heights": heights, + "lons": lons, + "lats": lats, + "signal_confidence": signal_confidence, + "mask_seaice": mask_seaice, + "delta_time": delta_time, + "along_track_distance": along_track_distance, #'delta_time_granule':delta_time_granule, + "across_track_distance": across_track_distance, + "ph_id_count": ph_id_count, + } + ) # , + #'year':year, 'month':month, 'day':day, 'hour':hour,'minute':minute , 'second':second}) + + dF_seg = pd.DataFrame( + { + "delta_time": delta_time_geolocation, + "segment_dist_x": segment_dist_x, + "segment_length": segment_length, + "segment_id": segment_id, + "reference_photon_index": reference_photon_index, + "ph_index_beg": ph_index_beg, + } + ) # Filter out high elevation values - print('seg_dist shape ', segment_dist_x.shape) - print('df shape ',dF.shape) + print("seg_dist shape ", segment_dist_x.shape) + print("df shape ", dF.shape) dF = dF[mask_total] - #dF_seg = dF_seg[mask_total] - #print('df[mask] shape ',dF.shape) + # dF_seg = dF_seg[mask_total] + # print('df[mask] shape ',dF.shape) # Reset row indexing - #dF=dF#.reset_index(drop=True) + # dF=dF#.reset_index(drop=True) return dF, dF_seg -def getATL03_height_correction(fileT, beam='gt1r'): + +def getATL03_height_correction(fileT, beam="gt1r"): """ This method returns relevant data for wave estimates from ALT 07 tracks. returns: Pandas data frame @@ -486,25 +593,27 @@ def getATL03_height_correction(fileT, beam='gt1r'): import h5py import pandas as pd + # Open the file - ATL03 = h5py.File(fileT, 'r') + ATL03 = h5py.File(fileT, "r") ### bulk positions and statistics vars_bulk = [ - 'delta_time', # referenc time since equator crossing - 'dem_h', # best giod approxiamtion - ] + "delta_time", # referenc time since equator crossing + "dem_h", # best giod approxiamtion + ] - D_bulk= dict() + D_bulk = dict() for var in vars_bulk: - D_bulk[var] = ATL03[beam+'/geophys_corr/'+var][:] + D_bulk[var] = ATL03[beam + "/geophys_corr/" + var][:] dF_bulk = pd.DataFrame(D_bulk) ATL03.close() return dF_bulk -def getATL07_beam(fileT, beam='gt1r', maxElev=1e6): + +def getATL07_beam(fileT, beam="gt1r", maxElev=1e6): """ This method returns relevant data for wave estimates from ALT 07 tracks. returns: Pandas data frame @@ -513,81 +622,82 @@ def getATL07_beam(fileT, beam='gt1r', maxElev=1e6): import h5py import pandas as pd + # Open the file - ATL07 = h5py.File(fileT, 'r') + ATL07 = h5py.File(fileT, "r") ### bulk positions and statistics vars_bulk = [ - 'longitude', - 'latitude', - 'height_segment_id',# Height segment ID (10 km segments) - 'seg_dist_x' # Along track distance from the equator crossing to the segment center. - ] + "longitude", + "latitude", + "height_segment_id", # Height segment ID (10 km segments) + "seg_dist_x", # Along track distance from the equator crossing to the segment center. + ] - D_bulk= dict() + D_bulk = dict() for var in vars_bulk: - D_bulk[var] = ATL07[beam+'/sea_ice_segments/'+var] + D_bulk[var] = ATL07[beam + "/sea_ice_segments/" + var] dF_bulk = pd.DataFrame(D_bulk) # Nathan says it's the number of seconds since the GPS epoch on midnight Jan. 6, 1980 - delta_time=ATL07[beam+'/sea_ice_segments/delta_time'][:] + delta_time = ATL07[beam + "/sea_ice_segments/delta_time"][:] # #Add this value to delta time parameters to compute full gps_seconds - atlas_epoch=ATL07['/ancillary_data/atlas_sdp_gps_epoch'][:] - dF_time = get_time_for_track(delta_time,atlas_epoch) - dF_time['delta_time'] = delta_time + atlas_epoch = ATL07["/ancillary_data/atlas_sdp_gps_epoch"][:] + dF_time = get_time_for_track(delta_time, atlas_epoch) + dF_time["delta_time"] = delta_time ### Primary variables of interest - vars = [ - 'across_track_distance', #Across track distance of photons averaged over the sea ice height segment. - 'height_segment_asr_calc', #Computed apparent surface reflectance for the sea ice segment. - 'height_segment_confidence',# # Height segment confidence flag - 'height_segment_fit_quality_flag', # Flag Values: ['-1', '1', '2', '3', '4', '5'] - #Flag Meanings: ['invalid', 'best', 'high', 'med', 'low', 'poor'] - 'height_segment_height', # Beam segment height - 'height_segment_length_seg', # Along track length of segment - 'height_segment_ssh_flag', #Flag for potential leads, 0=sea ice, 1 = sea surface - 'height_segment_surface_error_est', #Error estimate of the surface height - 'height_segment_type',# Flag Values: ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9'] - # Flag Meanings: ['cloud_covered', 'other', 'specular_lead_low_w_bkg', 'specular_lead_low', 'specular_lead_high_w_bkg', 'specular_lead_high', 'dark_lead_smooth_w_bkg', 'dark_lead_smooth' - 'height_segment_w_gaussian', # Width of Gaussian fit - 'height_segment_quality', # Height quality flag, 1 for good fit, 0 for bad - ] - #vars = ['beam_fb_height', 'beam_fb_sigma' , 'beam_fb_confidence' , 'beam_fb_quality_flag'] - - D_heights=dict() + "across_track_distance", # Across track distance of photons averaged over the sea ice height segment. + "height_segment_asr_calc", # Computed apparent surface reflectance for the sea ice segment. + "height_segment_confidence", # # Height segment confidence flag + "height_segment_fit_quality_flag", # Flag Values: ['-1', '1', '2', '3', '4', '5'] + # Flag Meanings: ['invalid', 'best', 'high', 'med', 'low', 'poor'] + "height_segment_height", # Beam segment height + "height_segment_length_seg", # Along track length of segment + "height_segment_ssh_flag", # Flag for potential leads, 0=sea ice, 1 = sea surface + "height_segment_surface_error_est", # Error estimate of the surface height + "height_segment_type", # Flag Values: ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9'] + # Flag Meanings: ['cloud_covered', 'other', 'specular_lead_low_w_bkg', 'specular_lead_low', 'specular_lead_high_w_bkg', 'specular_lead_high', 'dark_lead_smooth_w_bkg', 'dark_lead_smooth' + "height_segment_w_gaussian", # Width of Gaussian fit + "height_segment_quality", # Height quality flag, 1 for good fit, 0 for bad + ] + # vars = ['beam_fb_height', 'beam_fb_sigma' , 'beam_fb_confidence' , 'beam_fb_quality_flag'] + + D_heights = dict() for var in vars: - D_heights[var] = ATL07[beam+'/sea_ice_segments/heights/' +var][:] + D_heights[var] = ATL07[beam + "/sea_ice_segments/heights/" + var][:] dF_heights = pd.DataFrame(D_heights) - vars_env = { - 'mss':'geophysical/height_segment_mss', # Mean sea surface height above WGS-84 reference ellipsoid (range: -105 to 87m), based on the DTU13 model. - 't2m':'geophysical/height_segment_t2m',#Temperature at 2m above the displacement height (K) - 'u2m':'geophysical/height_segment_u2m',#Eastward wind at 2m above the displacement height (m/s-1) - 'v2m':'geophysical/height_segment_v2m',#Northward wind at 2m above the displacement height (m/s-1) - 'n_photons_actual':'stats/n_photons_actual', # Number of photons gathered - 'photon_rate':'stats/photon_rate', #photon_rate + "mss": "geophysical/height_segment_mss", # Mean sea surface height above WGS-84 reference ellipsoid (range: -105 to 87m), based on the DTU13 model. + "t2m": "geophysical/height_segment_t2m", # Temperature at 2m above the displacement height (K) + "u2m": "geophysical/height_segment_u2m", # Eastward wind at 2m above the displacement height (m/s-1) + "v2m": "geophysical/height_segment_v2m", # Northward wind at 2m above the displacement height (m/s-1) + "n_photons_actual": "stats/n_photons_actual", # Number of photons gathered + "photon_rate": "stats/photon_rate", # photon_rate } - D_env=dict() - for var,I in vars_env.items(): - D_env[ var] = ATL07[beam+'/sea_ice_segments/' +I][:] + D_env = dict() + for var, I in vars_env.items(): + D_env[var] = ATL07[beam + "/sea_ice_segments/" + I][:] dF_env = pd.DataFrame(D_env) - - #Df = pd.concat({k: pd.DataFrame(v).T for k, v in data.items()}, axis=0) - DF = pd.concat({ 'time': dF_time, 'ref': dF_bulk, 'heights': dF_heights, 'env': dF_env }, axis=1) + # Df = pd.concat({k: pd.DataFrame(v).T for k, v in data.items()}, axis=0) + DF = pd.concat( + {"time": dF_time, "ref": dF_bulk, "heights": dF_heights, "env": dF_env}, axis=1 + ) ATL07.close() # Filter out high elevation values - DF = DF[(DF['heights']['height_segment_height'] 0) | np.isnan(G_beam.ice) + + lats = list(ice_mask.latitude.data) + lats.sort(reverse=True) + + # find 1st latitude that is completely full with sea ice. + ice_lat_pos = next( + ( + i + for i, j in enumerate( + (ice_mask.sum("longitude") == ice_mask.longitude.size).sel( + latitude=lats + ) + ) + if j + ), + None, + ) + # recreate lat mask based on this criteria + lat_mask = lats < lats[ice_lat_pos] + lat_mask = xr.DataArray( + lat_mask.repeat(ice_mask.longitude.size).reshape(ice_mask.shape), + dims=ice_mask.dims, + coords=ice_mask.coords, + ) + lat_mask["latitude"] = lats + + # combine ice mask and new lat mask + ice_mask = ice_mask + lat_mask + + else: + ice_mask = np.isnan(G_beam.ice) + lats = ice_mask.latitude + + # find closed latituyde with with non-nan data + ice_lat_pos = ( + abs( + lats.where(ice_mask.sum("longitude") > 4, np.nan) + - np.array(lat_range).mean() + ) + .argmin() + .data + ) + + # redefine lat-range + lat_range = lats[ice_lat_pos].data - 2, lats[ice_lat_pos].data + 2 + lat_flag2 = (lat_range[0] < lats.data) & (lats.data < lat_range[1]) + + lat_mask = xr.DataArray( + lat_flag2.repeat(ice_mask.longitude.size).reshape(ice_mask.shape), + dims=ice_mask.dims, + coords=ice_mask.coords, + ) + lat_mask["latitude"] = lats + + # plot 1st figure + def draw_range(lon_range, lat_range, *args, **kargs): + plt.plot( + [lon_range[0], lon_range[1], lon_range[1], lon_range[0], lon_range[0]], + [lat_range[0], lat_range[0], lat_range[1], lat_range[1], lat_range[0]], + *args, + **kargs, + ) + + dir_clev = np.arange(0, 380, 20) + f_clev = np.arange(1 / 40, 1 / 5, 0.01) + fvar = ["ice", "dir", "dp", "spr", "fp", "hs"] + fcmap = [ + plt.cm.Blues_r, + color_schemes.circle_medium_triple, + color_schemes.circle_medium_triple, + plt.cm.Blues, + plt.cm.Blues, + plt.cm.Blues, + ] + fpos = [0, 1, 2, 3, 4, 5] + clevs = [ + np.arange(0, 1, 0.2), + dir_clev, + dir_clev, + np.arange(0, 90, 10), + f_clev, + np.arange(0.5, 9, 0.5), + ] + + font_for_print() + + F = M.figure_axis_xy(4, 3.5, view_scale=0.9, container=True) + plt.suptitle( + track_name_short + " | " + file_name_base[0:-1].replace("_", " "), y=1.3 + ) + lon, lat = G_beam.longitude, G_beam.latitude + + gs = GridSpec(9, 6, wspace=0.1, hspace=0.4) + + for fv, fp, fc, cl in zip(fvar, fpos, fcmap, clevs): + ax1 = F.fig.add_subplot(gs[0:7, fp]) + if fp == 0: + ax1.spines["bottom"].set_visible(False) + ax1.spines["left"].set_visible(False) + ax1.tick_params(labelbottom=True, bottom=True) + + else: + ax1.axis("off") + + plt.plot(G1["lons"], G1["lats"], ".r", markersize=5) + draw_range(lon_range, lat_range_prior, c="red", linewidth=1, zorder=12) + draw_range(lon_range, lat_range, c="blue", linewidth=0.7, zorder=10) + + if fv != "ice": + cm = plt.pcolor(lon, lat, G_beam[fv], vmin=cl[0], vmax=cl[-1], cmap=fc) + if G_beam.ice.shape[0] > 1: + plt.contour(lon, lat, G_beam.ice, colors="black", linewidths=0.6) + else: + cm = plt.pcolor(lon, lat, G_beam[fv], vmin=cl[0], vmax=cl[-1], cmap=fc) + + plt.title(G_beam[fv].long_name.replace(" ", "\n") + "\n" + fv, loc="left") + ax1.axis("equal") + + ax2 = F.fig.add_subplot(gs[-1, fp]) + cbar = plt.colorbar(cm, cax=ax2, orientation="horizontal", aspect=1, fraction=1) + cl_ticks = np.linspace(cl[0], cl[-1], 3) + + cbar.set_ticks(np.round(cl_ticks, 3)) + cbar.set_ticklabels(np.round(cl_ticks, 2)) + + F.save_pup(path=plot_path, name=plot_name + "_hindcast_data") + + G_beam_masked = G_beam.where(~ice_mask, np.nan) + ice_mask_prior = ice_mask.sel(latitude=G_prior.latitude) + G_prior_masked = G_prior.where(~ice_mask_prior, np.nan) + + def test_nan_frac(imask): + "test if False is less then 0.3" + return ((~imask).sum() / imask.size).data < 0.3 + + while test_nan_frac(ice_mask_prior): + print(lat_range_prior) + lat_range_prior = lat_range_prior[0] + 0.5, lat_range_prior[1] + 0.5 + G_prior = sel_data(G_beam, lon_range, lat_range_prior) + ice_mask_prior = ice_mask.sel(latitude=G_prior.latitude) + + G_prior_masked = G_prior.where(~ice_mask_prior, np.nan) + + ### make pandas table with obs track end postitions + + key_list = list(G_prior_masked.keys()) + # define directional and amplitude pairs + # pack as (amp, angle) + key_list_pairs = { + "mean": ("hs", "dir"), + "peak": ("hs", "dp"), + "partion0": ("phs0", "pdp0"), + "partion1": ("phs1", "pdp1"), + "partion2": ("phs2", "pdp2"), + "partion3": ("phs3", "pdp3"), + "partion4": ("phs4", "pdp4"), + } + + key_list_pairs2 = list() + for k in key_list_pairs.values(): + key_list_pairs2.append(k[0]) + key_list_pairs2.append(k[1]) + + key_list_scaler = set(key_list) - set(key_list_pairs2) + + ### derive angle avearge + Tend = pd.DataFrame(index=key_list, columns=["mean", "std", "name"]) + + for k, pair in key_list_pairs.items(): + ave_amp, ave_deg, std_amp, std_deg = waves.get_ave_amp_angle( + G_prior_masked[pair[0]].data, G_prior_masked[pair[1]].data + ) + Tend.loc[pair[0]] = ave_amp, std_amp, G_prior_masked[pair[0]].long_name + Tend.loc[pair[1]] = ave_deg, std_deg, G_prior_masked[pair[1]].long_name + + for k in key_list_scaler: + Tend.loc[k] = ( + G_prior_masked[k].mean().data, + G_prior_masked[k].std().data, + G_prior_masked[k].long_name, + ) + + Tend = Tend.T + Tend["lon"] = [ + ice_mask_prior.longitude.mean().data, + ice_mask_prior.longitude.std().data, + "lontigude", + ] + Tend["lat"] = [ + ice_mask_prior.latitude[ice_mask_prior.sum("longitude") == 0].mean().data, + ice_mask_prior.latitude[ice_mask_prior.sum("longitude") == 0].std().data, + "latitude", + ] + Tend = Tend.T + + Prior = dict() + Prior["incident_angle"] = { + "value": Tend["mean"]["dp"].astype("float"), + "name": Tend["name"]["dp"], + } + Prior["spread"] = { + "value": Tend["mean"]["spr"].astype("float"), + "name": Tend["name"]["spr"], + } + Prior["Hs"] = { + "value": Tend["mean"]["hs"].astype("float"), + "name": Tend["name"]["hs"], + } + Prior["peak_period"] = { + "value": 1 / Tend["mean"]["fp"].astype("float"), + "name": "1/" + Tend["name"]["fp"], + } + + Prior["center_lon"] = { + "value": Tend["mean"]["lon"].astype("float"), + "name": Tend["name"]["lon"], + } + Prior["center_lat"] = { + "value": Tend["mean"]["lat"].astype("float"), + "name": Tend["name"]["lat"], + } + + target_name = "A02_" + track_name + "_hindcast_success" + + MT.save_pandas_table({"priors_hindcast": Tend}, save_name, save_path) +except: + target_name = "A02_" + track_name + "_hindcast_fail" + +def plot_prior(Prior, axx): + angle = Prior["incident_angle"][ + "value" + ] # incident direction in degrees from North clockwise (Meerological convention) + # use + angle_plot = -angle - 90 + axx.quiver( + Prior["center_lon"]["value"], + Prior["center_lat"]["value"], + -np.cos(angle_plot * np.pi / 180), + -np.sin(angle_plot * np.pi / 180), + scale=4.5, + zorder=12, + width=0.1, + headlength=4.5, + minshaft=2, + alpha=0.6, + color="black", + ) + axx.plot( + Prior["center_lon"]["value"], + Prior["center_lat"]["value"], + ".", + markersize=6, + zorder=12, + alpha=1, + color="black", + ) + tstring = ( + " " + + str(np.round(Prior["peak_period"]["value"], 1)) + + "sec \n " + + str(np.round(Prior["Hs"]["value"], 1)) + + "m\n " + + str(np.round(angle, 1)) + + "deg" + ) + plt.text(lon_range[1], Prior["center_lat"]["value"], tstring) + + +try: + # plot 2nd figure + + font_for_print() + F = M.figure_axis_xy(2, 4.5, view_scale=0.9, container=False) + + ax1 = F.ax + lon, lat = G_beam.longitude, G_beam.latitude + ax1.spines["bottom"].set_visible(False) + ax1.spines["left"].set_visible(False) + ax1.tick_params(labelbottom=True, bottom=True) + + plot_prior(Prior, ax1) + + str_list = list() + for i in np.arange(0, 6): + str_list.append( + " " + + str(np.round(Tend.loc["ptp" + str(i)]["mean"], 1)) + + "sec\n " + + str(np.round(Tend.loc["phs" + str(i)]["mean"], 1)) + + "m " + + str(np.round(Tend.loc["pdp" + str(i)]["mean"], 1)) + + "d" + ) + + plt.text(lon_range[1], lat_range[0], "\n ".join(str_list)) + + for vv in zip( + ["pdp0", "pdp1", "pdp2", "pdp3", "pdp4", "pdp5"], + ["phs0", "phs1", "phs3", "phs4", "phs5"], + ): + angle_plot = -Tend.loc[vv[0]]["mean"] - 90 + vsize = (1 / Tend.loc[vv[1]]["mean"]) ** (1 / 2) * 5 + print(vsize) + ax1.quiver( + Prior["center_lon"]["value"], + Prior["center_lat"]["value"], + -np.cos(angle_plot * np.pi / 180), + -np.sin(angle_plot * np.pi / 180), + scale=vsize, + zorder=5, + width=0.1, + headlength=4.5, + minshaft=4, + alpha=0.6, + color="green", + ) + + plt.plot(G1["lons"], G1["lats"], ".r", markersize=5) + draw_range(lon_range, lat_range_prior, c="red", linewidth=1, zorder=11) + draw_range(lon_range, lat_range, c="blue", linewidth=0.7, zorder=10) + + fv = "ice" + if fv != "ice": + plt.pcolor(lon, lat, G_beam[fv].where(~ice_mask, np.nan), cmap=fc) + plt.contour(lon, lat, G_beam.ice, colors="black", linewidths=0.6) + else: + plt.pcolor(lon, lat, G_beam[fv], cmap=fc) + + plt.title( + "Prior\n" + + file_name_base[0:-1].replace("_", " ") + + "\n" + + track_name_short + + "\nIncident angle", + loc="left", + ) + ax1.axis("equal") + + F.save_pup(path=plot_path, name=plot_name + "_hindcast_prior") +except Exception as e: + print(e) + print("print 2nd figure failed") + +MT.json_save( + target_name, save_path, str(datetime.datetime.now().strftime("%Y-%m-%d %H:%M")) +) + +print("done") diff --git a/src/icesat2_tracks/analysis_db/B01_SL_load_single_file.py b/src/icesat2_tracks/analysis_db/B01_SL_load_single_file.py index 37fbfc8d..22c914f9 100644 --- a/src/icesat2_tracks/analysis_db/B01_SL_load_single_file.py +++ b/src/icesat2_tracks/analysis_db/B01_SL_load_single_file.py @@ -1,62 +1,48 @@ -# """ This file open a ICEsat2 tbeam_stats.pyrack applied filters and corections and returns smoothed photon heights on a regular grid in an .nc file. This is python 3.11 """ -import os, sys +import sys +import datetime +import copy +from pathlib import Path + +import xarray as xr +from sliderule import sliderule, icesat2 + from icesat2_tracks.config.IceSAT2_startup import ( mconfig, - xr, color_schemes, font_for_pres, plt, ) - -import geopandas as gpd - -from sliderule import sliderule, icesat2, earthdata - -import shapely -from ipyleaflet import basemaps, Map, GeoData - -import icesat2_tracks.ICEsat2_SI_tools.sliderule_converter_tools as sct -import icesat2_tracks.ICEsat2_SI_tools.io as io -import icesat2_tracks.ICEsat2_SI_tools.beam_stats as beam_stats -import icesat2_tracks.local_modules.m_tools_ph3 as MT -from icesat2_tracks.local_modules import m_general_ph3 as M - - -import h5py, imp, copy -import datetime - +from icesat2_tracks.ICEsat2_SI_tools import ( + sliderule_converter_tools as sct, + io, + beam_stats, +) +from icesat2_tracks.local_modules import m_tools_ph3 as MT, m_general_ph3 as M xr.set_options(display_style="text") - # Select region and retrive batch of tracks - track_name, batch_key, ID_flag = io.init_from_input( sys.argv ) # loads standard experiment -# 20190502052058_05180312_005_01 plot_flag = True hemis = batch_key.split("_")[0] - -save_path = mconfig["paths"]["work"] + "/" + batch_key + "/B01_regrid/" -MT.mkdirs_r(save_path) - -save_path_json = mconfig["paths"]["work"] + "/" + batch_key + "/A01b_ID/" -MT.mkdirs_r(save_path_json) +# Make target directories +basedir = Path(mconfig["paths"]["work"], batch_key) +save_path, save_path_json = Path(basedir, "B01_regrid"), Path(basedir, "A01b_ID") +for p in [save_path, save_path_json]: + MT.mkdirs_r(p) ATL03_track_name = "ATL03_" + track_name + ".h5" # Configure SL Session -sliderule.authenticate("brown", ps_username="mhell", ps_password="Oijaeth9quuh") -icesat2.init( - "slideruleearth.io", organization="brown", desired_nodes=1, time_to_live=90 -) # minutes +icesat2.init("slideruleearth.io") # plot the ground tracks in geographic location @@ -99,10 +85,8 @@ } maximum_height = 30 # (meters) maximum height past dem_h correction -print("STARTS") -gdf = icesat2.atl06p(params_yapc, resources=[ATL03_track_name]) -print("ENDS") -gdf = sct.correct_and_remove_height(gdf, maximum_height) + +gdf = io.get_gdf(ATL03_track_name, params_yapc, maximum_height) cdict = dict() @@ -118,8 +102,6 @@ # main routine for defining the x coordinate and sacing table data - - def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): """ converts a GeoDataFrame from Sliderule to GeoDataFrames for each beam witht the correct columns and names @@ -166,7 +148,7 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): # define reference point and then define 'x' table_data = copy.copy(gdf) -imp.reload(sct) + # the reference point is defined as the most equatorward point of the polygon. # It's distance from the equator is subtracted from the distance of each photon. table_data = sct.define_x_coordinate_from_data(table_data) @@ -181,7 +163,6 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): Ti[kk]["heights_c_weighted_mean"] = Ti[kk]["h_mean"].copy() Ti[kk]["heights_c_std"] = Ti[kk]["h_sigma"].copy() - segment = track_name.split("_")[1][-2:] ID_name = sct.create_ID_name(gdf.iloc[0], segment=segment) print(ID_name) @@ -203,7 +184,6 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): F_atl06.save_light(path=plot_path, name="B01b_ATL06_corrected.png") plt.close() -imp.reload(beam_stats) if plot_flag: font_for_pres() F = M.figure_axis_xy(8, 4.3, view_scale=0.6) diff --git a/src/icesat2_tracks/analysis_db/B02_make_spectra_gFT.py b/src/icesat2_tracks/analysis_db/B02_make_spectra_gFT.py index 28536c01..b23dcfea 100644 --- a/src/icesat2_tracks/analysis_db/B02_make_spectra_gFT.py +++ b/src/icesat2_tracks/analysis_db/B02_make_spectra_gFT.py @@ -1,31 +1,27 @@ -import sys - - """ This file open a ICEsat2 track applied filters and corections and returns smoothed photon heights on a regular grid in an .nc file. This is python 3 """ -from icesat2_tracks.config.IceSAT2_startup import mconfig, xr, plt, np - +import copy +import datetime +import h5py +import time +import sys -from threadpoolctl import threadpool_info, threadpool_limits +import numpy as np +import xarray as xr from pprint import pprint +from scipy.ndimage.measurements import label +from threadpoolctl import threadpool_info, threadpool_limits - -import h5py +import icesat2_tracks.ICEsat2_SI_tools.generalized_FT as gFT import icesat2_tracks.ICEsat2_SI_tools.io as io import icesat2_tracks.ICEsat2_SI_tools.spectral_estimates as spec - -import time -import imp -import copy -from icesat2_tracks.local_modules.m_spectrum_ph3 import spicke_remover -import datetime -import icesat2_tracks.ICEsat2_SI_tools.generalized_FT as gFT -from scipy.ndimage.measurements import label +import icesat2_tracks.local_modules.m_general_ph3 as M +import icesat2_tracks.local_modules.m_spectrum_ph3 as spicke_remover import icesat2_tracks.local_modules.m_tools_ph3 as MT -from icesat2_tracks.local_modules import m_general_ph3 as M +from icesat2_tracks.config.IceSAT2_startup import mconfig, plt import tracemalloc @@ -99,7 +95,7 @@ def linear_gap_fill(F, key_lead, key_int): Ib = Gd[group[1]] ratio = Ia["x"][:].size / Ib["x"][:].size if (ratio > 10) | (ratio < 0.1): - print("bad data ratio ", ratio, 1 / ratio) + # print("bad data ratio ", ratio, 1 / ratio) # TODO: add logger bad_ratio_flag = True if (np.array(nan_fraction).mean() > 0.95) | bad_ratio_flag: @@ -118,7 +114,6 @@ def linear_gap_fill(F, key_lead, key_int): exit() # test LS with an even grid where missing values are set to 0 -imp.reload(spec) print(Gd.keys()) Gi = Gd[list(Gd.keys())[0]] # to select a test beam dist = io.get_beam_var_hdf_store(Gd[list(Gd.keys())[0]], "dist") @@ -153,10 +148,10 @@ def linear_gap_fill(F, key_lead, key_int): print("define global xlims") dist_list = np.array([np.nan, np.nan]) for k in all_beams: - print(k) + # print(k) # TODO: add logger hkey = "heights_c_weighted_mean" x = Gd[k + "/dist"][:] - print(x[0], x[-1]) + # print(x[0], x[-1]) # TODO: add logger dist_list = np.vstack([dist_list, [x[0], x[-1]]]) xlims = np.nanmin(dist_list[:, 0]) - dx, np.nanmin(dist_list[:, 1]) @@ -172,7 +167,7 @@ def linear_gap_fill(F, key_lead, key_int): for k in all_beams: dist_i = io.get_beam_var_hdf_store(Gd[k], "dist") x_mask = (dist_i > xlims[0]) & (dist_i < xlims[1]) - print(k, sum(x_mask["dist"]) / (xlims[1] - xlims[0])) + # print(k, sum(x_mask["dist"]) / (xlims[1] - xlims[0])) # TODO: add logger print("-reduced frequency resolution") @@ -199,7 +194,7 @@ def linear_gap_fill(F, key_lead, key_int): Gi = io.get_beam_hdf_store(Gd[k]) x_mask = (Gi["dist"] > xlims[0]) & (Gi["dist"] < xlims[1]) if sum(x_mask) / (xlims[1] - xlims[0]) < 0.005: - print("------------------- not data in beam found; skip") + # print("------------------- not data in beam found; skip") # TODO: add logger continue Gd_cut = Gi[x_mask] @@ -219,7 +214,7 @@ def linear_gap_fill(F, key_lead, key_int): # compute slope spectra !! dd = np.gradient(dd) - dd, _ = spicke_remover(dd, spreed=10, verbose=False) + dd, _ = spicke_remover.spicke_remover(dd, spreed=10, verbose=False) dd_nans = (np.isnan(dd)) + (Gd_cut["N_photos"] <= 5) # using gappy data diff --git a/src/icesat2_tracks/analysis_db/B03_plot_spectra_ov.py b/src/icesat2_tracks/analysis_db/B03_plot_spectra_ov.py index e69f6728..ca7984d2 100644 --- a/src/icesat2_tracks/analysis_db/B03_plot_spectra_ov.py +++ b/src/icesat2_tracks/analysis_db/B03_plot_spectra_ov.py @@ -3,45 +3,31 @@ This is python 3 """ import sys - -from icesat2_tracks.config.IceSAT2_startup import ( - mconfig, - xr, - color_schemes, - plt, - np, - font_for_print, -) - -import icesat2_tracks.ICEsat2_SI_tools.io as io -import icesat2_tracks.ICEsat2_SI_tools.spectral_estimates as spec - -import time +import numpy as np +import xarray as xr from matplotlib.gridspec import GridSpec +import icesat2_tracks.ICEsat2_SI_tools.io as io import icesat2_tracks.ICEsat2_SI_tools.generalized_FT as gFT -from scipy.ndimage.measurements import label import icesat2_tracks.local_modules.m_tools_ph3 as MT from icesat2_tracks.local_modules import m_general_ph3 as M - +from icesat2_tracks.config.IceSAT2_startup import mconfig, color_schemes, plt, font_for_print track_name, batch_key, test_flag = io.init_from_input( - sys.argv + sys.argv # TODO: Handle via CLI ) # loads standard experiment hemis, batch = batch_key.split("_") load_path = mconfig["paths"]["work"] + batch_key + "/B02_spectra/" -load_file = load_path + "B02_" + track_name # + '.nc' +load_file = load_path + "B02_" + track_name plot_path = ( mconfig["paths"]["plot"] + "/" + hemis + "/" + batch_key + "/" + track_name + "/" -) +) # TODO: Update with pathlib MT.mkdirs_r(plot_path) Gk = xr.open_dataset(load_file + "_gFT_k.nc") Gx = xr.open_dataset(load_file + "_gFT_x.nc") Gfft = xr.open_dataset(load_file + "_FFT.nc") -time.sleep(2) - all_beams = mconfig["beams"]["all_beams"] high_beams = mconfig["beams"]["high_beams"] diff --git a/src/icesat2_tracks/config/IceSAT2_startup.py b/src/icesat2_tracks/config/IceSAT2_startup.py index b14cc5b4..d1e27938 100644 --- a/src/icesat2_tracks/config/IceSAT2_startup.py +++ b/src/icesat2_tracks/config/IceSAT2_startup.py @@ -1,26 +1,11 @@ import os import pathlib -### -# THIS FILE IS A LOCAL FILE -# it is not maintained via git, it contains configs specific to the machine -### +import string -#os.environ["DISPLAY"] = "localhost:10.0" -# 14, 16, work -#standart libraries: -import numpy as np -import matplotlib import matplotlib.pyplot as plt -import matplotlib.colors as colors -import pandas as pd from icesat2_tracks.local_modules import m_colormanager_ph3 as M_color from icesat2_tracks.local_modules import m_tools_ph3 as MT -from icesat2_tracks.local_modules import m_general_ph3 as M - -import string - -import xarray as xr ## Read folders and configuration paths config_dir_path = os.path.dirname(__file__) @@ -37,12 +22,10 @@ #load colorscheme color_schemes=M_color.color(path=mconfig['paths']['config'], name='color_def') - lstrings =iter([i+') ' for i in list(string.ascii_lowercase)]) # define journal fig sizes fig_sizes = mconfig['fig_sizes']['AMS'] - SMALL_SIZE = 8 MEDIUM_SIZE = 10 BIGGER_SIZE = 12 @@ -52,11 +35,11 @@ plt.rc('font', size=SMALL_SIZE, serif='Helvetica Neue', weight='normal') # controls default text sizes plt.rc('text', usetex='false') plt.rc('axes', titlesize=MEDIUM_SIZE, labelweight='normal') # fontsize of the axes title -plt.rc('axes', labelsize=SMALL_SIZE, labelweight='normal') #, family='bold') # fontsize of the x and y labels +plt.rc('axes', labelsize=SMALL_SIZE, labelweight='normal') # fontsize of the x and y labels plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels plt.rc('legend', fontsize=SMALL_SIZE, frameon=False) # legend fontsize -plt.rc('figure', titlesize=MEDIUM_SIZE, titleweight='bold', autolayout=True) #, family='bold') # fontsize of the figure title +plt.rc('figure', titlesize=MEDIUM_SIZE, titleweight='bold', autolayout=True) # fontsize of the figure title plt.rc('path', simplify=True) plt.rcParams['figure.figsize'] = (10, 8) plt.rcParams['pcolor.shading'] = 'auto' @@ -66,15 +49,7 @@ plt.rc('axes', labelsize= MEDIUM_SIZE, labelweight='normal') plt.rc('axes.spines', top= False, right=False ) - - -def font_for_print(): - - SMALL_SIZE = 6 - MEDIUM_SIZE = 8 - BIGGER_SIZE = 10 # not used. CP - legend_properties = {'weight':'bold'} # not used. CP - +def font_for_print(SMALL_SIZE = 6, MEDIUM_SIZE = 8): plt.rc('font', size=SMALL_SIZE, serif='Helvetica Neue', weight='normal') # controls default text sizes plt.rc('text', usetex='false') plt.rc('axes', titlesize=MEDIUM_SIZE, labelweight='normal') # fontsize of the axes title @@ -85,21 +60,13 @@ def font_for_print(): plt.rc('figure', titlesize=MEDIUM_SIZE, titleweight='bold', autolayout=True) #, family='bold') # fontsize of the figure title plt.rc('axes', labelsize= SMALL_SIZE, labelweight='normal') -def font_for_pres(): - - SMALL_SIZE = 10 - MEDIUM_SIZE = 12 - BIGGER_SIZE = 14 # not used. CP - legend_properties = {'weight':'bold'} # not used. CP - - - plt.rc('font', size=SMALL_SIZE, serif='Helvetica Neue', weight='normal') # controls default text sizes +def font_for_pres(SMALL_SIZE = 10, MEDIUM_SIZE = 12): + plt.rc('font', size=SMALL_SIZE, serif='Helvetica Neue', weight='normal') # controls default text sizes plt.rc('text', usetex='false') plt.rc('axes', titlesize=MEDIUM_SIZE, labelweight='normal') # fontsize of the axes title plt.rc('axes', labelsize=SMALL_SIZE, labelweight='normal') #, family='bold') # fontsize of the x and y labels plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels plt.rc('legend', fontsize=SMALL_SIZE, frameon=False) # legend fontsize - plt.rc('figure', titlesize=MEDIUM_SIZE, titleweight='bold', autolayout=True) #, family='bold') # fontsize of the figure title - + plt.rc('figure', titlesize=MEDIUM_SIZE, titleweight='bold', autolayout=True) # fontsize of the figure title plt.rc('axes', labelsize= SMALL_SIZE, labelweight='normal')