# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python [conda env:notebook] *
#     language: python
#     name: conda-env-notebook-py
# ---

# # demo downloading data
#

# following icepyx tutorial:

from icepyx import icesat2data as ipd
import os
import shutil
from pprint import pprint
# %matplotlib inline

# Whillans Ice Plain
short_name = 'ATL06'
spatial_extent = [-162, -85, -147, -83]
date_range = ['2019-06-18','2019-06-25'] # leave some short range

# +
region_a = ipd.Icesat2Data(short_name, spatial_extent, date_range, \
                           start_time='03:30:00', end_time='21:30:00')#, version='001')

print(region_a.dataset)
print(region_a.dates)
print(region_a.start_time)
print(region_a.end_time)
print(region_a.dataset_version)
print(region_a.spatial_extent)
# -

region_a = ipd.Icesat2Data(short_name, spatial_extent, date_range)

print(region_a.dataset)
print(region_a.dates)
print(region_a.start_time)
print(region_a.end_time)
print(region_a.dataset_version)
print(region_a.spatial_extent)

region_a.visualize_spatial_extent()

print(region_a.latest_version())
region_a.dataset_summary_info()

# + jupyter={"outputs_hidden": true}
region_a.dataset_all_info()
# -

#build and view the parameters that will be submitted in our query
region_a.CMRparams

#get a list of the available granule IDs that meet your search criteria
region_a.avail_granules(ids=True)

# + jupyter={"outputs_hidden": true}
#print detailed information about the returned search results
region_a.granules.avail

# +
# Fill this in with your Earthdata Login user name and associated email

earthdata_uid = 'cgbarcheck'
email = '<EMAIL>'
region_a.earthdata_login(earthdata_uid, email)
# -

print(region_a.reqparams)

region_a.subsetparams()

# +
# region_a.show_custom_options(dictview=True)

# + jupyter={"outputs_hidden": true}
region_a.order_vars.avail(options=True)
# -

# Default variables to order:
#

# + jupyter={"outputs_hidden": true}
region_a.order_vars.remove(all=True)
region_a.order_vars.append(defaults=True)
pprint(region_a.order_vars.wanted)
# print(region_a.order_vars.wanted['latitude'])
# print(region_a.order_vars.wanted['longitude'])
# print(region_a.order_vars.wanted['h_li'])

# + [markdown] jupyter={"source_hidden": true}
# lets stick with default order variables for now
# -

# download data

# +
# Without variable subsetting, or with variable subsetting if you have run region_a.subsetparams(Coverage=region_a.order_vars.wanted):
# region_a.order_granules()

# # With variable subsetting, if this is your first use of the Coverage keyword argument (kwarg):
region_a.order_granules(Coverage=region_a.order_vars.wanted)
# -

# # Download data

# +
# # !mkdir ~/downloads

path = '/home/jovyan/downloads'

# without variable subsetting, or with variable subsetting if you have run region_a.order_granules(Coverage=region_a.order_vars.wanted)
region_a.download_granules(path)

# with variable subsetting, if this is your first use of the Coverage keyword argument (kwarg)
# region_a.download_granules(Coverage=region_a.order_vars.wanted)
# -

# # look at the data, starting with the demo for hdf5

# pull out x,y,z,t and whatever else was needed for Johan's tutorial; try and follow it



# +
print('product:    ', region_a.dataset)
print('dates:      ', region_a.dates)
print('start time: ', region_a.start_time)
print('end time:   ', region_a.end_time)
print('version:    ', region_a.dataset_version)
print('extent:     ', region_a.spatial_extent)

print('\nDATA:')
print('\n'.join([str(item) for item in region_a.avail_granules().items()]))

region_a.visualize_spatial_extent()

# +
from pathlib import Path
import h5py


# look at a file
data_home = Path(path)
files = list(data_home.glob('*.h5'))

for f in files[:10]: print(f)
# print('Total number of files:', len(files))

# !h5ls -r {files[4]} 
# -



# utility funcitons:
#

# +
import pyproj
from astropy.time import Time

def gps2dyr(time):
    """Converte GPS time to decimal years."""
    return Time(time, format='gps').decimalyear


def orbit_type(time, lat, tmax=1):
    """Separate tracks into ascending and descending.
    
    Defines tracks as segments with time breaks > tmax,
    and tests whether lat increases or decreases w/time.
    """
    tracks = np.zeros(lat.shape)  # generate track segment
    tracks[0:np.argmax(np.abs(lat))] = 1  # set values for segment
    is_asc = np.zeros(tracks.shape, dtype=bool)  # output index array

    # Loop trough individual secments
    for track in np.unique(tracks):
    
        i_track, = np.where(track == tracks)  # get all pts from seg
    
        if len(i_track) < 2: continue
    
        # Test if lat increases (asc) or decreases (des) w/time
        i_min = time[i_track].argmin()
        i_max = time[i_track].argmax()
        lat_diff = lat[i_track][i_max] - lat[i_track][i_min]
    
        # Determine track type
        if lat_diff > 0:  is_asc[i_track] = True
    
    return is_asc


