diff --git a/src/icesat2_tracks/analysis_db/A02c_IOWAGA_thredds_prior.py b/src/icesat2_tracks/analysis_db/A02c_IOWAGA_thredds_prior.py index 4ca09890..bc24735f 100644 --- a/src/icesat2_tracks/analysis_db/A02c_IOWAGA_thredds_prior.py +++ b/src/icesat2_tracks/analysis_db/A02c_IOWAGA_thredds_prior.py @@ -649,4 +649,5 @@ def plot_prior(Prior, axx): make_iowaga_threads_prior_app = makeapp(run_A02c_IOWAGA_thredds_prior, name="threads-prior") if __name__ == "__main__": + logging.basicConfig(level=logging.DEBUG) make_iowaga_threads_prior_app() 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 1cad5357..68a58712 100644 --- a/src/icesat2_tracks/analysis_db/B02_make_spectra_gFT.py +++ b/src/icesat2_tracks/analysis_db/B02_make_spectra_gFT.py @@ -110,9 +110,7 @@ def run_B02_make_spectra_gFT( nan_fraction = list() for k in all_beams: 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] - ) + nan_fraction.append(np.sum(np.isnan(heights_c_std)) / heights_c_std.shape[0]) del heights_c_std @@ -575,4 +573,5 @@ def make_dummy_beam(GG, beam): make_spectra_app = makeapp(run_B02_make_spectra_gFT, name="makespectra") if __name__ == "__main__": + logging.basicConfig(level=logging.DEBUG) make_spectra_app() 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 f0d65327..e5816c5a 100644 --- a/src/icesat2_tracks/analysis_db/B03_plot_spectra_ov.py +++ b/src/icesat2_tracks/analysis_db/B03_plot_spectra_ov.py @@ -34,6 +34,7 @@ _logger = logging.getLogger(__name__) + def plot_wavenumber_spectrogram(ax, Gi, clev, title=None, plot_photon_density=True): if Gi.k[0] == 0: Gi = Gi.sel(k=Gi.k[1:]) @@ -138,9 +139,7 @@ def run_B03_plot_spectra_ov( for k in all_beams: I = Gk.sel(beam=k) I2 = Gx.sel(beam=k) - plt.plot( - I["lon"], I["lat"], ".", c=col_dict[k], markersize=0.7, linewidth=0.3 - ) + plt.plot(I["lon"], I["lat"], ".", c=col_dict[k], markersize=0.7, linewidth=0.3) plt.plot(I2["lon"], I2["lat"], "|", c=col_dict[k], markersize=0.7) plt.xlabel("lon") @@ -184,8 +183,7 @@ def run_B03_plot_spectra_ov( # TODO: refactor to make more readable. CP G_gFT_wmean = ( - Gk["gFT_PSD_data"].where(~np.isnan(Gk["gFT_PSD_data"]), 0) - * Gk["N_per_stancil"] + Gk["gFT_PSD_data"].where(~np.isnan(Gk["gFT_PSD_data"]), 0) * Gk["N_per_stancil"] ).sum("beam") / Gk["N_per_stancil"].sum("beam") G_gFT_wmean["N_per_stancil"] = Gk["N_per_stancil"].sum("beam") @@ -218,9 +216,7 @@ def run_B03_plot_spectra_ov( ) dd = 10 * np.log10(Gplot) dd = dd.where(~np.isinf(dd), np.nan) - clev_log = ( - M.clevels([dd.quantile(0.01).data, dd.quantile(0.98).data * 1.2], 31) * 1 - ) + clev_log = M.clevels([dd.quantile(0.01).data, dd.quantile(0.98).data * 1.2], 31) * 1 xlims = Gmean.x[0] / 1e3, Gmean.x[-1] / 1e3 @@ -312,9 +308,7 @@ def run_B03_plot_spectra_ov( I = Gfft.sel(beam=k) plt.scatter( (x0 + I.x.data) / 1e3, - I.power_spec.sel(k=slice(k_max_range[0], k_max_range[2])).integrate( - "k" - ), + I.power_spec.sel(k=slice(k_max_range[0], k_max_range[2])).integrate("k"), s=0.5, marker=".", c="blue", @@ -375,19 +369,14 @@ def run_B03_plot_spectra_ov( .argmax() .data ) - xpp = x_pos_sel[ - [int(i) for i in np.round(np.linspace(0, x_pos_sel.size - 1, 4))] - ] + xpp = x_pos_sel[[int(i) for i in np.round(np.linspace(0, x_pos_sel.size - 1, 4))]] xpp = np.insert(xpp, 0, x_pos_max) for i in xpp: F = M.figure_axis_xy(6, 8, container=True, view_scale=0.8) plt.suptitle( - "gFT Model and Spectrograms | x=" - + str(Gk.x[i].data) - + " \n" - + track_name, + "gFT Model and Spectrograms | x=" + str(Gk.x[i].data) + " \n" + track_name, y=0.95, ) gs = GridSpec(5, 6, wspace=0.2, hspace=0.7) @@ -456,9 +445,7 @@ def run_B03_plot_spectra_ov( plt.ylabel("relative slope (m/m)") # TODO: compute xlabel as fstring. CP plt.xlabel( - "segment distance $\eta$ (km) @ x=" - + fltostr(Gx_1.x.data / 1e3, 2) - + "km" + "segment distance $\eta$ (km) @ x=" + fltostr(Gx_1.x.data / 1e3, 2) + "km" ) # spectra @@ -484,9 +471,7 @@ def run_B03_plot_spectra_ov( dd = Gk_1.gFT_PSD_data plt.plot(Gk_1.k, dd, color="gray", linewidth=0.5, alpha=0.5) - dd = Gk_1.gFT_PSD_data.rolling( - k=10, min_periods=1, center=True - ).mean() + dd = Gk_1.gFT_PSD_data.rolling(k=10, min_periods=1, center=True).mean() plt.plot(Gk_1.k, dd, color=col_d[k], linewidth=0.8) # handle the 'All-NaN slice encountered' warning if np.all(np.isnan(dd.data)): @@ -566,14 +551,10 @@ def run_B03_plot_spectra_ov( plt.ylabel("relative slope (m/m)") # TODO: compute xlabel as fstring. CP plt.xlabel( - "segment distance $\eta$ (km) @ x=" - + fltostr(Gx_1.x.data / 1e3, 2) - + "km" + "segment distance $\eta$ (km) @ x=" + fltostr(Gx_1.x.data / 1e3, 2) + "km" ) - F.save_pup( - path=str(plot_path / "B03_spectra"), name=f"B03_freq_reconst_x{i}" - ) + F.save_pup(path=str(plot_path / "B03_spectra"), name=f"B03_freq_reconst_x{i}") MT.json_save( "B03_success", @@ -587,4 +568,5 @@ def run_B03_plot_spectra_ov( plot_spectra = makeapp(run_B03_plot_spectra_ov, name="plotspectra") if __name__ == "__main__": + logging.basicConfig(level=logging.DEBUG) plot_spectra() diff --git a/src/icesat2_tracks/analysis_db/B04_angle.py b/src/icesat2_tracks/analysis_db/B04_angle.py index 838a8b7b..4066c557 100644 --- a/src/icesat2_tracks/analysis_db/B04_angle.py +++ b/src/icesat2_tracks/analysis_db/B04_angle.py @@ -50,6 +50,7 @@ _logger = logging.getLogger(__name__) + def run_B04_angle( track_name: str = Option(..., callback=validate_track_name_steps_gt_1), batch_key: str = Option(..., callback=validate_batch_key), @@ -141,12 +142,10 @@ def run_B04_angle( #### Define Prior Pperiod = Prior.loc[["ptp0", "ptp1", "ptp2", "ptp3", "ptp4", "ptp5"]]["mean"] - Pdir = Prior.loc[["pdp0", "pdp1", "pdp2", "pdp3", "pdp4", "pdp5"]][ - "mean" - ].astype("float") - Pspread = Prior.loc[["pspr0", "pspr1", "pspr2", "pspr3", "pspr4", "pspr5"]][ - "mean" - ] + Pdir = Prior.loc[["pdp0", "pdp1", "pdp2", "pdp3", "pdp4", "pdp5"]]["mean"].astype( + "float" + ) + Pspread = Prior.loc[["pspr0", "pspr1", "pspr2", "pspr3", "pspr4", "pspr5"]]["mean"] Pperiod = Pperiod[~np.isnan(list(Pspread))] Pdir = Pdir[~np.isnan(list(Pspread))] @@ -265,10 +264,11 @@ def run_B04_angle( prior_angle = Prior_smth.Prior_direction * 180 / np.pi if (abs(prior_angle) > 80).all(): _logger.critical( - "Prior angle is " + - prior_angle.min().data + " " + - prior_angle.max().data + - ". quit.", + "Prior angle is " + + prior_angle.min().data + + " " + + prior_angle.max().data + + ". quit.", ) dd_save = { "time": time.asctime(time.localtime(time.time())), @@ -454,9 +454,7 @@ def plot_instance( "-", c=next(col_list), ) - plt.xlim( - x_concat[y_concat == y_pos][0], x_concat[y_concat == y_pos][-1] - ) + plt.xlim(x_concat[y_concat == y_pos][0], x_concat[y_concat == y_pos][-1]) plt.xlabel("meter") F.ax3 = F.fig.add_subplot(gs[2:, 0:-1]) @@ -467,9 +465,7 @@ def plot_instance( ) if optimze is True: - SM.plot_optimze( - color="r", markersize=10, zorder=12, label="Dual Annealing" - ) + SM.plot_optimze(color="r", markersize=10, zorder=12, label="Dual Annealing") if sample is True: SM.plot_sample( @@ -626,9 +622,7 @@ def plot_instance( amp_Z = 1 prior_sel = { "alpha": ( - Prior_smth.sel( - k=k_prime_max, method="nearest" - ).Prior_direction.data, + Prior_smth.sel(k=k_prime_max, method="nearest").Prior_direction.data, Prior_smth.sel(k=k_prime_max, method="nearest").Prior_spread.data, ) } @@ -659,9 +653,7 @@ def get_instance(k_pair): Prior_smth.sel( k=k_prime_max, method="nearest" ).Prior_direction.data, - Prior_smth.sel( - k=k_prime_max, method="nearest" - ).Prior_spread.data, + Prior_smth.sel(k=k_prime_max, method="nearest").Prior_spread.data, ) } @@ -678,9 +670,7 @@ def get_instance(k_pair): min=k_prime_max * 0.5, max=k_prime_max * 1.5, ) - SM.params.add( - "K_amp", amp_Z, vary=False, min=amp_Z * 0.0, max=amp_Z * 5 - ) + SM.params.add("K_amp", amp_Z, vary=False, min=amp_Z * 0.0, max=amp_Z * 5) L_sample_i = None L_optimize_i = None @@ -716,9 +706,7 @@ def get_instance(k_pair): "alpha", alpha_dx, burn=N_sample_chain_burn, plot_flag=False ) fitter = SM.fitter # MCMC results - z_model = SM.objective_func( - fitter.params, *fitting_args, test_flag=True - ) + z_model = SM.objective_func(fitter.params, *fitting_args, test_flag=True) cost = (fitter.residual**2).sum() / (z_concat**2).sum() if plot_flag: @@ -824,9 +812,7 @@ def get_instance(k_pair): cost_stack = dict() marginal_stack = dict() L_sample = pd.DataFrame(index=["alpha", "group_phase", "K_prime", "K_amp"]) - L_optimize = pd.DataFrame( - index=["alpha", "group_phase", "K_prime", "K_amp"] - ) + L_optimize = pd.DataFrame(index=["alpha", "group_phase", "K_prime", "K_amp"]) L_brute = pd.DataFrame(index=["alpha", "group_phase", "K_prime", "K_amp"]) for kk, I in A.items(): @@ -1016,4 +1002,5 @@ def get_instance(k_pair): make_b04_angle_app = makeapp(run_B04_angle, name="B04_angle") if __name__ == "__main__": + logging.basicConfig(level=logging.DEBUG) make_b04_angle_app() diff --git a/src/icesat2_tracks/analysis_db/B05_define_angle.py b/src/icesat2_tracks/analysis_db/B05_define_angle.py index c214a50f..7c1060c4 100644 --- a/src/icesat2_tracks/analysis_db/B05_define_angle.py +++ b/src/icesat2_tracks/analysis_db/B05_define_angle.py @@ -41,6 +41,7 @@ _logger = logging.getLogger(__name__) + def derive_weights(weights): """ Normalize weights to have a minimum of 0 @@ -324,9 +325,7 @@ def define_angle( ax_list[group] = ax0 data = corrected_marginals.isel(x=xi).sel(beam_group=group) - weights = derive_weights( - Marginals.weight.isel(x=xi).sel(beam_group=group) - ) + weights = derive_weights(Marginals.weight.isel(x=xi).sel(beam_group=group)) weights = weights**2 # derive angle axis @@ -348,9 +347,9 @@ def define_angle( data_collect[group] = data_wmean data_collect = xr.concat(data_collect.values(), dim="beam_group") - final_data = (group_weight * data_collect).sum( + final_data = (group_weight * data_collect).sum("beam_group") / group_weight.sum( "beam_group" - ) / group_weight.sum("beam_group").data + ).data plt.sca(ax_sum) plt.stairs(final_data, x_angle, color="k", alpha=1, linewidth=0.8) @@ -377,9 +376,7 @@ def define_angle( priors_k = Marginals.Prior_direction[~np.isnan(k_mask.isel(x=xi))] for pk in priors_k: - ax_final.axvline( - pk, color=color_schemes.cascade2, linewidth=1, alpha=0.7 - ) + ax_final.axvline(pk, color=color_schemes.cascade2, linewidth=1, alpha=0.7) plt.stairs(final_data, x_angle, color="k", alpha=0.5, linewidth=0.8) @@ -465,9 +462,7 @@ def define_angle( _title = f"{track_name}\nPower Spectra (m/m)$^2$ k$^{{-1}}$" plt.title(_title, loc="left") - cbar = plt.colorbar( - fraction=0.018, pad=0.01, orientation="vertical", label="Power" - ) + cbar = plt.colorbar(fraction=0.018, pad=0.01, orientation="vertical", label="Power") cbar.outline.set_visible(False) clev_ticks = np.round(clev_spec[::3], 0) cbar.set_ticks(clev_ticks) @@ -528,9 +523,9 @@ def define_angle( i_spec = weighted_spec.sel(x=slice(x_range[0], x_range[-1])) i_dir = corrected_marginals.sel(x=slice(x_range[0], x_range[-1])) - dir_data = (i_dir * i_dir.N_data).sum( + dir_data = (i_dir * i_dir.N_data).sum(["beam_group", "x"]) / i_dir.N_data.sum( ["beam_group", "x"] - ) / i_dir.N_data.sum(["beam_group", "x"]) + ) lims = get_first_and_last_nonzero_data(dir_data) plot_data = build_plot_data(dir_data, i_spec, lims) @@ -570,4 +565,5 @@ def define_angle( define_angle_app = makeapp(define_angle, name="B04_angle") if __name__ == "__main__": + logging.basicConfig(level=logging.DEBUG) define_angle_app() diff --git a/src/icesat2_tracks/analysis_db/B06_correct_separate_var.py b/src/icesat2_tracks/analysis_db/B06_correct_separate_var.py index 4c8e688f..7a24c415 100644 --- a/src/icesat2_tracks/analysis_db/B06_correct_separate_var.py +++ b/src/icesat2_tracks/analysis_db/B06_correct_separate_var.py @@ -42,6 +42,7 @@ _logger = logging.getLogger(__name__) + def get_correct_breakpoint(pw_results): br_points = [i for i in pw_results.keys() if "breakpoint" in i] @@ -99,7 +100,6 @@ def get_correct_breakpoint(pw_results): def get_breakingpoints(xx, dd): - convergence_flag = True n_breakpoints = 3 while convergence_flag: @@ -124,7 +124,6 @@ def get_breakingpoints(xx, dd): def define_noise_wavenumber_piecewise(data_xr, plot_flag=False): - data_log = np.log(data_xr) k = data_log.k.data @@ -183,7 +182,6 @@ def tanh_fitler(x, x_cutoff, sigma_g=0.01): def get_k_x_corrected(Gk, theta=0, theta_flag=False): - if not theta_flag: return np.nan, np.nan @@ -248,14 +246,15 @@ def save_table(data, tablename, save_path): io.save_pandas_table(data, tablename, save_path) except Exception as e: tabletoremove = save_path + tablename + ".h5" - _logger.warning(f"Failed to save table with error {e}. " - f"Removing {tabletoremove} and re-trying..") + _logger.warning( + f"Failed to save table with error {e}. " + f"Removing {tabletoremove} and re-trying.." + ) os.remove(tabletoremove) io.save_pandas_table(data, tablename, save_path) def buil_G_error(Gk_sel, PSD_list, list_name): - return xr.DataArray( data=PSD_list, coords=Gk_sel.drop("N_per_stancil").coords, name=list_name ).expand_dims("beam") @@ -267,7 +266,6 @@ def run_B06_correct_separate_var( ID_flag: bool = True, output_dir: str = typer.Option(..., callback=validate_output_dir), ): - color_schemes.colormaps2(31, gamma=1) col_dict = color_schemes.rels @@ -307,9 +305,7 @@ def run_B06_correct_separate_var( file_suffixes = ["_gFT_k.nc", "_gFT_x.nc"] Gk, Gx = [ - xr.open_dataset( - load_path_work / "B02_spectra" / f"B02_{track_name}{suffix}" - ) + xr.open_dataset(load_path_work / "B02_spectra" / f"B02_{track_name}{suffix}") for suffix in file_suffixes ] @@ -512,9 +508,7 @@ def run_B06_correct_separate_var( + " to " + str(np.round(Gx.isel(x=-1).lat.mean().data, 2)) ) - plt.title( - next(fn) + "Mean Displacement Spectra\n(lat=" + lat_str + ")", loc="left" - ) + plt.title(next(fn) + "Mean Displacement Spectra\n(lat=" + lat_str + ")", loc="left") dd = 10 * np.log((G_gFT_smth / G_gFT_smth.k).isel(x=slice(0, -1))) dd = dd.where(~np.isinf(dd), np.nan) @@ -680,17 +674,13 @@ def reconstruct_displacement(Gx_1, Gk_1, T3, k_thresh): T3["heights_c_model_err"] = np.interp( T3["dist"], continous_x_grid, concented_err ) - T3["heights_c_residual"] = ( - T3["heights_c_weighted_mean"] - T3["heights_c_model"] - ) + T3["heights_c_residual"] = T3["heights_c_weighted_mean"] - T3["heights_c_model"] B3_v2[bb] = T3 Gx_v2[bb] = Gx_k try: - G_angle = xr.open_dataset( - load_path_angle / (f"B05_{track_name}_angle_pdf.nc") - ) + G_angle = xr.open_dataset(load_path_angle / (f"B05_{track_name}_angle_pdf.nc")) except ValueError as e: _logger.warning(f"{e} no angle data found, skip angle corretion") @@ -776,4 +766,5 @@ def reconstruct_displacement(Gx_1, Gk_1, T3, k_thresh): correct_separate_app = makeapp(run_B06_correct_separate_var, name="correct-separate") if __name__ == "__main__": + logging.basicConfig(level=logging.DEBUG) correct_separate_app()