From 65cc3da775ab1c9bbba42edfe86ec3266b0dd26e Mon Sep 17 00:00:00 2001 From: Camilo Diaz Date: Tue, 16 Jan 2024 00:15:47 -0500 Subject: [PATCH 1/5] step5 not working here --- src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py | 2 +- src/icesat2_tracks/analysis_db/B04_angle.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py b/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py index 32280284..1f916da1 100644 --- a/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py +++ b/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py @@ -1002,4 +1002,4 @@ def create_weight(self, freq=None, plot_flag=True, flim=None, max_nfev=None): if self.freq_cut_flag and freq is None: result = np.insert(result, 0, 0) - return result + return result \ No newline at end of file diff --git a/src/icesat2_tracks/analysis_db/B04_angle.py b/src/icesat2_tracks/analysis_db/B04_angle.py index 473bf1a4..54b9dbac 100644 --- a/src/icesat2_tracks/analysis_db/B04_angle.py +++ b/src/icesat2_tracks/analysis_db/B04_angle.py @@ -784,7 +784,8 @@ def get_instance(k_pair): try: LL = pd.concat(L_collect) MT.save_pandas_table({"L_sample": LL}, save_name + "_res_table", save_path) -except: +except Exception as e: + print(e) pass # plot From f48f904517c433029774c95a16b756285a45c74e Mon Sep 17 00:00:00 2001 From: Camilo Diaz Date: Tue, 16 Jan 2024 03:40:20 -0500 Subject: [PATCH 2/5] updated make_spectra and generalize_FT to fix step2 and step5 integration --- .../ICEsat2_SI_tools/generalized_FT.py | 183 ++++++++++-------- .../analysis_db/B02_make_spectra_gFT.py | 121 ++++++------ 2 files changed, 166 insertions(+), 138 deletions(-) diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py b/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py index 1f916da1..3e88f22a 100644 --- a/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py +++ b/src/icesat2_tracks/ICEsat2_SI_tools/generalized_FT.py @@ -2,6 +2,7 @@ import icesat2_tracks.ICEsat2_SI_tools.spectral_estimates as spec import icesat2_tracks.ICEsat2_SI_tools.lanczos as lanczos +import matplotlib.pyplot as plt def rebin(data, dk): @@ -14,7 +15,6 @@ def rebin(data, dk): Gmean["k_bins"] = k_low[0:-1] Gmean = Gmean.rename({"k_bins": "k"}) return Gmean, k_low_limits - # define weight function @@ -29,9 +29,13 @@ def smooth_data_to_weight(dd, m=150): dd_fake[2 * m : -2 * m] = dd weight = lanczos.lanczos_filter_1d_wrapping(np.arange(dd_fake.size), dd_fake, m) + weight = weight[2 * m : -2 * m] weight = weight / weight.max() + # ensure non-negative weights + weight[weight < 0.02] = 0.02 + return weight @@ -65,7 +69,7 @@ def get_weights_from_data( pars = Spec_fft.set_parameters(flim=np.sqrt(9.81 * k[-1]) / 2 / np.pi) k_max = (pars["f_max"].value * 2 * np.pi) ** 2 / 9.81 - if method == "gaussian": + if method is "gaussian": # simple gaussian weight def gaus(x, x_0, amp, sigma_g): return amp * np.exp(-0.5 * ((x - x_0) / sigma_g) ** 2) @@ -73,7 +77,7 @@ def gaus(x, x_0, amp, sigma_g): weight = gaus(k, k_max, 1, 0.02) ** (1 / 2) params = None - elif method == "parametric": + elif method is "parametric": # JONSWAP weight f = np.sqrt(9.81 * k) / (2 * np.pi) weight = Spec_fft.create_weight(freq=f, plot_flag=False, max_nfev=max_nfev) @@ -82,14 +86,11 @@ def gaus(x, x_0, amp, sigma_g): Spec_fft.fitter.params.pretty_print() params = Spec_fft.fitter.params - + else: raise ValueError(" 'method' must be either 'gaussian' or 'parametric' ") if plot_flag: - import matplotlib.pyplot as plt - - plt.plot( k_fft[1:], Spec_fft.data, c="gray", label="FFT for Prior", linewidth=0.5 ) @@ -97,7 +98,6 @@ def gaus(x, x_0, amp, sigma_g): k, weight, zorder=12, c="black", label="Fitted model to FFT", linewidth=0.5 ) plt.xlim(k[0], k[-1]) - # add pemnalty floor weight = weight + weight.max() * 0.1 @@ -108,6 +108,18 @@ def gaus(x, x_0, amp, sigma_g): return weight, params +def define_weight_shutter(weight, k, Ncut=3): + "creates masking function to lower high wavenumber weights, Ncut is the numnber by which the spectral peak is multiplied" + # Limit high wavenumbers weights + weight_shutter = weight * 0 + 1 + ki_cut = weight.argmax() * Ncut # k of peak + N_res = weight.size - ki_cut + if N_res < 1: + return weight_shutter + weight_shutter[ki_cut:] = np.linspace(1, 0, N_res) + return weight_shutter + + def make_xarray_from_dict(D, name, dims, coords): import xarray as xr @@ -124,26 +136,26 @@ def define_weights(stancil, prior, x, y, dx, k, max_nfev, plot_flag=False): return weights normalized to 1, prior_pars used for the next iteration """ - if (type(prior[0]) is bool) and not prior[ - 0 - ]: # prior = (False, None), this is the first iteration + if (type(prior[0]) is bool) and not prior[0]: # fit function to data weight, prior_pars = get_weights_from_data( x, y, dx, stancil, k, max_nfev, plot_flag=plot_flag, method="parametric" ) + weight_name = "$P_{init}$ from FFT" - elif ( - type(prior) is tuple - ): + elif type(prior) is tuple: + # combine old and new weights weight = 0.2 * smooth_data_to_weight(prior[0]) + 0.8 * prior[1] - weight_name = "smth. $P_{i-1}$" prior_pars = {"alpha": None, "amp": None, "f_max": None, "gamma": None} - else: + else: # prior = weight, this is all other iterations weight = smooth_data_to_weight(prior) weight_name = "smth. from data" prior_pars = {"alpha": None, "amp": None, "f_max": None, "gamma": None} + # Limit high wavenumbers weights + weight = weight * define_weight_shutter(weight, k, Ncut=3) + if plot_flag: import matplotlib.pyplot as plt @@ -176,9 +188,7 @@ def __init__(self, x, data, L, dx, wavenumber, data_error=None, ov=None): """ self.Lmeters = L - self.ov = ( - int(L / 2) if ov is None else ov - ) # when not defined in create_chunk_boundaries then L/2 + self.ov = int(L / 2) if ov is None else ov self.x = x self.dx = dx @@ -228,9 +238,9 @@ def cal_spectrogram( Lmeters, dk = self.Lmeters, self.dk Lpoints = self.Lpoints Lpoints_full = int(Lmeters / self.dx) + self.xlims = (np.round(X.min()), X.max()) if xlims is None else xlims - # init Lomb scargle object with noise as nummy data () def calc_gFT_apply(stancil, prior): """ windows the data accoding to stencil and applies LS spectrogram @@ -243,11 +253,11 @@ def calc_gFT_apply(stancil, prior): ta = time.perf_counter() x_mask = (stancil[0] <= X) & (X <= stancil[-1]) - print(stancil[1]) x = X[x_mask] if ( - x.size / Lpoints < 0.1 + x.size / Lpoints < 0.40 ): # if there are not enough photos set results to nan + print(" -- data density to low, skip stancil") return { "stancil_center": stancil[1], "p_hat": np.concatenate([self.k * np.nan, self.k * np.nan]), @@ -265,7 +275,6 @@ def calc_gFT_apply(stancil, prior): y_var = y.var() FT = generalized_Fourier(x, y, self.k) - if plot_flag: import matplotlib.pyplot as plt @@ -282,14 +291,39 @@ def calc_gFT_apply(stancil, prior): # define error err = ERR[x_mask] if ERR is not None else 1 - print("weights : ", time.perf_counter() - ta) + print("compute time weights : ", time.perf_counter() - ta) + ta = time.perf_counter() - FT.define_problem(weight, err) # 1st arg is Penalty, 2nd is error + FT.define_problem(weight, err) # solve problem: p_hat = FT.solve() - print("solve : ", time.perf_counter() - ta) + if np.isnan(np.mean(p_hat)): + print(" -- inversion nan!") + print(" -- data fraction", x.size / Lpoints) + print( + " -- weights:", + np.mean(weight), + "err:", + np.mean(err), + "y:", + np.mean(y), + ) + print(" -- skip stancil") + return { + "stancil_center": stancil[1], + "p_hat": np.concatenate([self.k * np.nan, self.k * np.nan]), + "inverse_stats": np.nan, + "y_model_grid": np.nan, + "y_data_grid": np.nan, + "x_size": x.size, + "PSD": False, + "weight": False, + "spec_adjust": np.nan, + } + + print("compute time solve : ", time.perf_counter() - ta) ta = time.perf_counter() x_pos = (np.round((x - stancil[0]) / self.dx, 0)).astype("int") @@ -301,7 +335,7 @@ def calc_gFT_apply(stancil, prior): y_data_grid = np.copy(eta) * np.nan y_data_grid[x_pos] = y - inverse_stats = FT.get_stats(self.dk, Lpoints_full, print_flag=True) + inverse_stats = FT.get_stats(self.dk, Lpoints_full, print_flag=plot_flag) # add fitting parameters of Prior to stats dict for k, I in prior_pars.items(): try: @@ -309,19 +343,17 @@ def calc_gFT_apply(stancil, prior): except: inverse_stats[k] = np.nan - print("stats : ", time.perf_counter() - ta) - - # multiply with the standard deviation of the data to get dimensions right - PSD = power_from_model( - p_hat, dk, self.k.size, x.size, Lpoints - ) + print("compute time stats : ", time.perf_counter() - ta) - if self.k.size * 2 > x.size: - col = "red" - else: - col = "blue" + # multiply with the standard deviation of the data to get dimensions right + PSD = power_from_model(p_hat, dk, self.k.size, x.size, Lpoints) if plot_flag: + if self.k.size * 2 > x.size: + col = "red" + else: + col = "blue" + plt.plot(self.k, PSD, color=col, label="GFT fit", linewidth=0.5) plt.title( "non-dim Spectral Segment Models, 2M=" @@ -357,28 +389,34 @@ def calc_gFT_apply(stancil, prior): return return_dict - # % derive L2 stancil - self.stancil_iter = spec.create_chunk_boundaries_unit_lengths( - Lmeters, self.xlims, ov=self.ov, iter_flag=True + # derive L2 stancil + self.stancil_iter_list = spec.create_chunk_boundaries_unit_lengths( + Lmeters, self.xlims, ov=self.ov, iter_flag=False ) - + self.stancil_iter = iter(self.stancil_iter_list.T.tolist()) + # apply func to all stancils Spec_returns = list() + # form: PSD_from_GFT, weight_used in inversion prior = False, False + N_stencil = len(self.stancil_iter_list.T) + Ni = 1 for ss in copy.copy(self.stancil_iter): + print(Ni, "/", N_stencil, "Stancils") + # prior step if prior[0] is False: # make NL fit of piors do not exist - print("1st step with NL-fit") + print("1st step: with NL-fit") I_return = calc_gFT_apply(ss, prior=prior) - prior = I_return["PSD"], I_return["weight"] - + prior = I_return["PSD"], I_return["weight"] # 2nd step if prior[0] is False: - print("priors still false skip 2nd step") + print("1st GFD failed (priors[0]=false), skip 2nd step") else: - print("2nd step use set priors:", type(prior[0]), type(prior[0])) + print("2nd step: use set priors:", type(prior[0]), type(prior[1])) + print(prior[0][0:3], prior[1][0:3]) I_return = calc_gFT_apply(ss, prior=prior) - prior = I_return["PSD"], I_return["weight"] + prior = I_return["PSD"], I_return["weight"] Spec_returns.append( dict( @@ -395,7 +433,9 @@ def calc_gFT_apply(stancil, prior): ) ) ) - + + Ni += 1 + # unpack resutls of the mapping: GFT_model = dict() Z_model = dict() @@ -441,28 +481,26 @@ def calc_gFT_apply(stancil, prior): self.N_per_stancil = N_per_stancil chunk_positions = np.array(list(D_specs.keys())) - self.N_stancils = len(chunk_positions) # number of spectral realizations + self.N_stancils = len(chunk_positions) # number of spectral realizatiobs # repack data, create xarray # 1st LS spectal estimates - G_power_data = make_xarray_from_dict( D_specs, "gFT_PSD_data", ["k"], {"k": self.k} ) - G_power_data = xr.concat(G_power_data.values(), dim="x").T # .to_dataset() + G_power_data = xr.concat(G_power_data.values(), dim="x").T G_power_model = make_xarray_from_dict( D_specs_model, "gFT_PSD_model", ["k"], {"k": self.k} ) - G_power_model = xr.concat(G_power_model.values(), dim="x").T # .to_dataset() + G_power_model = xr.concat(G_power_model.values(), dim="x").T self.G = G_power_model self.G.name = "gFT_PSD_model" - # 2nd FFT(Y_model) G_model_Z = make_xarray_from_dict(Z_model, "Z_hat", ["k"], {"k": self.k}) - G_model_Z = xr.concat(G_model_Z.values(), dim="x").T + G_model_Z = xr.concat(G_model_Z.values(), dim="x").T # 3rd GFT_model_coeff_A = dict() @@ -475,24 +513,19 @@ def calc_gFT_apply(stancil, prior): I[1], dims=["k"], coords={"k": self.k, "x": xi}, name="gFT_sin_coeff" ) - GFT_model_coeff_A = xr.concat( - GFT_model_coeff_A.values(), dim="x" - ).T - GFT_model_coeff_B = xr.concat( - GFT_model_coeff_B.values(), dim="x" - ).T + GFT_model_coeff_A = xr.concat(GFT_model_coeff_A.values(), dim="x").T + GFT_model_coeff_B = xr.concat(GFT_model_coeff_B.values(), dim="x").T # add weights to the data weights_k = make_xarray_from_dict(weights, "weight", ["k"], {"k": self.k}) - weights_k = xr.concat(weights_k.values(), dim="x").T + weights_k = xr.concat(weights_k.values(), dim="x").T # .to_dataset() - # 4th: model in real space eta = np.arange(0, self.Lmeters + self.dx, self.dx) - self.Lmeters / 2 y_model_eta = make_xarray_from_dict(y_model, "y_model", ["eta"], {"eta": eta}) y_data_eta = make_xarray_from_dict(y_data, "y_data", ["eta"], {"eta": eta}) - y_model_eta = xr.concat(y_model_eta.values(), dim="x").T - y_data_eta = xr.concat(y_data_eta.values(), dim="x").T + y_model_eta = xr.concat(y_model_eta.values(), dim="x").T + y_data_eta = xr.concat(y_data_eta.values(), dim="x").T # merge wavenumber datasets self.GG = xr.merge( @@ -517,7 +550,7 @@ def calc_gFT_apply(stancil, prior): for k, I in Pars.items(): if I is not np.nan: PP2[k] = I - + keys = Pars[next(iter(PP2))].keys() keys_DF = list(set(keys) - set(["model_error_k", "model_error_x"])) params_dataframe = pd.DataFrame(index=keys_DF) @@ -542,19 +575,18 @@ def calc_gFT_apply(stancil, prior): ) sta, ste = xi - self.Lmeters / 2, xi + self.Lmeters / 2 + x_pos = (np.round((X[(sta <= X) & (X <= ste)] - sta) / self.dx)).astype( "int" ) x_err = np.copy(eta) * np.nan - # check sizes and adjust if necessary. + # check sizes and adjust if necessary. if x_pos.size > I["model_error_x"].size: x_pos = x_pos[0 : I["model_error_x"].size] print("adjust x") elif x_pos.size < I["model_error_x"].size: - I["model_error_x"] = I["model_error_x"][ - 0:-1 - ] + I["model_error_x"] = I["model_error_x"][0:-1] print("adjust y") x_err[x_pos] = I["model_error_x"] @@ -605,6 +637,7 @@ def parceval(self, add_attrs=True, weight_data=False): import copy DATA = self.data + L = self.Lmeters X = self.x def get_stancil_var_apply(stancil): @@ -712,7 +745,7 @@ def power_from_model(p_hat, dk, M, N_x, N_x_full): """ Z = complex_represenation(p_hat, M, N_x_full) - spec, _ = Z_to_power_gFT(Z, dk, N_x, N_x_full) + spec, _ = Z_to_power_gFT(Z, dk, N_x, N_x_full) # use spec_incomplete # spectral density respesenting the incomplete data return spec @@ -783,7 +816,6 @@ def solve(self): return p_hat - def model(self): "returns the model dimensional units" if "p_hat" not in self.__dict__: @@ -831,8 +863,9 @@ def parceval(self, dk, Nx_full): return pars def get_stats(self, dk, Nx_full, print_flag=False): - residual = self.ydata - self.model() + + Lmeters = self.x[-1] - self.x[0] pars = { "data_var": self.ydata.var(), "model_var": self.model().var(), @@ -863,8 +896,6 @@ def get_stats(self, dk, Nx_full, print_flag=False): class get_prior_spec: def __init__(self, freq, data): - """ """ - import numpy as np import lmfit as LM self.LM = LM @@ -969,7 +1000,7 @@ def runningmean(self, var, m, tailcopy=False): print("0 Dimension is smaller then averaging length") return rr = np.asarray(var) * np.nan - + var_range = np.arange(m, int(s[0]) - m - 1, 1) for i in var_range[np.isfinite(var[m : int(s[0]) - m - 1])]: rr[int(i)] = np.nanmean(var[i - m : i + m]) @@ -1002,4 +1033,4 @@ def create_weight(self, freq=None, plot_flag=True, flim=None, max_nfev=None): if self.freq_cut_flag and freq is None: result = np.insert(result, 0, 0) - return result \ No newline at end of file + return result 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 b23dcfea..c2e3e1d1 100644 --- a/src/icesat2_tracks/analysis_db/B02_make_spectra_gFT.py +++ b/src/icesat2_tracks/analysis_db/B02_make_spectra_gFT.py @@ -1,32 +1,40 @@ +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 """ -import copy -import datetime -import h5py -import time -import sys +from icesat2_tracks.config.IceSAT2_startup import ( + mconfig, + plt, +) +from threadpoolctl import threadpool_info, threadpool_limits +from pprint import pprint 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 icesat2_tracks.ICEsat2_SI_tools.generalized_FT as gFT +import h5py import icesat2_tracks.ICEsat2_SI_tools.io as io import icesat2_tracks.ICEsat2_SI_tools.spectral_estimates as spec -import icesat2_tracks.local_modules.m_general_ph3 as M -import icesat2_tracks.local_modules.m_spectrum_ph3 as spicke_remover + +import time +import imp +import copy +import icesat2_tracks.ICEsat2_SI_tools.spicke_remover as 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_tools_ph3 as MT -from icesat2_tracks.config.IceSAT2_startup import mconfig, plt +from icesat2_tracks.local_modules import m_general_ph3 as M + import tracemalloc -def linear_gap_fill(F, key_lead, key_int): +def linear_gap_fill(F,key_lead, key_int): """ F pd.DataFrame key_lead key in F that determined the independent coordindate @@ -43,9 +51,7 @@ def linear_gap_fill(F, key_lead, key_int): track_name, batch_key, test_flag = io.init_from_input( sys.argv ) # loads standard experiment - hemis, batch = batch_key.split("_") - ATlevel = "ATL03" load_path = mconfig["paths"]["work"] + "/" + batch_key + "/B01_regrid/" @@ -62,28 +68,26 @@ def linear_gap_fill(F, key_lead, key_int): + batch_key + "/" + track_name - + "/B_spectra/" + + "/B03_spectra/" ) MT.mkdirs_r(plot_path) MT.mkdirs_r(save_path) bad_track_path = mconfig["paths"]["work"] + "bad_tracks/" + batch_key + "/" - all_beams = mconfig["beams"]["all_beams"] high_beams = mconfig["beams"]["high_beams"] low_beams = mconfig["beams"]["low_beams"] N_process = 4 +print("N_process=", N_process) - -# open with hdf5 Gd = h5py.File(load_path + "/" + track_name + "_B01_binned.h5", "r") # test amount of nans in the data nan_fraction = list() for k in all_beams: - heights_c_std = io.get_beam_var_hdf_store(Gd[k], "dist") + heights_c_std = io.get_beam_var_hdf_store(Gd[k], "x") nan_fraction.append(np.sum(np.isnan(heights_c_std)) / heights_c_std.shape[0]) del heights_c_std @@ -95,7 +99,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) # TODO: add logger + print("bad data ratio ", ratio, 1 / ratio) bad_ratio_flag = True if (np.array(nan_fraction).mean() > 0.95) | bad_ratio_flag: @@ -114,10 +118,12 @@ 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") +dist = io.get_beam_var_hdf_store(Gd[list(Gd.keys())[0]], "x") +# make dataframe form hdf5 # derive spectal limits # Longest deserved period: T_max = 40 # sec @@ -129,6 +135,7 @@ def linear_gap_fill(F, key_lead, key_int): Lpoints = int(np.round(min_datapoint) * 10) Lmeters = Lpoints * dx + print("L number of gridpoint:", Lpoints) print("L length in km:", Lmeters / 1e3) print("approx number windows", 2 * dist.iloc[-1] / Lmeters - 1) @@ -137,6 +144,7 @@ def linear_gap_fill(F, key_lead, key_int): lambda_min = 9.81 * T_min**2 / (2 * np.pi) flim = 1 / T_min + oversample = 2 dlambda = Lmeters * oversample dk = 2 * np.pi / dlambda @@ -148,31 +156,23 @@ 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) # TODO: add logger - hkey = "heights_c_weighted_mean" - x = Gd[k + "/dist"][:] - # print(x[0], x[-1]) # TODO: add logger + print(k) + hkey = "h_mean" + x = Gd[k + "/x"][:] + print(x[0], x[-1]) dist_list = np.vstack([dist_list, [x[0], x[-1]]]) xlims = np.nanmin(dist_list[:, 0]) - dx, np.nanmin(dist_list[:, 1]) -dist_lim = 2000e3 # maximum distanc in the sea ice tha tis analysed: - - -if (xlims[1] - xlims[0]) > dist_lim: - xlims = xlims[0], xlims[0] + dist_lim - print("-reduced xlims length to ", xlims[0] + dist_lim, "m") - for k in all_beams: - dist_i = io.get_beam_var_hdf_store(Gd[k], "dist") + dist_i = io.get_beam_var_hdf_store(Gd[k], "x") x_mask = (dist_i > xlims[0]) & (dist_i < xlims[1]) - # print(k, sum(x_mask["dist"]) / (xlims[1] - xlims[0])) # TODO: add logger + print(k, sum(x_mask["x"]) / (xlims[1] - xlims[0])) print("-reduced frequency resolution") kk = kk[::2] - print("set xlims: ", xlims) print( "Loop start: ", @@ -180,37 +180,37 @@ def linear_gap_fill(F, key_lead, key_int): tracemalloc.get_traced_memory()[1] / 1e6, ) - G_gFT = dict() G_gFT_x = dict() G_rar_fft = dict() Pars_optm = dict() -k = all_beams[0] + +k = all_beams[1] for k in all_beams: tracemalloc.start() # ------------------------------- use gridded data - hkey = "heights_c_weighted_mean" + + # sliderule version + hkey = "h_mean" + hkey_sigma = "h_sigma" + Gi = io.get_beam_hdf_store(Gd[k]) - x_mask = (Gi["dist"] > xlims[0]) & (Gi["dist"] < xlims[1]) + x_mask = (Gi["x"] > xlims[0]) & (Gi["x"] < xlims[1]) if sum(x_mask) / (xlims[1] - xlims[0]) < 0.005: - # print("------------------- not data in beam found; skip") # TODO: add logger - continue + print("------------------- not data in beam found; skip") Gd_cut = Gi[x_mask] - x = Gd_cut["dist"] + x = Gd_cut["x"] del Gi # cut data: x_mask = (x >= xlims[0]) & (x <= xlims[1]) x = x[x_mask] - dd = np.copy(Gd_cut[hkey]) - dd_error = np.copy(Gd_cut["heights_c_std"]) + dd_error = np.copy(Gd_cut[hkey_sigma]) + dd_error[np.isnan(dd_error)] = 100 - # plt.hist(1/dd_weight, bins=40) - F = M.figure_axis_xy(6, 3) - plt.subplot(2, 1, 1) # compute slope spectra !! dd = np.gradient(dd) @@ -222,12 +222,7 @@ def linear_gap_fill(F, key_lead, key_int): x_no_nans = x[~dd_nans] dd_error_no_nans = dd_error[~dd_nans] - plt.plot( - x_no_nans, dd_no_nans, ".", color="black", markersize=1, label="slope (m/m)" - ) - plt.legend() - plt.show() - + print("gFT") with threadpool_limits(limits=N_process, user_api="blas"): pprint(threadpool_info()) @@ -330,6 +325,7 @@ def linear_gap_fill(F, key_lead, key_int): plt.ylim(ylims[0], ylims[-1]) plt.show() + # add x-mean spectal error estimate to xarray S.parceval(add_attrs=True, weight_data=False) # assign beam coordinate @@ -340,10 +336,10 @@ def linear_gap_fill(F, key_lead, key_int): GG.coords["spec_adjust"] = (("x", "beam"), np.expand_dims(GG["spec_adjust"], 1)) # add more coodindates to the Dataset - x_coord_no_gaps = linear_gap_fill(Gd_cut, "dist", "x") - y_coord_no_gaps = linear_gap_fill(Gd_cut, "dist", "y") + x_coord_no_gaps = linear_gap_fill(Gd_cut, "x", "x") + y_coord_no_gaps = linear_gap_fill(Gd_cut, "x", "y") mapped_coords = spec.sub_sample_coords( - Gd_cut["dist"], x_coord_no_gaps, y_coord_no_gaps, S.stancil_iter, map_func=None + Gd_cut["x"], x_coord_no_gaps, y_coord_no_gaps, S.stancil_iter, map_func=None ) GG.coords["x_coord"] = GG_x.coords["x_coord"] = ( @@ -362,10 +358,10 @@ def linear_gap_fill(F, key_lead, key_int): GG.coords["x_coord"][nan_mask] = np.nan GG.coords["y_coord"][nan_mask] = np.nan - lons_no_gaps = linear_gap_fill(Gd_cut, "dist", "lons") - lats_no_gaps = linear_gap_fill(Gd_cut, "dist", "lats") + lons_no_gaps = linear_gap_fill(Gd_cut, "x", "lons") + lats_no_gaps = linear_gap_fill(Gd_cut, "x", "lats") mapped_coords = spec.sub_sample_coords( - Gd_cut["dist"], lons_no_gaps, lats_no_gaps, S.stancil_iter, map_func=None + Gd_cut["x"], lons_no_gaps, lats_no_gaps, S.stancil_iter, map_func=None ) GG.coords["lon"] = GG_x.coords["lon"] = ( @@ -445,10 +441,11 @@ def get_stancil_nans(stancil): label="mean FFT", ) plt.legend() - plt.show() + except: pass time.sleep(3) + plt.close("all") del Gd_cut @@ -458,7 +455,7 @@ def get_stancil_nans(stancil): MT.save_pandas_table(Pars_optm, save_name + "_params", save_path) -# repack data +# repack data def repack_attributes(DD): attr_dim_list = list(DD.keys()) for k in attr_dim_list: From e78b79d1f56588908602f46b3fff68ae29101a64 Mon Sep 17 00:00:00 2001 From: Camilo Diaz Date: Tue, 16 Jan 2024 03:42:08 -0500 Subject: [PATCH 3/5] commented step5 test in the workflow --- .github/workflows/test-B01_SL_load_single_file.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/test-B01_SL_load_single_file.yml b/.github/workflows/test-B01_SL_load_single_file.yml index 9f070a7a..fa227d1d 100644 --- a/.github/workflows/test-B01_SL_load_single_file.yml +++ b/.github/workflows/test-B01_SL_load_single_file.yml @@ -29,5 +29,5 @@ jobs: run: python src/icesat2_tracks/analysis_db/B03_plot_spectra_ov.py SH_20190502_05180312 SH_testSLsinglefile2 True - name: fourth step IOWAGA thredds run: python src/icesat2_tracks/analysis_db/A02c_IOWAGA_thredds_prior.py SH_20190502_05180312 SH_testSLsinglefile2 True - - name: Fifth step B04_angle - run: python src/icesat2_tracks/analysis_db/B04_angle.py SH_20190502_05180312 SH_testSLsinglefile2 True + # - name: Fifth step B04_angle + # run: python src/icesat2_tracks/analysis_db/B04_angle.py SH_20190502_05180312 SH_testSLsinglefile2 True From 53c9b9be91e6f5fc8c24b0305a71e158686f69ba Mon Sep 17 00:00:00 2001 From: Camilo Diaz Date: Thu, 18 Jan 2024 13:27:15 -0500 Subject: [PATCH 4/5] Apply suggestions from code review Co-authored-by: Carlos Paniagua --- src/icesat2_tracks/analysis_db/B04_angle.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/icesat2_tracks/analysis_db/B04_angle.py b/src/icesat2_tracks/analysis_db/B04_angle.py index 54b9dbac..11449b53 100644 --- a/src/icesat2_tracks/analysis_db/B04_angle.py +++ b/src/icesat2_tracks/analysis_db/B04_angle.py @@ -785,8 +785,7 @@ def get_instance(k_pair): LL = pd.concat(L_collect) MT.save_pandas_table({"L_sample": LL}, save_name + "_res_table", save_path) except Exception as e: - print(e) - pass + print(f"This is a warning: {e}) # plot font_for_print() From 9f0ebf5e94fc7140a8b80b7d1737095ba7899374 Mon Sep 17 00:00:00 2001 From: Carlos Paniagua Date: Thu, 18 Jan 2024 13:32:09 -0500 Subject: [PATCH 5/5] fix: typo --- src/icesat2_tracks/analysis_db/B04_angle.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/icesat2_tracks/analysis_db/B04_angle.py b/src/icesat2_tracks/analysis_db/B04_angle.py index 11449b53..4ad3dc7f 100644 --- a/src/icesat2_tracks/analysis_db/B04_angle.py +++ b/src/icesat2_tracks/analysis_db/B04_angle.py @@ -785,7 +785,7 @@ def get_instance(k_pair): LL = pd.concat(L_collect) MT.save_pandas_table({"L_sample": LL}, save_name + "_res_table", save_path) except Exception as e: - print(f"This is a warning: {e}) + print(f"This is a warning: {e}") # plot font_for_print()