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') -# %% 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..7a146e8b --- /dev/null +++ b/src/icesat2_tracks/ICEsat2_SI_tools/beam_stats.py @@ -0,0 +1,198 @@ +import numpy as np +import pandas as pd +import icesat2_tracks.ICEsat2_SI_tools.spectral_estimates as spec +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.ylabel('segment std log(m)') + + 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') + + 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() + +## 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..31b274bf 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,13 @@ 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) + 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,20 +135,13 @@ 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'] - # 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: @@ -205,7 +213,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 +324,9 @@ 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 +350,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 +367,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 +387,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..b279e89c --- /dev/null +++ b/src/icesat2_tracks/analysis_db/B01_SL_load_single_file.py @@ -0,0 +1,229 @@ + +# +""" +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 +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 + + +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) + +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 + + +# 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) + + +cdict = dict() +for s,b in zip(gdf['spot'].unique(), ['gt1l', 'gt1r', 'gt2l', 'gt2r', 'gt3l', 'gt3r']): + cdict[s] = color_schemes.rels[b] + + +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.rename(columns= { + 'n_fit_photons':'N_photos', + 'w_surface_window_final':'signal_confidence', + 'y_atc':'y', + 'x_atc': 'distance', + } + , 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) + +# renames columns and splits beams +Ti = make_B01_dict(table_data, split_by_beam=True, to_hdf5=True) + +for kk in Ti.keys(): + Ti[kk]['dist'] = Ti[kk]['x'].copy() + 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 ) +io.write_track_to_HDF5(Ti, ID_name + '_B01_binned' , save_path) # regridding heights + +# plot the ground tracks in geographic location + +all_beams = mconfig['beams']['all_beams'] +high_beams = mconfig['beams']['high_beams'] +low_beams = mconfig['beams']['low_beams'] + +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) +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 ) + 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() + + # 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') + M.save_anyfig(plt.gcf(), path = plot_path ,name = 'B01_track.png') + plt.close() + + + +print('write A01b .json') +DD= {'case_ID': ID_name , 'tracks' : {} } + +DD['tracks']['ATL03'] = 'ATL10-' +track_name + + +start_pos = abs(table_data.lats).argmin() +end_pos = abs(table_data.lats).argmax() + + + +# add other pars: +DD['pars'] ={ +'poleward': sct.ascending_test(gdf), 'region': '0', +'start': {'longitude': table_data.lons[start_pos], 'latitude': table_data.lats[start_pos] +, 'seg_dist_x': float(table_data.x[start_pos]) +, 'delta_time': datetime.datetime.timestamp(table_time[start_pos]) +}, +'end': {'longitude': table_data.lons[end_pos], 'latitude': table_data.lats[end_pos] +, 'seg_dist_x': float(table_data.x[end_pos]) +, 'delta_time': datetime.datetime.timestamp(table_time[end_pos]) +}, + } + +MT.json_save2(name='A01b_ID_'+ID_name, path=save_path_json, data= DD) + +print('done') 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)])