def transform_coord(proj1, proj2, x, y):
    """Transform coordinates from proj1 to proj2 (EPSG num).

    Example EPSG projections:
        Geodetic (lon/lat): 4326
        Polar Stereo AnIS (x/y): 3031
        Polar Stereo GrIS (x/y): 3413
    """
    # Set full EPSG projection strings
    proj1 = pyproj.Proj("+init=EPSG:"+str(proj1))
    proj2 = pyproj.Proj("+init=EPSG:"+str(proj2))
    return pyproj.transform(proj1, proj2, x, y)  # convert


# -

# The simple reader

# +
import h5py
import numpy as np

def read_atl06(fname, outdir='data', bbox=None):
    """Read one ATL06 file and output 6 reduced files. 
    
    Extract variables of interest and separate the ATL06 file 
    into each beam (ground track) and ascending/descending orbits.
    """

    # Each beam is a group
    group = ['/gt1l', '/gt1r', '/gt2l', '/gt2r', '/gt3l', '/gt3r']

    # Loop trough beams
    for k, g in enumerate(group):
    
        #-----------------------------------#
        # 1) Read in data for a single beam #
        #-----------------------------------#
        
        data = {}
    
        try:
            # Load vars into memory (include as many as you want)
            with h5py.File(fname, 'r') as fi:
                
                data['lat'] = fi[g+'/land_ice_segments/latitude'][:]
                data['lon'] = fi[g+'/land_ice_segments/longitude'][:]
                data['h_li'] = fi[g+'/land_ice_segments/h_li'][:]
                data['s_li'] = fi[g+'/land_ice_segments/h_li_sigma'][:]
                data['t_dt'] = fi[g+'/land_ice_segments/delta_time'][:] # time elapsed from a reference time (below, t_ref)
                data['q_flag'] = fi[g+'/land_ice_segments/atl06_quality_summary'][:]
                data['s_fg'] = fi[g+'/land_ice_segments/fit_statistics/signal_selection_source'][:]
                data['snr'] = fi[g+'/land_ice_segments/fit_statistics/snr_significance'][:]
                data['h_rb'] = fi[g+'/land_ice_segments/fit_statistics/h_robust_sprd'][:]
                data['dac'] = fi[g+'/land_ice_segments/geophysical/dac'][:]
                data['f_sn'] = fi[g+'/land_ice_segments/geophysical/bsnow_conf'][:]
                data['dh_fit_dx'] = fi[g+'/land_ice_segments/fit_statistics/dh_fit_dx'][:]
                data['tide_earth'] = fi[g+'/land_ice_segments/geophysical/tide_earth'][:]
                data['tide_load'] = fi[g+'/land_ice_segments/geophysical/tide_load'][:]
                data['tide_ocean'] = fi[g+'/land_ice_segments/geophysical/tide_ocean'][:]
                data['tide_pole'] = fi[g+'/land_ice_segments/geophysical/tide_pole'][:]
                
                rgt = fi['/orbit_info/rgt'][:]                           # single value
                t_ref = fi['/ancillary_data/atlas_sdp_gps_epoch'][:]     # single value # reference time
                beam_type = fi[g].attrs["atlas_beam_type"].decode()      # strong/weak (str)
                spot_number = fi[g].attrs["atlas_spot_number"].decode()  # number (str)
                
        except:
            print('skeeping group:', g)
            print('in file:', fname)
            continue
            
        #---------------------------------------------#
        # 2) Filter data according region and quality #
        #---------------------------------------------#
        
        # Select a region of interest
        if bbox:
            lonmin, lonmax, latmin, latmax = bbox
            bbox_mask = (data['lon'] >= lonmin) & (data['lon'] <= lonmax) & \
                        (data['lat'] >= latmin) & (data['lat'] <= latmax)
        else:
            bbox_mask = np.ones_like(data['lat'], dtype=bool)  # get all
            
        # Only keep good data (quality flag + threshold + bbox) # keep only data less than 10000m high (if you are into surface features and not clouds)
        mask = (data['q_flag'] == 0) & (np.abs(data['h_li']) < 10e3) & (bbox_mask == 1)
        
        # If no data left, skeep
        if not any(mask): continue
        
        # Update data variables
        for k, v in data.items(): data[k] = v[mask]
            
        #----------------------------------------------------#
        # 3) Convert time, separate tracks, reproject coords #
        #----------------------------------------------------#
        
        # Time in GPS seconds (secs sinde Jan 5, 1980)
        t_gps = t_ref + data['t_dt']

        # Time in decimal years
        t_year = gps2dyr(t_gps)

        # Determine orbit type
        is_asc = orbit_type(t_year, data['lat'])
        
        # Geodetic lon/lat -> Polar Stereo x/y
        x, y = transform_coord(4326, 3031, data['lon'], data['lat'])
        
        data['x'] = x
        data['y'] = y
        data['t_gps'] = t_gps
        data['t_year'] = t_year
        data['is_asc'] = is_asc
        
        #-----------------------#
        # 4) Save selected data #
        #-----------------------#
        
        # Define output dir and file
        outdir = Path(outdir)    
        fname = Path(fname)
        outdir.mkdir(exist_ok=True)
        outfile = outdir / fname.name.replace('.h5', '_' + g[1:] + '.h5')
        
        # Save variables
        with h5py.File(outfile, 'w') as fo:
            for k, v in data.items(): fo[k] = v
            print('out ->', outfile)
                
# -









StarcoderdataJupyterScriptsDedupFiltered