From 3543a53ddf545d5e08eb968732c8ac705f365f14 Mon Sep 17 00:00:00 2001 From: Camilo Diaz Date: Wed, 8 Nov 2023 10:29:36 -0500 Subject: [PATCH 1/7] changes in the B01_SL_load_single_file.py file to run sliderule integration --- .../ICEsat2_SI_tools/beam_stats.py | 249 ++++++++++++++++ .../sliderule_converter_tools.py | 55 ++-- .../analysis_db/B01_SL_load_single_file.py | 270 ++++++++++++++++++ 3 files changed, 555 insertions(+), 19 deletions(-) create mode 100644 src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py create mode 100644 src/icesat2_tracks/analysis_db/B01_SL_load_single_file.py diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py b/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py new file mode 100644 index 00000000..66950e2b --- /dev/null +++ b/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py @@ -0,0 +1,249 @@ +import numpy as np +import pandas as pd +import icesat2_tracks.ICEsat2_SI_tools.spectral_estimates as spec +#import io as io_local +import icesat2_tracks.ICEsat2_SI_tools.io as io_local + +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec + +def derive_beam_statistics(Gd, all_beams, Lmeter=10e3, dx =10): + """ + this method returns a dict of dataframes with the beam statistics + Gd is a dict of beam tables or a hdf5 file + all_beams is a list of all beams + Lemter is the length of the segment in meters for the statistics + dx is the nominal resolution of the ATL06 data in meters + """ + import h5py + + D=dict() + for k in all_beams: + if isinstance(Gd, dict): + Gi = Gd[k] + elif isinstance(Gd, h5py.File): + Gi = io_local.get_beam_hdf_store(Gd[k]) + else: + print('Gd is neither dict nor hdf5 file') + break + + dd = Gi['h_mean'] + xx = Gi['dist'] + + def get_var(sti): + mask = (sti[0] < xx) & (xx <= sti[1]) + return np.nanvar(dd[mask]) + + def get_N(sti): + mask = (sti[0] < xx) & (xx <= sti[1]) + return dd[mask].size + + def get_lat(sti): + mask = (sti[0] < xx) & (xx <= sti[1]) + return np.nanmean(Gi['lats'][mask]) + + iter_x = spec.create_chunk_boundaries_unit_lengths( Lmeter, [ xx.min(), xx.max()],ov =0, iter_flag= False)[1,:] + + stencil_iter = spec.create_chunk_boundaries_unit_lengths( Lmeter, [ xx.min(), xx.max()],ov =0, iter_flag= True) + var_list = np.array(list(map(get_var, stencil_iter))) + + stencil_iter = spec.create_chunk_boundaries_unit_lengths( Lmeter, [ xx.min(), xx.max()],ov =0, iter_flag= True) + N_list = np.array(list(map(get_N, stencil_iter))) + + stencil_iter = spec.create_chunk_boundaries_unit_lengths( Lmeter, [ xx.min(), xx.max()],ov =0, iter_flag= True) + lat_list = np.array(list(map(get_lat, stencil_iter))) + + # make Dataframe + df = pd.DataFrame() + df['x'] = iter_x + df['lat'] = lat_list + df['var'] = var_list + df['N'] = N_list * 2* dx / Lmeter + + D[k] = df + + return D + +def plot_beam_statistics(D, high_beams, low_beams, col_dict, track_name =None): + """ + Plots the beam statistics in a 2 x 2 plot + D is a dict of dataframes with the beam statistics + high_beams is a list of high beams + low_beams is a list of low beams + col_dict is a dict with the colors for the beams + track_name is the name of the track + """ + import matplotlib.pyplot as plt + + if track_name is not None: + plt.suptitle(track_name, fontsize=10) + + import matplotlib.gridspec as gridspec + + gs = gridspec.GridSpec(2, 3) + + # make 2 x 2 plot + ax1 = plt.subplot(gs[0, 0]) + for k in high_beams: + plt.plot(D[k]['x']/1e3, np.sqrt(D[k]['var']), '.', color= col_dict[k], markersize=4, label=k) + + plt.title('high beams std', loc='left') + #plt.xlabel('along track distance (km)') + plt.ylabel('segment std log(m)') + #plt.ylim(0, 20) + ax1.set_yscale('log') + + ax2 = plt.subplot(gs[1, 0]) + for k in high_beams: + Di = D[k]['N'] + Di[Di ==0] =np.nan + plt.plot(D[k]['x']/1e3, D[k]['N'], '.', color= col_dict[k], markersize=4, label=k) + + plt.title('high beams N', loc='left') + plt.xlabel('along track distance (km)') + plt.ylabel('Point Density (m)') + + ax3 = plt.subplot(gs[0, 1]) + for k in low_beams: + plt.plot(D[k]['x']/1e3, np.sqrt(D[k]['var']), '.', color= col_dict[k], markersize=4, label=k) + + plt.title('low beams std', loc='left') + #plt.xlabel('along track distance (km)') + #plt.ylabel('segment std (m)') + #plt.ylim(0, 20) + ax3.set_yscale('log') + + ax4 = plt.subplot(gs[1, 1]) + for k in low_beams: + Di = D[k]['N'] + Di[Di ==0] =np.nan + plt.plot(D[k]['x']/1e3, D[k]['N'], '.', color= col_dict[k], markersize=4, label=k) + + plt.title('low beams N', loc='left') + plt.xlabel('along track distance (km)') + #plt.ylabel('Point density (m)') + + ax5 = plt.subplot(gs[0:2, 2]) + + lat_shift = 0 + for k in low_beams: + Di = D[k] + plt.scatter(Di['x']/1e3, Di['lat']+lat_shift, s= np.exp(Di['N'] *5) , marker='.', color= col_dict[k], label=k, alpha = 0.3) + lat_shift = lat_shift + 2 + + for k in high_beams: + Di = D[k] + plt.scatter(Di['x']/1e3, Di['lat']+lat_shift, s= np.exp(Di['N'] *5) , marker='.', color= col_dict[k], label=k, alpha = 0.3) + lat_shift = lat_shift + 2 + + plt.title('Density in space', loc='left') + plt.ylabel('Latitude (deg)') + plt.xlabel('along track distance (km)') + plt.legend() + plt.show() + + + # # make 2 x 2 plot + # plt.subplot(2, 3, 1) + # for k in high_beams: + # plt.plot(D[k]['x']/1e3, np.sqrt(D[k]['var']), '.', color= col_dict[k], markersize=4, label=k) + + # plt.title('high beams std', loc='left') + # #plt.xlabel('along track distance (km)') + # plt.ylabel('segment std (m)') + # plt.ylim(0, 20) + + # plt.subplot(2, 3, 3) + # for k in high_beams: + # Di = D[k]['N'] + # Di[Di ==0] =np.nan + # plt.plot(D[k]['lat'], D[k]['N'], '.', color= col_dict[k], markersize=4, label=k) + + # plt.title('high beams N', loc='left') + # plt.xlabel('along track distance (km) or Lat') + # plt.ylabel('Point Density (m)') + + + # plt.subplot(2, 3, 2) + # for k in low_beams: + # plt.plot(D[k]['x']/1e3, np.sqrt(D[k]['var']), '.', color= col_dict[k], markersize=4, label=k) + + # plt.title('low beams std', loc='left') + # #plt.xlabel('along track distance (km)') + # #plt.ylabel('segment std (m)') + # plt.ylim(0, 20) + + # plt.subplot(2, 3, 4) + # for k in low_beams: + # Di = D[k]['N'] + # Di[Di ==0] =np.nan + # plt.plot(D[k]['lat'], D[k]['N'], '.', color= col_dict[k], markersize=4, label=k) + + # plt.title('low beams N', loc='left') + # plt.xlabel('along track distance (km) or Lat') + # #plt.ylabel('Point density (m)') + # #return F + + # plt.subplot(2,3, 5) + # plt.scatter(D['gt1l']['x']/1e3, D['gt1l']['lat'], s= np.exp(D['gt1l']['N'] *10) , marker='.', color= col_dict['gt1l'], label='gt1l', alpha = 0.3) + + # plt.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, wspace=0.3, hspace=0.3) + +## plot track stats basics for sliderules ATL06 output + +def plot_ATL06_track_data( G2, cdict): + """ + Plots the beam statistics in a 3 x 3 plot + G2 is a GeoDataFrame from SL (ATL06) + title is the title of the plot + cdict is a dict with the colors for each 'spot' + returns a figure object + """ + gs = gridspec.GridSpec(3, 3) + + ax1 = plt.subplot(gs[0, 0:2]) + ax2 = plt.subplot(gs[1, 0:2]) + ax3 = plt.subplot(gs[2, 0:2]) + ax4 = plt.subplot(gs[0, 2]) + ax5 = plt.subplot(gs[1, 2]) + ax6 = plt.subplot(gs[2, 2]) + + for sp in G2['spot'].unique(): + Gc = G2[G2['spot'] == 1] + + Gc['h_mean_gradient'] = np.gradient(Gc['h_mean']) + ts_config = {'marker': '.', 'markersize': 0.2, 'linestyle': 'none', 'color': cdict[sp], 'alpha': 0.3} + hist_confit = {'density': True, 'color': cdict[sp], 'alpha': 0.3} + + ax1.plot(Gc.geometry.y, Gc['h_mean'], **ts_config) + ax2.plot(Gc.geometry.y, Gc['h_mean_gradient'], **ts_config) + ax3.plot(Gc.geometry.y, Gc['n_fit_photons'], **ts_config) + + Gc['h_mean'].plot.hist(ax=ax4, bins=30, **hist_confit) + Gc['h_mean_gradient'].plot.hist(ax=ax5, bins=np.linspace(-5, 5, 30), **hist_confit) + Gc['rms_misfit'].plot.hist(ax=ax6, bins=30, **hist_confit) + + ax1.set_ylabel('h_mean (m)') + ax2.set_ylabel('slope (m/m)') + ax3.set_ylabel('N Photons') + ax3.set_xlabel('Latitude (degree)') + ax1.set_xticklabels([]) + ax2.set_xticklabels([]) + + ax1.axhline(0, color='k', linestyle='-', linewidth=0.8) + ax2.axhline(0, color='k', linestyle='-', linewidth=0.8) + + ax1.set_title('Height', loc='left') + ax2.set_title('Slope', loc='left') + ax3.set_title('Photons per extend', loc='left') + + ax4.set_title('Histograms', loc='left') + ax5.set_title('Histograms', loc='left') + + ax6.set_title('Error Hist.', loc='left') + ax6.set_xlabel('rms_misfit (m)') + + for axi in [ax4, ax5, ax6]: + axi.set_ylabel('') + + return [ax1, ax2, ax3, ax4, ax5, ax6] \ No newline at end of file diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py b/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py index f543b888..c5af7029 100644 --- a/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py +++ b/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py @@ -1,6 +1,20 @@ from ipyleaflet import basemaps + +# height correction tools +def correct_and_remove_height(Gi, height_limit): + """ + corrects the height and removes the points with height +- height_limit + """ + h_mean_corrected = (Gi['h_mean'] - Gi['dem_h']) + false_height_mask = abs(h_mean_corrected) > height_limit + Gi['h_mean']=h_mean_corrected + Gi.drop(columns=['dem_h'], inplace=True) + return Gi[~false_height_mask] + + + # polygon tools def make_plot_polygon(poly_test, color="green"): """ create a plot polygon from the given coordinates""" @@ -60,12 +74,14 @@ def get_min_eq_dist(ppoly): """ min_eq_dist= list() for point in ppoly: - point['distance'] = haversine(point['lon'], 0, point['lon'], point['lat'], arc=False) - #print(point['lat'], point['distance']) - min_eq_dist.append(point['distance']) + point['x_atc'] = haversine(point['lon'], 0, point['lon'], point['lat'], arc=False) + #print(point['lat'], point['x_atc']) + min_eq_dist.append(point['x_atc']) return min(min_eq_dist)*1e3 # in meters. needed for defining the x-axis + + def create_polygons(latR, lonR): from shapely.geometry import Polygon, Point latR.sort() @@ -120,16 +136,16 @@ def ascending_test(track_data): def ascending_test_distance(track_data): """ - test if the track is ascending or descending based on 'distance' column + test if the track is ascending or descending based on 'x_atc' column """ # get the first and last point - first_point = abs(track_data.iloc[track_data['distance'].argmin()].geometry.y) - last_point = abs(track_data.iloc[track_data['distance'].argmax()].geometry.y) + first_point = abs(track_data.iloc[track_data['x_atc'].argmin()].geometry.y) + last_point = abs(track_data.iloc[track_data['x_atc'].argmax()].geometry.y) - # first_point = track_data.iloc[track_data.index.argmin()]['distance'] - # last_point = track_data.iloc[track_data.index.argmax()]['distance'] + # first_point = track_data.iloc[track_data.index.argmin()]['x_atc'] + # last_point = track_data.iloc[track_data.index.argmax()]['x_atc'] # get the time # first_time = cGPS.convert_GPS_time(first_point.index, first_point['rgt'], first_point['orbit_number']) # last_time = cGPS.convert_GPS_time(last_point['delta_time'], last_point['rgt'], last_point['orbit_number']) @@ -205,7 +221,7 @@ def plot_reference_point_coordinates(tmp, start_point_dist, start_point): ax.set_title('RGT: ' + str(rgt), loc='left') ax = axx[1] - ax.plot(tmp['distance'], tmp.geometry.y ,'.', markersize=0.5, c=data_col, label ='data') + ax.plot(tmp['x_atc'], tmp.geometry.y ,'.', markersize=0.5, c=data_col, label ='data') ax.plot(start_point_dist, start_point.geometry.y ,'.', c =spoint_color, label='reference point') ax.set_xlabel('Distance from Equator') @@ -316,10 +332,10 @@ def define_x_coordinate_in_polygon(table_data, polygon, round=True): min_eq_dist = np.round(get_min_eq_dist(polygon)) if ascending_test(table_data): - table_data['x'] = table_data['distance'] - min_eq_dist + table_data['x'] = table_data['x_atc'] - min_eq_dist else: #print('descending') - table_data['x'] = ((np.pi * 6371*1e3) - min_eq_dist) - table_data['distance'] + table_data['x'] = ((np.pi * 6371*1e3) - min_eq_dist) - table_data['x_atc'] return table_data @@ -343,9 +359,9 @@ def define_x_coordinate_with_RGT(table_data, Gtrack_lowest): start_point_dist, start_point = define_reference_distance_with_RGT(Gtrack_lowest, rgt, acending) if acending: - table_data['x'] = table_data['distance'] - start_point_dist + table_data['x'] = table_data['x_atc'] - start_point_dist else: - table_data['x'] = start_point_dist - table_data['distance'] + table_data['x'] = start_point_dist - table_data['x_atc'] return table_data @@ -360,12 +376,13 @@ def define_x_coordinate_from_data(table_data): acending = ascending_test_distance(table_data) if acending: - start_point_dist = table_data['distance'].min() - table_data['x'] = table_data['distance'] - start_point_dist + start_point_dist = table_data['x_atc'].min() + table_data['x'] = table_data['x_atc'] - start_point_dist else: - start_point_dist = table_data['distance'].max() - table_data['x'] = start_point_dist - table_data['distance'] - + start_point_dist = table_data['x_atc'].max() + table_data['x'] = start_point_dist - table_data['x_atc'] + table_data.sort_values(by='x', inplace=True) + table_data.reset_index(inplace=True) return table_data @@ -379,4 +396,4 @@ def create_ID_name(tt, segment=None): """ hemis= 'NH' if tt.geometry.y > 0 else 'SH' segment = '00' if segment is None else str(segment) - return hemis + '_' + str(tt.name.year) + str(tt.name.month).zfill(2) + str(tt.name.day).zfill(2) + '_' + str(tt.rgt).zfill(4) + str(tt.cycle).zfill(2) + segment + return hemis + '_' + str(tt.name.year) + str(tt.name.month).zfill(2) + str(tt.name.day).zfill(2) + '_' + str(tt.rgt).zfill(4) + str(tt.cycle).zfill(2) + segment \ No newline at end of file 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 new file mode 100644 index 00000000..a449f493 --- /dev/null +++ b/src/icesat2_tracks/analysis_db/B01_SL_load_single_file.py @@ -0,0 +1,270 @@ + +# %% +""" +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.11 +""" +import os, sys +from icesat2_tracks.config.IceSAT2_startup import mconfig,xr + +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 + +#import spicke_remover + +import h5py, imp, copy +import datetime + + +xr.set_options(display_style='text') + +# %matplotlib inline +# %matplotlib widget + + +# Select region and retrive batch of tracks + +track_name, batch_key, ID_flag = io.init_from_input(sys.argv) # loads standard experiment +# define file with ID: +#track_name, batch_key , ID_flag = '20190219073735_08070210_005_01', 'SH_testSLsinglefile2' , False +#track_name, batch_key , ID_flag = '20190219075059_08070212_005_01', 'SH_testSLsinglefile2' , False +#track_name, batch_key , ID_flag = '20190502052058_05180312_005_01', 'SH_testSLsinglefile2' , False + +#track_name, batch_key , ID_flag = '20190504201233_05580312_005_01', 'SH_testSLsinglefile2' , False + + +#20190502052058_05180312_005_01 +plot_flag = True +hemis = batch_key.split('_')[0] +#plot_path = mconfig['paths']['plot'] + '/'+hemis+'/'+batch_key+'/' + track_name +'/' + +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) + +#ID, _, _, _ = io.init_data(track_name, batch_key, True, mconfig['paths']['work'], ) +ATL03_track_name = 'ATL03_'+track_name+'.h5' +#track_name = ID['tracks']['ATL03'][0] +'.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 + + +# %% plot the ground tracks in geographic location +# Generate ATL06-type segments using the ATL03-native photon classification +# Use the ocean classification for photons with a confidence parmeter to 2 or higher (low confidence or better) + +params={'srt': 1, # Ocean classification + 'len': 25, # 10-meter segments + 'ats':3, # require that each segment contain photons separated by at least 5 m + 'res':10, # return one photon every 10 m + 'dist_in_seg': False, # if False units of len and res are in meters + 'track': 0, # return all ground tracks + 'pass_invalid': False, + 'cnt': 20, + 'sigma_r_max': 4, # maximum standard deviation of photon in extend + 'cnf': 2, # require classification confidence of 2 or more + 'atl03_geo_fields' : ['dem_h'] + } + + +# YAPC alternatibe +params_yapc={'srt': 1, + 'len': 20, + 'ats':3, + 'res':10, + 'dist_in_seg': False, # if False units of len and res are in meters + 'track': 0, + 'pass_invalid': False, + 'cnf': 2, + 'cnt': 20, + 'sigma_r_max': 4, # maximum standard deviation of photon in extend + 'maxi':10, + 'yapc': dict(knn=0, win_h=6, win_x=11, min_ph=4, score=100), # use the YAPC photon classifier; these are the recommended parameters, but the results might be more specific with a smaller win_h value, or a higher score cutoff +# "yapc": dict(knn=0, win_h=3, win_x=11, min_ph=4, score=50), # use the YAPC photon classifier; these are the recommended parameters, but the results might be more specific with a smaller win_h value, or a higher score cutoff +'atl03_geo_fields' : ['dem_h'] +} + +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)#.plot(markersize=0.2) + + +# %% +cdict = dict() +for s,b in zip(gdf['spot'].unique(), ['gt1l', 'gt1r', 'gt2l', 'gt2r', 'gt3l', 'gt3r']): + cdict[s] = col.rels[b] + +#imp.reload(beam_stats) + +font_for_pres() +F_atl06 = M.figure_axis_xy(6.5, 5, view_scale=0.6) +F_atl06.fig.suptitle(track_name) + +beam_stats.plot_ATL06_track_data(gdf, cdict) + + +# %% 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 + inputs: + table_data: GeoDataFrame with the data + split_by_beam: True/False. If True the data is split by beam + returns: + if split_by_beam: + table_data: dict of GeoDataFrame with the data for each beam + else: + table_data: GeoDataFrame with the data for all beams in one table + """ + #table_data = copy.copy(B01b_atl06_native) + + table_data.rename(columns= { + 'n_fit_photons':'N_photos', + 'w_surface_window_final':'signal_confidence', + 'y_atc':'y', + 'x_atc': 'distance', + #'spot':'beam', + } + , inplace=True) + + table_data['lons'] = table_data['geometry'].x + table_data['lats'] = table_data['geometry'].y + + drop_columns = ['cycle','gt','rgt', 'pflags'] + if to_hdf5: + drop_columns.append('geometry') + table_data.drop(columns=drop_columns, inplace=True) + + if split_by_beam: + B01b = dict() + # this is not tested + for spot,beam in zip([1,2,3, 4, 5, 6],['gt1l', 'gt1r', 'gt2l', 'gt2r', 'gt3l', 'gt3r']) : + ii = table_data.spot==spot + B01b[beam] = table_data[ii] + return B01b + else: + return table_data + +# %% +# 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) +table_time = table_data['time'] +table_data.drop(columns=['time'], inplace=True) +#table_data['time'] = np.datetime_as_string(table_data['time'])#.view(' Date: Wed, 8 Nov 2023 15:38:57 -0500 Subject: [PATCH 2/7] all necessary changes to run step 1 --- src/icesat2_tracks/analysis_db/B01_SL_load_single_file.py | 7 ++++--- src/icesat2_tracks/config/IceSAT2_startup.py | 2 +- 2 files changed, 5 insertions(+), 4 deletions(-) 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 a449f493..f6686900 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 @@ -5,7 +5,7 @@ This is python 3.11 """ import os, sys -from icesat2_tracks.config.IceSAT2_startup import mconfig,xr +from icesat2_tracks.config.IceSAT2_startup import mconfig,xr,color_schemes,font_for_pres,plt import geopandas as gpd @@ -18,6 +18,7 @@ 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 spicke_remover @@ -107,7 +108,7 @@ # %% cdict = dict() for s,b in zip(gdf['spot'].unique(), ['gt1l', 'gt1r', 'gt2l', 'gt2r', 'gt3l', 'gt3r']): - cdict[s] = col.rels[b] + cdict[s] = color_schemes.rels[b] #imp.reload(beam_stats) @@ -215,7 +216,7 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): font_for_pres() F = M.figure_axis_xy(8, 4.3, view_scale= 0.6 ) - beam_stats.plot_beam_statistics(D, high_beams, low_beams, col.rels, track_name = track_name + '| ascending =' + str(sct.ascending_test_distance(gdf)) ) + beam_stats.plot_beam_statistics(D, high_beams, low_beams, color_schemes.rels, track_name = track_name + '| ascending =' + str(sct.ascending_test_distance(gdf)) ) F.save_light(path = plot_path , name = 'B01b_beam_statistics.png') plt.close() diff --git a/src/icesat2_tracks/config/IceSAT2_startup.py b/src/icesat2_tracks/config/IceSAT2_startup.py index 2614093c..30ffaffb 100644 --- a/src/icesat2_tracks/config/IceSAT2_startup.py +++ b/src/icesat2_tracks/config/IceSAT2_startup.py @@ -36,7 +36,7 @@ mconfig["paths"].update({"config": config_dir_path}) #load colorscheme -col=M_color.color(path=mconfig['paths']['config'], name='color_def') +color_schemes=M_color.color(path=mconfig['paths']['config'], name='color_def') lstrings =iter([i+') ' for i in list(string.ascii_lowercase)]) From cc91553cf3c3e6f645072b16269c3b71388c3c6b Mon Sep 17 00:00:00 2001 From: Camilo Diaz Date: Wed, 8 Nov 2023 15:44:45 -0500 Subject: [PATCH 3/7] removed B01_SL_load_single_file.py from root folder --- analysis_db/B01_SL_load_single_file.py | 144 ------------------------- 1 file changed, 144 deletions(-) delete mode 100644 analysis_db/B01_SL_load_single_file.py diff --git a/analysis_db/B01_SL_load_single_file.py b/analysis_db/B01_SL_load_single_file.py deleted file mode 100644 index cb800a31..00000000 --- a/analysis_db/B01_SL_load_single_file.py +++ /dev/null @@ -1,144 +0,0 @@ - -# %% -""" -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.11 -""" -import os, sys -exec(open(os.environ["PYTHONSTARTUP"]).read()) -exec(open(STARTUP_2021_IceSAT2).read()) - -import geopandas as gpd - -from sliderule import icesat2 -from sliderule import sliderule - -import shapely -from ipyleaflet import basemaps, Map, GeoData - -import ICEsat2_SI_tools.sliderule_converter_tools as sct -import ICEsat2_SI_tools.io as io - -import spicke_remover - -import h5py, imp, copy - -xr.set_options(display_style='text') - -# %matplotlib inline -# %matplotlib widget - - -# %% Select region and retrive batch of tracks - -# define file with ID: -track_name, batch_key , ID_flag = '20190219073735_08070210_005_01', 'SH_testSLsinglefile' , False - -plot_flag = True - - -save_path = mconfig['paths']['work'] +'/'+batch_key+'/B01_regrid/' -MT.mkdirs_r(save_path) - -#ID, _, _, _ = io.init_data(track_name, batch_key, True, mconfig['paths']['work'], ) -ATL03_track_name = 'ATL03_'+track_name+'.h5' -#track_name = ID['tracks']['ATL03'][0] +'.h5' - -# %% Configure SL Session # -icesat2.init("slideruleearth.io", True) -asset = 'nsidc-s3' - -# %% plot the ground tracks in geographic location -# Generate ATL06-type segments using the ATL03-native photon classification -# Use the ocean classification for photons with a confidence parmeter to 2 or higher (low confidence or better) - -params={'srt': 1, # Ocean classification - 'len': 25, # 10-meter segments - 'ats':3, # require that each segment contain photons separated by at least 5 m - 'res':10, # return one photon every 10 m - 'track': 0, # return all ground tracks - 'pass_invalid': True, - 'cnf': 2, # require classification confidence of 2 or more - #'iterations':10, # iterate the fit -} - -# YAPC alternatibe -params_yapc={'srt': 1, - 'len': 20, - 'ats':3, - 'res':10, - 'track': 0, - 'pass_invalid': True, - 'cnf': -2, - 'iterations':10,} -# # "yapc": dict(knn=0, win_h=6, win_x=11, min_ph=4, score=100), # use the YAPC photon classifier; these are the recommended parameters, but the results might be more specific with a smaller win_h value, or a higher score cutoff -# "yapc": dict(knn=0, win_h=3, win_x=11, min_ph=4, score=50), # use the YAPC photon classifier; these are the recommended parameters, but the results might be more specific with a smaller win_h value, or a higher score cutoff - -gdf = icesat2.atl06p(params, asset="nsidc-s3", resources=[ATL03_track_name]) - -gdf.plot() -# %% 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 - inputs: - table_data: GeoDataFrame with the data - split_by_beam: True/False. If True the data is split by beam - returns: - if split_by_beam: - table_data: dict of GeoDataFrame with the data for each beam - else: - table_data: GeoDataFrame with the data for all beams in one table - """ - #table_data = copy.copy(B01b_atl06_native) - - table_data.rename(columns= { - 'n_fit_photons':'N_photos', - 'w_surface_window_final':'signal_confidence', - 'across':'y', - #'spot':'beam', - } - , inplace=True) - - table_data['lons'] = table_data['geometry'].x - table_data['lats'] = table_data['geometry'].y - - drop_columns = ['cycle','gt','rgt', 'pflags'] - if to_hdf5: - drop_columns.append('geometry') - table_data.drop(columns=drop_columns, inplace=True) - - if split_by_beam: - B01b = dict() - # this is not tested - for spot,beam in zip([1,2,3, 4, 5, 6],['gt1l', 'gt1r', 'gt2l', 'gt2r', 'gt3l', 'gt3r']) : - ii = table_data.spot==spot - B01b[beam] = table_data[ii] - return B01b - else: - return table_data - -# define reference point and then define 'x' -table_data = copy.copy(gdf) - -# 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) - -# fake 'across' -table_data['across'] = table_data['x']*0 +table_data['spot'] - -# add spike remover - -# renames columns and splits beams -Ti = make_B01_dict(table_data, split_by_beam=True, to_hdf5=True) -#Ti['gt1l'].drop('geometry', axis=1, inplace=True) - -segment = track_name.split('_')[1][-2:] -ID_name = sct.create_ID_name(gdf.iloc[0], segment=segment) -print( ID_name ) -io.write_track_to_HDF5(Ti, ID_name + '_B01_binned' , save_path) # regridding heights - -print('done') -# %% From a59ca4f21d2c31d4094e33d3646def61c6db99dd Mon Sep 17 00:00:00 2001 From: Camilo Diaz Date: Fri, 10 Nov 2023 09:59:30 -0500 Subject: [PATCH 4/7] removing commented lines --- .../ICEsat2_SI_tools/beam_stats.py | 54 +------------------ .../sliderule_converter_tools.py | 2 - .../analysis_db/B01_SL_load_single_file.py | 38 ++----------- 3 files changed, 5 insertions(+), 89 deletions(-) diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py b/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py index 66950e2b..257aad1b 100644 --- a/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py +++ b/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py @@ -88,9 +88,8 @@ def plot_beam_statistics(D, high_beams, low_beams, col_dict, track_name =None): plt.plot(D[k]['x']/1e3, np.sqrt(D[k]['var']), '.', color= col_dict[k], markersize=4, label=k) plt.title('high beams std', loc='left') - #plt.xlabel('along track distance (km)') plt.ylabel('segment std log(m)') - #plt.ylim(0, 20) + ax1.set_yscale('log') ax2 = plt.subplot(gs[1, 0]) @@ -108,9 +107,7 @@ def plot_beam_statistics(D, high_beams, low_beams, col_dict, track_name =None): plt.plot(D[k]['x']/1e3, np.sqrt(D[k]['var']), '.', color= col_dict[k], markersize=4, label=k) plt.title('low beams std', loc='left') - #plt.xlabel('along track distance (km)') - #plt.ylabel('segment std (m)') - #plt.ylim(0, 20) + ax3.set_yscale('log') ax4 = plt.subplot(gs[1, 1]) @@ -142,53 +139,6 @@ def plot_beam_statistics(D, high_beams, low_beams, col_dict, track_name =None): plt.legend() plt.show() - - # # make 2 x 2 plot - # plt.subplot(2, 3, 1) - # for k in high_beams: - # plt.plot(D[k]['x']/1e3, np.sqrt(D[k]['var']), '.', color= col_dict[k], markersize=4, label=k) - - # plt.title('high beams std', loc='left') - # #plt.xlabel('along track distance (km)') - # plt.ylabel('segment std (m)') - # plt.ylim(0, 20) - - # plt.subplot(2, 3, 3) - # for k in high_beams: - # Di = D[k]['N'] - # Di[Di ==0] =np.nan - # plt.plot(D[k]['lat'], D[k]['N'], '.', color= col_dict[k], markersize=4, label=k) - - # plt.title('high beams N', loc='left') - # plt.xlabel('along track distance (km) or Lat') - # plt.ylabel('Point Density (m)') - - - # plt.subplot(2, 3, 2) - # for k in low_beams: - # plt.plot(D[k]['x']/1e3, np.sqrt(D[k]['var']), '.', color= col_dict[k], markersize=4, label=k) - - # plt.title('low beams std', loc='left') - # #plt.xlabel('along track distance (km)') - # #plt.ylabel('segment std (m)') - # plt.ylim(0, 20) - - # plt.subplot(2, 3, 4) - # for k in low_beams: - # Di = D[k]['N'] - # Di[Di ==0] =np.nan - # plt.plot(D[k]['lat'], D[k]['N'], '.', color= col_dict[k], markersize=4, label=k) - - # plt.title('low beams N', loc='left') - # plt.xlabel('along track distance (km) or Lat') - # #plt.ylabel('Point density (m)') - # #return F - - # plt.subplot(2,3, 5) - # plt.scatter(D['gt1l']['x']/1e3, D['gt1l']['lat'], s= np.exp(D['gt1l']['N'] *10) , marker='.', color= col_dict['gt1l'], label='gt1l', alpha = 0.3) - - # plt.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, wspace=0.3, hspace=0.3) - ## plot track stats basics for sliderules ATL06 output def plot_ATL06_track_data( G2, cdict): diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py b/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py index c5af7029..fb1d34b1 100644 --- a/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py +++ b/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py @@ -75,7 +75,6 @@ def get_min_eq_dist(ppoly): min_eq_dist= list() for point in ppoly: point['x_atc'] = haversine(point['lon'], 0, point['lon'], point['lat'], arc=False) - #print(point['lat'], point['x_atc']) min_eq_dist.append(point['x_atc']) return min(min_eq_dist)*1e3 # in meters. needed for defining the x-axis @@ -334,7 +333,6 @@ def define_x_coordinate_in_polygon(table_data, polygon, round=True): if ascending_test(table_data): table_data['x'] = table_data['x_atc'] - min_eq_dist else: - #print('descending') table_data['x'] = ((np.pi * 6371*1e3) - min_eq_dist) - table_data['x_atc'] return table_data 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 f6686900..3c48c152 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,7 +1,7 @@ # %% """ -This file open a ICEsat2 track applied filters and corections and returns smoothed photon heights on a regular grid in an .nc file. +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 @@ -20,7 +20,6 @@ import icesat2_tracks.local_modules.m_tools_ph3 as MT from icesat2_tracks.local_modules import m_general_ph3 as M -#import spicke_remover import h5py, imp, copy import datetime @@ -28,25 +27,15 @@ xr.set_options(display_style='text') -# %matplotlib inline -# %matplotlib widget - # Select region and retrive batch of tracks track_name, batch_key, ID_flag = io.init_from_input(sys.argv) # loads standard experiment -# define file with ID: -#track_name, batch_key , ID_flag = '20190219073735_08070210_005_01', 'SH_testSLsinglefile2' , False -#track_name, batch_key , ID_flag = '20190219075059_08070212_005_01', 'SH_testSLsinglefile2' , False -#track_name, batch_key , ID_flag = '20190502052058_05180312_005_01', 'SH_testSLsinglefile2' , False - -#track_name, batch_key , ID_flag = '20190504201233_05580312_005_01', 'SH_testSLsinglefile2' , False - #20190502052058_05180312_005_01 plot_flag = True hemis = batch_key.split('_')[0] -#plot_path = mconfig['paths']['plot'] + '/'+hemis+'/'+batch_key+'/' + track_name +'/' + save_path = mconfig['paths']['work'] +'/'+batch_key+'/B01_regrid/' MT.mkdirs_r(save_path) @@ -54,9 +43,7 @@ save_path_json = mconfig['paths']['work'] +'/'+ batch_key +'/A01b_ID/' MT.mkdirs_r(save_path_json) -#ID, _, _, _ = io.init_data(track_name, batch_key, True, mconfig['paths']['work'], ) ATL03_track_name = 'ATL03_'+track_name+'.h5' -#track_name = ID['tracks']['ATL03'][0] +'.h5' # %% Configure SL Session # sliderule.authenticate("brown", ps_username="mhell", ps_password="Oijaeth9quuh") @@ -110,7 +97,6 @@ for s,b in zip(gdf['spot'].unique(), ['gt1l', 'gt1r', 'gt2l', 'gt2r', 'gt3l', 'gt3r']): cdict[s] = color_schemes.rels[b] -#imp.reload(beam_stats) font_for_pres() F_atl06 = M.figure_axis_xy(6.5, 5, view_scale=0.6) @@ -133,14 +119,12 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): else: table_data: GeoDataFrame with the data for all beams in one table """ - #table_data = copy.copy(B01b_atl06_native) - + table_data.rename(columns= { 'n_fit_photons':'N_photos', 'w_surface_window_final':'signal_confidence', 'y_atc':'y', 'x_atc': 'distance', - #'spot':'beam', } , inplace=True) @@ -171,13 +155,6 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): table_data = sct.define_x_coordinate_from_data(table_data) table_time = table_data['time'] table_data.drop(columns=['time'], inplace=True) -#table_data['time'] = np.datetime_as_string(table_data['time'])#.view(' Date: Fri, 10 Nov 2023 15:25:35 -0500 Subject: [PATCH 5/7] Apply suggestions from code review Co-authored-by: Carlos Paniagua --- src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py | 1 - src/icesat2_tracks/analysis_db/B01_SL_load_single_file.py | 4 +--- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py b/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py index 257aad1b..7a146e8b 100644 --- a/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py +++ b/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py @@ -1,7 +1,6 @@ import numpy as np import pandas as pd import icesat2_tracks.ICEsat2_SI_tools.spectral_estimates as spec -#import io as io_local import icesat2_tracks.ICEsat2_SI_tools.io as io_local import matplotlib.pyplot as plt 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 3c48c152..daebccdf 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 @@ -89,7 +89,7 @@ print("STARTS") gdf = icesat2.atl06p(params_yapc, resources=[ATL03_track_name]) print("ENDS") -gdf = sct.correct_and_remove_height(gdf, maximum_height)#.plot(markersize=0.2) +gdf = sct.correct_and_remove_height(gdf, maximum_height) # %% @@ -199,9 +199,7 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): # plot the ground tracks in geographic location gdf[::100].plot(markersize=0.1, figsize=(4,6)) plt.title(track_name + '\nascending =' + str(sct.ascending_test_distance(gdf)) , loc ='left') - # gdf_atl03 = icesat2.atl03s(params, [ATL03_track_name]) M.save_anyfig(plt.gcf(), path = plot_path ,name = 'B01_track.png') - # gdf_atl03.plot() plt.close() From 71ad3c7e29db77e0a78e21d38723a754b30e1276 Mon Sep 17 00:00:00 2001 From: Camilo Diaz Date: Tue, 14 Nov 2023 12:33:28 -0500 Subject: [PATCH 6/7] removed matlab marks --- .../analysis_db/B01_SL_load_single_file.py | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) 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 3c48c152..dedb18de 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,5 +1,5 @@ -# %% +# """ 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 @@ -45,12 +45,12 @@ ATL03_track_name = 'ATL03_'+track_name+'.h5' -# %% Configure SL Session # +# 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 -# %% plot the ground tracks in geographic location +# plot the ground tracks in geographic location # Generate ATL06-type segments using the ATL03-native photon classification # Use the ocean classification for photons with a confidence parmeter to 2 or higher (low confidence or better) @@ -92,7 +92,6 @@ gdf = sct.correct_and_remove_height(gdf, maximum_height)#.plot(markersize=0.2) -# %% cdict = dict() for s,b in zip(gdf['spot'].unique(), ['gt1l', 'gt1r', 'gt2l', 'gt2r', 'gt3l', 'gt3r']): cdict[s] = color_schemes.rels[b] @@ -105,7 +104,7 @@ beam_stats.plot_ATL06_track_data(gdf, cdict) -# %% main routine for defining the x coordinate and sacing table data +# main routine for defining the x coordinate and sacing table data def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): """ @@ -146,7 +145,6 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): else: return table_data -# %% # define reference point and then define 'x' table_data = copy.copy(gdf) imp.reload(sct) @@ -170,7 +168,7 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): print( ID_name ) io.write_track_to_HDF5(Ti, ID_name + '_B01_binned' , save_path) # regridding heights -# %% plot the ground tracks in geographic location +# plot the ground tracks in geographic location all_beams = mconfig['beams']['all_beams'] high_beams = mconfig['beams']['high_beams'] @@ -178,8 +176,6 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): D = beam_stats.derive_beam_statistics(Ti, all_beams, Lmeter=12.5e3, dx =10) -# %% - # save figure from above: plot_path = mconfig['paths']['plot'] + '/'+hemis+'/'+batch_key+'/' + ID_name +'/' MT.mkdirs_r(plot_path) @@ -206,8 +202,6 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): - -# %% print('write A01b .json') DD= {'case_ID': ID_name , 'tracks' : {} } @@ -235,5 +229,3 @@ def make_B01_dict(table_data, split_by_beam=True, to_hdf5=False): MT.json_save2(name='A01b_ID_'+ID_name, path=save_path_json, data= DD) print('done') - -# %% \ No newline at end of file From 33520579fee53e7f5d5191916ac6829d97f1eb9d Mon Sep 17 00:00:00 2001 From: Camilo Diaz Date: Tue, 14 Nov 2023 12:49:50 -0500 Subject: [PATCH 7/7] removed commented code --- .../ICEsat2_SI_tools/sliderule_converter_tools.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py b/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py index fb1d34b1..31b274bf 100644 --- a/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py +++ b/src/icesat2_tracks/ICEsat2_SI_tools/sliderule_converter_tools.py @@ -138,17 +138,10 @@ def ascending_test_distance(track_data): test if the track is ascending or descending based on 'x_atc' column """ # get the first and last point - first_point = abs(track_data.iloc[track_data['x_atc'].argmin()].geometry.y) last_point = abs(track_data.iloc[track_data['x_atc'].argmax()].geometry.y) - # first_point = track_data.iloc[track_data.index.argmin()]['x_atc'] - # last_point = track_data.iloc[track_data.index.argmax()]['x_atc'] - # get the time - # first_time = cGPS.convert_GPS_time(first_point.index, first_point['rgt'], first_point['orbit_number']) - # last_time = cGPS.convert_GPS_time(last_point['delta_time'], last_point['rgt'], last_point['orbit_number']) - # test if ascending or descending if first_point < last_point: return True else: