Skip to content

Commit 6a58cbb

Browse files
jpswinskijpswinski
committed
ancillary data support for both extents and photons (atl03 and atl06)
1 parent d5f3d3e commit 6a58cbb

File tree

5 files changed

+197
-114
lines changed

5 files changed

+197
-114
lines changed

‎sliderule/icesat2.py

Lines changed: 119 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -745,25 +745,25 @@ def atl06p(parm, asset=DEFAULT_ASSET, version=DEFAULT_ICESAT2_SDP_VERSION, callb
745745
columns = {}
746746
elevation_records = []
747747
num_elevations = 0
748-
field_dictionary = {}
748+
field_dictionary = {} # ['field_name'] = {"extent_id": [], field_name: []}
749749
if len(rsps) > 0:
750750
# Sort Records
751751
for rsp in rsps:
752752
if 'atl06rec' in rsp['__rectype']:
753753
elevation_records += rsp,
754754
num_elevations += len(rsp['elevation'])
755-
elif 'ga3rec' == rsp['__rectype']:
755+
elif 'extrec' == rsp['__rectype']:
756756
field_name = parm['atl03_geo_fields'][rsp['field_index']]
757757
if field_name not in field_dictionary:
758758
field_dictionary[field_name] = {"extent_id": [], field_name: []}
759759
# Parse Ancillary Data
760760
data = __get_values(rsp['data'], rsp['data_type'], len(rsp['data']))
761761
# Add Left Pair Track Entry
762762
field_dictionary[field_name]['extent_id'] += rsp['extent_id'] | 0x2,
763-
field_dictionary[field_name][field_name] += data[0],
763+
field_dictionary[field_name][field_name] += data[LEFT_PAIR],
764764
# Add Right Pair Track Entry
765765
field_dictionary[field_name]['extent_id'] += rsp['extent_id'] | 0x3,
766-
field_dictionary[field_name][field_name] += data[1],
766+
field_dictionary[field_name][field_name] += data[RIGHT_PAIR],
767767
# Build Elevation Columns
768768
if num_elevations > 0:
769769
# Initialize Columns
@@ -781,6 +781,7 @@ def atl06p(parm, asset=DEFAULT_ASSET, version=DEFAULT_ICESAT2_SDP_VERSION, callb
781781
elev_cnt += 1
782782
else:
783783
logger.debug("No response returned")
784+
784785
profiles["flatten"] = time.perf_counter() - tstart_flatten
785786

786787
# Build GeoDataFrame
@@ -878,71 +879,125 @@ def atl03sp(parm, asset=DEFAULT_ASSET, version=DEFAULT_ICESAT2_SDP_VERSION, call
878879
# Flatten Responses
879880
tstart_flatten = time.perf_counter()
880881
columns = {}
881-
if len(rsps) <= 0:
882-
logger.debug("no response returned")
883-
elif rsps[0]['__rectype'] != 'atl03rec':
884-
logger.debug("invalid response returned: %s", rsps[0]['__rectype'])
885-
else:
886-
# Count Rows
887-
num_rows = 0
888-
for rsp in rsps:
889-
num_rows += len(rsp["data"])
890-
# Build Columns
891-
for rsp in rsps:
892-
if len(rsp["data"]) > 0:
893-
# Allocate Columns
894-
for field in rsp.keys():
895-
fielddef = sliderule.get_definition("atl03rec", field)
896-
if len(fielddef) > 0:
897-
columns[field] = numpy.empty(num_rows, fielddef["nptype"])
898-
for field in rsp["data"][0].keys():
899-
fielddef = sliderule.get_definition("atl03rec.photons", field)
900-
if len(fielddef) > 0:
901-
columns[field] = numpy.empty(num_rows, fielddef["nptype"])
902-
break
903-
# Populate Columns
904-
ph_cnt = 0
882+
sample_photon_record = None
883+
photon_records = []
884+
num_photons = 0
885+
extent_dictionary = {}
886+
extent_field_types = {} # ['field_name'] = nptype
887+
photon_dictionary = {}
888+
photon_field_types = {} # ['field_name'] = nptype
889+
if len(rsps) > 0:
890+
# Sort Records
905891
for rsp in rsps:
906-
ph_index = 0
907-
pair = 0
908-
left_cnt = rsp["count"][0]
909-
for photon in rsp["data"]:
910-
if ph_index >= left_cnt:
911-
pair = 1
912-
for field in rsp.keys():
913-
if field in columns:
914-
if field == "count":
915-
columns[field][ph_cnt] = pair
916-
elif type(rsp[field]) is tuple:
917-
columns[field][ph_cnt] = rsp[field][pair]
918-
else:
919-
columns[field][ph_cnt] = rsp[field]
920-
for field in photon.keys():
921-
if field in columns:
922-
columns[field][ph_cnt] = photon[field]
923-
ph_cnt += 1
924-
ph_index += 1
925-
# Rename Count Column to Pair Column
926-
columns["pair"] = columns.pop("count")
927-
profiles["flatten"] = time.perf_counter() - tstart_flatten
928-
929-
# Create DataFrame
930-
df = __todataframe(columns, "delta_time", "longitude", "latitude")
931-
932-
# Calculate Spot Column
933-
df['spot'] = df.apply(lambda row: __calcspot(row["sc_orient"], row["track"], row["pair"]), axis=1)
934-
935-
# Return Response
936-
profiles[atl03sp.__name__] = time.perf_counter() - tstart
937-
return df
938-
939-
# Error Case
940-
return __emptyframe()
892+
extent_id = rsp['extent_id']
893+
if 'atl03rec' in rsp['__rectype']:
894+
photon_records += rsp,
895+
num_photons += len(rsp['data'])
896+
if sample_photon_record == None and len(rsp['data']) > 0:
897+
sample_photon_record = rsp
898+
elif 'extrec' == rsp['__rectype']:
899+
# Get Field Type
900+
field_name = parm['atl03_geo_fields'][rsp['field_index']]
901+
if field_name not in extent_field_types:
902+
extent_field_types[field_name] = sliderule.basictypes[sliderule.codedtype2str[rsp['data_type']]]["nptype"]
903+
# Initialize Extent Dictionary Entry
904+
if extent_id not in extent_dictionary:
905+
extent_dictionary[extent_id] = {}
906+
# Save of Values per Extent ID per Field Name
907+
data = __get_values(rsp['data'], rsp['data_type'], len(rsp['data']))
908+
extent_dictionary[extent_id][field_name] = data
909+
elif 'phrec' == rsp['__rectype']:
910+
# Get Field Type
911+
field_name = parm['atl03_ph_fields'][rsp['field_index']]
912+
if field_name not in photon_field_types:
913+
photon_field_types[field_name] = sliderule.basictypes[sliderule.codedtype2str[rsp['data_type']]]["nptype"]
914+
# Initialize Extent Dictionary Entry
915+
if extent_id not in photon_dictionary:
916+
photon_dictionary[extent_id] = {}
917+
# Save of Values per Extent ID per Field Name
918+
data = __get_values(rsp['data'], rsp['data_type'], len(rsp['data']))
919+
photon_dictionary[extent_id][field_name] = data
920+
# Build Elevation Columns
921+
if num_photons > 0:
922+
# Initialize Columns
923+
for field in sample_photon_record.keys():
924+
fielddef = sliderule.get_definition("atl03rec", field)
925+
if len(fielddef) > 0:
926+
columns[field] = numpy.empty(num_photons, fielddef["nptype"])
927+
for field in sample_photon_record["data"][0].keys():
928+
fielddef = sliderule.get_definition("atl03rec.photons", field)
929+
if len(fielddef) > 0:
930+
columns[field] = numpy.empty(num_photons, fielddef["nptype"])
931+
for field in extent_field_types.keys():
932+
columns[field] = numpy.empty(num_photons, extent_field_types[field])
933+
for field in photon_field_types.keys():
934+
columns[field] = numpy.empty(num_photons, photon_field_types[field])
935+
# Populate Columns
936+
ph_cnt = 0
937+
for record in photon_records:
938+
ph_index = 0
939+
pair = 0
940+
left_cnt = record["count"][0]
941+
extent_id = record['extent_id']
942+
# Get Extent Fields to Add to Extent
943+
extent_field_dictionary = {}
944+
if extent_id in extent_dictionary:
945+
extent_field_dictionary = extent_dictionary[extent_id]
946+
# Get Photon Fields to Add to Extent
947+
photon_field_dictionary = {}
948+
if extent_id in photon_dictionary:
949+
photon_field_dictionary = photon_dictionary[extent_id]
950+
# For Each Photon in Extent
951+
for photon in record["data"]:
952+
if ph_index >= left_cnt:
953+
pair = 1
954+
# Add per Extent Fields
955+
for field in record.keys():
956+
if field in columns:
957+
if field == "count":
958+
columns[field][ph_cnt] = pair # count gets changed to pair id
959+
elif type(record[field]) is tuple:
960+
columns[field][ph_cnt] = record[field][pair]
961+
else:
962+
columns[field][ph_cnt] = record[field]
963+
# Add per Photon Fields
964+
for field in photon.keys():
965+
if field in columns:
966+
columns[field][ph_cnt] = photon[field]
967+
# Add Ancillary Extent Fields
968+
for field in extent_field_dictionary:
969+
columns[field][ph_cnt] = extent_field_dictionary[field][pair]
970+
# Add Ancillary Extent Fields
971+
for field in photon_field_dictionary:
972+
columns[field][ph_cnt] = photon_field_dictionary[field][ph_index]
973+
# Goto Next Photon
974+
ph_cnt += 1
975+
ph_index += 1
976+
# Rename Count Column to Pair Column
977+
columns["pair"] = columns.pop("count")
978+
979+
profiles["flatten"] = time.perf_counter() - tstart_flatten
980+
981+
# Create DataFrame
982+
gdf = __todataframe(columns, "delta_time", "longitude", "latitude")
983+
984+
# Calculate Spot Column
985+
gdf['spot'] = gdf.apply(lambda row: __calcspot(row["sc_orient"], row["track"], row["pair"]), axis=1)
986+
987+
# Return Response
988+
profiles[atl03sp.__name__] = time.perf_counter() - tstart
989+
return gdf
990+
else:
991+
logger.debug("No photons returned")
992+
else:
993+
logger.debug("No response returned")
941994

942995
# Handle Runtime Errors
943996
except RuntimeError as e:
944997
logger.critical(e)
945-
return __emptyframe()
998+
999+
# Error or No Data
1000+
return __emptyframe()
9461001

9471002
#
9481003
# H5

‎utils/query_elevations.py

Lines changed: 2 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -4,9 +4,8 @@
44

55
import sys
66
import logging
7-
import time
87
from sliderule import icesat2
9-
from utils import initialize_client
8+
from utils import initialize_client, display_statistics
109

1110
###############################################################################
1211
# MAIN
@@ -21,20 +20,7 @@
2120
parms, cfg = initialize_client(sys.argv)
2221

2322
# Request ATL06 Data
24-
tstart = time.perf_counter()
2523
gdf = icesat2.atl06p(parms, asset=cfg["asset"], resources=[cfg["resource"]])
26-
perf_duration = time.perf_counter() - tstart
2724

2825
# Display Statistics
29-
print("Completed in {:.3f} seconds of wall-clock time".format(perf_duration))
30-
if len(gdf) > 0:
31-
print("Reference Ground Tracks: {}".format(gdf["rgt"].unique()))
32-
print("Cycles: {}".format(gdf["cycle"].unique()))
33-
print("Received {} elevations".format(len(gdf)))
34-
else:
35-
print("No elevations were returned")
36-
37-
# Display Profile
38-
print("\nTiming Profiles")
39-
for key in icesat2.profiles:
40-
print("{:16}: {:.6f} secs".format(key, icesat2.profiles[key]))
26+
display_statistics(gdf, "elevations")

‎utils/query_photons.py

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
#
2+
# Perform a proxy request for photons
3+
#
4+
5+
import sys
6+
import logging
7+
from sliderule import icesat2
8+
from utils import initialize_client, display_statistics
9+
10+
###############################################################################
11+
# MAIN
12+
###############################################################################
13+
14+
if __name__ == '__main__':
15+
16+
# Configure Logging
17+
logging.basicConfig(level=logging.INFO)
18+
19+
# Initialize Client
20+
parms, cfg = initialize_client(sys.argv)
21+
22+
# Request ATL06 Data
23+
gdf = icesat2.atl03sp(parms, asset=cfg["asset"], resources=[cfg["resource"]])
24+
25+
# Display Statistics
26+
display_statistics(gdf, "photons")

‎utils/region_of_interest.py

Lines changed: 2 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -16,11 +16,10 @@
1616

1717
import sys
1818
import logging
19-
import time
2019
import geopandas as gpd
2120
import matplotlib.pyplot as plt
2221
from sliderule import icesat2
23-
from utils import initialize_client
22+
from utils import initialize_client, display_statistics
2423

2524
###############################################################################
2625
# MAIN
@@ -35,24 +34,10 @@
3534
parms, cfg = initialize_client(sys.argv)
3635

3736
# Get ATL06 Elevations
38-
tstart = time.perf_counter()
3937
atl06 = icesat2.atl06p(parms, cfg["asset"])
40-
perf_duration = time.perf_counter() - tstart
4138

4239
# Display Statistics
43-
print("Completed in {:.3f} seconds of wall-clock time".format(perf_duration))
44-
if len(atl06) > 0:
45-
print("Reference Ground Tracks: {}".format(atl06["rgt"].unique()))
46-
print("Cycles: {}".format(atl06["cycle"].unique()))
47-
print("Received {} elevations".format(len(atl06)))
48-
else:
49-
print("No elevations were returned")
50-
sys.exit()
51-
52-
# Display Profile
53-
print("\nTiming Profiles")
54-
for key in icesat2.profiles:
55-
print("{:16}: {:.6f} secs".format(key, icesat2.profiles[key]))
40+
display_statistics(atl06, "elevations")
5641

5742
# Calculate Extent
5843
lons = [p["lon"] for p in parms["poly"]]

0 commit comments

Comments
 (0)