Skip to content

Commit

Permalink
working atl08 ancillary data
Browse files Browse the repository at this point in the history
  • Loading branch information
@jpswinski
jpswinski committed Dec 1, 2023
1 parent db48027 commit 389d4a6
Show file tree
Hide file tree
Showing 7 changed files with 87 additions and 50 deletions.
2 changes: 1 addition & 1 deletion clients/python/sliderule/icesat2.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,7 @@ def __flattenbatches(rsps, rectype, batch_column, parm, keep_id, as_numpy_array,
if field_name not in field_dictionary:
field_dictionary[field_name] = {'extent_id': [], field_name: []}
field_dictionary[field_name]['extent_id'] += extent_id,
field_dictionary[field_name][field_name] += sliderule.getvalues(field_rec['value'], field_rec['datatype'], len(field_rec['value']))[0],
field_dictionary[field_name][field_name] += sliderule.getvalues(field_rec['value'], field_rec['datatype'], len(field_rec['value']), num_elements=1)[0],
elif 'rsrec' == rsp['__rectype'] or 'zsrec' == rsp['__rectype']:
if rsp["num_samples"] <= 0:
continue
Expand Down
5 changes: 3 additions & 2 deletions clients/python/sliderule/sliderule.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -561,7 +561,7 @@ def procoutputfile(parm):
#
# Get Values from Raw Buffer
#
def getvalues(data, dtype, size):
def getvalues(data, dtype, size, num_elements=0):
"""
data: tuple of bytes
dtype: element of codedtype
Expand All @@ -570,7 +570,8 @@ def getvalues(data, dtype, size):

raw = bytes(data)
datatype = basictypes[codedtype2str[dtype]]["nptype"]
num_elements = int(size / numpy.dtype(datatype).itemsize)
if num_elements == 0: # dynamically determine number of elements
num_elements = int(size / numpy.dtype(datatype).itemsize)
slicesize = num_elements * numpy.dtype(datatype).itemsize # truncates partial bytes
values = numpy.frombuffer(raw[:slicesize], dtype=datatype, count=num_elements)

Expand Down
15 changes: 9 additions & 6 deletions clients/python/tests/test_ancillary.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,21 +54,24 @@ def test_atl06(self, init):
# parms = {
# "poly": region["poly"],
# "srt": icesat2.SRT_LAND,
# "atl08_fields": ["cloud_flag_atm"]
# "atl08_fields": ["h_dif_ref"]
# }
# gdf = icesat2.atl03s(parms, "ATL03_20181017222812_02950102_006_02.h5")
# assert init
# assert len(gdf) == 1029
# assert abs(gdf["cloud_flag_atm"].quantile(q=.75) - 56.234378814697266) < 0.000001
# assert len(gdf) == 403464
# assert abs(gdf["h_dif_ref"].quantile(q=.75) - 10.342773) < 0.000001

def test_atl08_phoreal(self, init):
region = sliderule.toregion(os.path.join(TESTDIR, "data/grandmesa.geojson"))
parms = {
"poly": region["poly"],
"srt": icesat2.SRT_LAND,
"atl08_fields": ["cloud_flag_atm"]
"atl08_fields": ["h_dif_ref", "rgt", "sigma_atlas_land%", "cloud_flag_atm%"]
}
gdf = icesat2.atl08(parms, "ATL03_20181017222812_02950102_006_02.h5")
assert init
assert len(gdf) == 1029
assert abs(gdf["cloud_flag_atm"].quantile(q=.75) - 56.234378814697266) < 0.000001
assert len(gdf) == 819
assert abs(gdf["h_dif_ref"].quantile(q=.50) - -0.9146728515625) < 0.000001
assert abs(gdf["rgt_y"].quantile(q=.50) - 295.0) < 0.000001
assert abs(gdf["sigma_atlas_land%"].quantile(q=.50) - 0.2402431309223175) < 0.000001
assert abs(gdf["cloud_flag_atm%"].quantile(q=.50) - 1.0) < 0.000001
62 changes: 42 additions & 20 deletions plugins/icesat2/plugin/Atl03Reader.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -580,6 +580,7 @@ Atl03Reader::Atl03Data::~Atl03Data (void)
Atl03Reader::Atl08Class::Atl08Class (info_t* info):
enabled (info->reader->parms->stages[Icesat2Parms::STAGE_ATL08]),
phoreal (info->reader->parms->stages[Icesat2Parms::STAGE_PHOREAL]),
ancillary (info->reader->parms->atl08_fields != NULL),
classification {NULL},
relief {NULL},
landcover {NULL},
Expand All @@ -588,20 +589,19 @@ Atl03Reader::Atl08Class::Atl08Class (info_t* info):
atl08_pc_indx (enabled ? info->reader->asset : NULL, info->reader->resource08, FString("%s/%s", info->prefix, "signal_photons/classed_pc_indx").c_str(), &info->reader->context08),
atl08_pc_flag (enabled ? info->reader->asset : NULL, info->reader->resource08, FString("%s/%s", info->prefix, "signal_photons/classed_pc_flag").c_str(), &info->reader->context08),
atl08_ph_h (phoreal ? info->reader->asset : NULL, info->reader->resource08, FString("%s/%s", info->prefix, "signal_photons/ph_h").c_str(), &info->reader->context08),
segment_id_beg (phoreal ? info->reader->asset : NULL, info->reader->resource08, FString("%s/%s", info->prefix, "land_segments/segment_id_beg").c_str(), &info->reader->context08),
segment_id_beg ((phoreal || ancillary) ? info->reader->asset : NULL, info->reader->resource08, FString("%s/%s", info->prefix, "land_segments/segment_id_beg").c_str(), &info->reader->context08),
segment_landcover (phoreal ? info->reader->asset : NULL, info->reader->resource08, FString("%s/%s", info->prefix, "land_segments/segment_landcover").c_str(), &info->reader->context08),
segment_snowcover (phoreal ? info->reader->asset : NULL, info->reader->resource08, FString("%s/%s", info->prefix, "land_segments/segment_snowcover").c_str(), &info->reader->context08),
anc_seg_data (NULL),
anc_seg_indices (NULL)
{
AncillaryFields::list_t* atl08_fields = info->reader->parms->atl08_fields;

if(atl08_fields)
if(ancillary)
{
/* Allocate Ancillary Data Dictionary */
anc_seg_data = new H5DArrayDictionary(Icesat2Parms::EXPECTED_NUM_FIELDS);
anc_seg_indices = new List<int32_t>;

/* Read Ancillary Fields */
AncillaryFields::list_t* atl08_fields = info->reader->parms->atl08_fields;
for(int i = 0; i < atl08_fields->length(); i++)
{
const char* field_name = (*atl08_fields)[i].field.c_str();
Expand Down Expand Up @@ -633,7 +633,7 @@ Atl03Reader::Atl08Class::~Atl08Class (void)
delete [] landcover;
delete [] snowcover;
delete anc_seg_data;
delete anc_seg_indices;
delete [] anc_seg_indices;
}

/*----------------------------------------------------------------------------
Expand All @@ -651,10 +651,13 @@ void Atl03Reader::Atl08Class::classify (info_t* info, const Region& region, cons
atl08_segment_id.join(info->reader->read_timeout_ms, true);
atl08_pc_indx.join(info->reader->read_timeout_ms, true);
atl08_pc_flag.join(info->reader->read_timeout_ms, true);
if(phoreal || ancillary)
{
segment_id_beg.join(info->reader->read_timeout_ms, true);
}
if(phoreal)
{
atl08_ph_h.join(info->reader->read_timeout_ms, true);
segment_id_beg.join(info->reader->read_timeout_ms, true);
segment_landcover.join(info->reader->read_timeout_ms, true);
segment_snowcover.join(info->reader->read_timeout_ms, true);
}
Expand All @@ -671,6 +674,11 @@ void Atl03Reader::Atl08Class::classify (info_t* info, const Region& region, cons
snowcover = new uint8_t [num_photons];
}

if(ancillary)
{
anc_seg_indices = new int32_t [num_photons];
}

/* Populate ATL08 Classifications */
int32_t atl03_photon = 0;
int32_t atl08_photon = 0;
Expand All @@ -680,7 +688,7 @@ void Atl03Reader::Atl08Class::classify (info_t* info, const Region& region, cons
int32_t atl03_segment = atl03.segment_id[atl03_segment_index];

/* Get Land and Snow Flags */
if(phoreal)
if(phoreal || ancillary)
{
while( (atl08_segment_index < segment_id_beg.size) &&
((segment_id_beg[atl08_segment_index] + NUM_ATL03_SEGS_IN_ATL08_SEG) <= atl03_segment) )
Expand All @@ -707,12 +715,6 @@ void Atl03Reader::Atl08Class::classify (info_t* info, const Region& region, cons
atl08_photon++;
}

/* Populate Ancillary Indices */
if(anc_seg_indices)
{
anc_seg_indices->add(atl08_segment_index);
}

/* Check Match */
if( (atl08_photon < atl08_segment_id.size) &&
(atl08_segment_id[atl08_photon] == atl03_segment) &&
Expand Down Expand Up @@ -744,6 +746,12 @@ void Atl03Reader::Atl08Class::classify (info_t* info, const Region& region, cons
}
}

/* Populate Ancillary Index */
if(ancillary)
{
anc_seg_indices[atl03_photon] = atl08_segment_index;
}

/* Go To Next ATL08 Photon */
atl08_photon++;
}
Expand All @@ -760,11 +768,10 @@ void Atl03Reader::Atl08Class::classify (info_t* info, const Region& region, cons
snowcover[atl03_photon] = INVALID_FLAG;
}

/* Set Ancillary Indices to Invalid */
if(anc_seg_indices)
/* Set Ancillary Index to Invalid */
if(ancillary)
{
int32_t tmp_invalid = Atl03Reader::INVALID_INDICE;
anc_seg_indices->add(tmp_invalid);
anc_seg_indices[atl03_photon] = Atl03Reader::INVALID_INDICE;
}
}

Expand Down Expand Up @@ -1164,6 +1171,7 @@ void* Atl03Reader::subsettingThread (void* parm)
stats_t local_stats = {0, 0, 0, 0, 0};
List<int32_t>* segment_indices = NULL; // used for ancillary data
List<int32_t>* photon_indices = NULL; // used for ancillary data
List<int32_t>* atl08_indices = NULL; // used for ancillary data

/* Start Trace */
uint32_t trace_id = start_trace(INFO, reader->traceId, "atl03_subsetter", "{\"asset\":\"%s\", \"resource\":\"%s\", \"track\":%d}", info->reader->asset->getName(), info->reader->resource, info->track);
Expand Down Expand Up @@ -1229,6 +1237,13 @@ void* Atl03Reader::subsettingThread (void* parm)
else photon_indices = new List<int32_t>;
}

/* Ancillary ATL08 Fields */
if(atl08.anc_seg_data)
{
if(atl08_indices) atl08_indices->clear();
else atl08_indices = new List<int32_t>;
}

/* Traverse Photons Until Desired Along Track Distance Reached */
while(!extent_complete || !step_complete)
{
Expand Down Expand Up @@ -1372,12 +1387,18 @@ void* Atl03Reader::subsettingThread (void* parm)
{
segment_indices->add(current_segment);
}
/* Index Photon for Ancillary Fields */

/* Index Photon for Ancillary Fields */
if(photon_indices)
{
photon_indices->add(current_photon);
}

/* Index ATL08 Segment for Photon for Ancillary Fields */
if(atl08_indices)
{
atl08_indices->add(atl08.anc_seg_indices[current_photon]);
}
} while(false);
}
else
Expand Down Expand Up @@ -1453,7 +1474,7 @@ void* Atl03Reader::subsettingThread (void* parm)
reader->generateExtentRecord(extent_id, info, state, atl03, rec_list, rec_total_size);
Atl03Reader::generateAncillaryRecords(extent_id, parms->atl03_ph_fields, atl03.anc_ph_data, AncillaryFields::PHOTON_ANC_TYPE, photon_indices, rec_list, rec_total_size);
Atl03Reader::generateAncillaryRecords(extent_id, parms->atl03_geo_fields, atl03.anc_geo_data, AncillaryFields::EXTENT_ANC_TYPE, segment_indices, rec_list, rec_total_size);
Atl03Reader::generateAncillaryRecords(extent_id, parms->atl08_fields, atl08.anc_seg_data, AncillaryFields::ATL08_ANC_TYPE, atl08.anc_seg_indices, rec_list, rec_total_size);
Atl03Reader::generateAncillaryRecords(extent_id, parms->atl08_fields, atl08.anc_seg_data, AncillaryFields::ATL08_ANC_TYPE, atl08_indices, rec_list, rec_total_size);

/* Send Records */
if(rec_list.size() == 1)
Expand Down Expand Up @@ -1524,6 +1545,7 @@ void* Atl03Reader::subsettingThread (void* parm)
/* Clean Up Indices */
delete segment_indices;
delete photon_indices;
delete atl08_indices;

/* Clean Up Info */
delete info;
Expand Down
3 changes: 2 additions & 1 deletion plugins/icesat2/plugin/Atl03Reader.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,6 +218,7 @@ class Atl03Reader: public LuaObject
/* Class Data */
bool enabled;
bool phoreal;
bool ancillary;

/* Generated Data */
uint8_t* classification; // [num_photons]
Expand All @@ -238,7 +239,7 @@ class Atl03Reader: public LuaObject

/* Ancillary Data */
H5DArrayDictionary* anc_seg_data;
List<int32_t>* anc_seg_indices;
int32_t* anc_seg_indices;
};

/* YAPC Score Subclass */
Expand Down
41 changes: 26 additions & 15 deletions plugins/icesat2/plugin/Atl08Dispatch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,10 +36,13 @@
#include <math.h>
#include <float.h>
#include <map>
#include <limits>

#include "core.h"
#include "icesat2.h"

using std::numeric_limits;

/******************************************************************************
* STATIC DATA
******************************************************************************/
Expand Down Expand Up @@ -190,9 +193,6 @@ bool Atl08Dispatch::processRecord (RecordObject* record, okey_t key, recVec_t* r
return true;
}

/* Build Ancillary Inputs */
RecordObject* atl08_anc_rec = buildAncillaryRecord(extent, records);

/* Initialize Results */
vegetation_t result;
result.pflags = 0;
Expand All @@ -204,6 +204,9 @@ bool Atl08Dispatch::processRecord (RecordObject* record, okey_t key, recVec_t* r
phorealAlgorithm(extent, result);
}

/* Build Ancillary Inputs */
RecordObject* atl08_anc_rec = buildAncillaryRecord(extent, records);

/* Post Results */
postResult(&result, atl08_anc_rec);

Expand Down Expand Up @@ -254,50 +257,58 @@ RecordObject* Atl08Dispatch::buildAncillaryRecord (Atl03Reader::extent_t* extent
AncillaryFields::field_t field;
field.anc_type = atl03_anc_rec->anc_type;
field.field_index = atl03_anc_rec->field_index;
field.data_type = atl03_anc_rec->data_type;

/* Calculate Estimation */
AncillaryFields::setValueAsInteger(&field, 0.0);
if((atl03_anc_rec->data_type == (uint8_t)RecordObject::DOUBLE) || (atl03_anc_rec->data_type == (uint8_t)RecordObject::FLOAT))
{
AncillaryFields::setValueAsDouble(&field, 0.0);
double* values = AncillaryFields::extractAsDoubles(atl03_anc_rec); // `new` memory allocated here
if(entry.estimation == AncillaryFields::NEAREST_NEIGHBOR)
{
std::map<double, int> counts;
for(unsigned int j = 0; j < atl03_anc_rec->num_elements; j++)
{
if(counts.count(values[j])) counts[values[j]] += 1;
else counts[values[j]] = 1;
if(values[j] < numeric_limits<float>::max())
{
if(counts.count(values[j])) counts[values[j]] += 1;
else counts[values[j]] = 1;
}
}
double nearest = 0.0;
int nearest_count = 0;
for (auto itr = counts.begin(); itr != counts.end(); ++itr)
for (auto itr: counts)
{
if(itr->second > nearest_count)
if(itr.second > nearest_count)
{
nearest_count = itr->second;
nearest = itr->first;
nearest_count = itr.second;
nearest = itr.first;
}
}
AncillaryFields::setValueAsDouble(&field, nearest);
}
else if(entry.estimation == AncillaryFields::INTERPOLATION)
{
double average = 0.0;
int samples = 0;
for(unsigned int j = 0; j < atl03_anc_rec->num_elements; j++)
{
average += values[j];
if(values[j] < numeric_limits<float>::max())
{
average += values[j];
samples++;
}
}
if(atl03_anc_rec->num_elements)
if(samples)
{
average /= atl03_anc_rec->num_elements;
average /= samples;
}
AncillaryFields::setValueAsDouble(&field, average);
}
delete [] values;
}
else // integer type
{
AncillaryFields::setValueAsInteger(&field, 0);
int64_t* values = AncillaryFields::extractAsIntegers(atl03_anc_rec); // `new` memory allocated here
if(entry.estimation == AncillaryFields::NEAREST_NEIGHBOR)
{
Expand Down Expand Up @@ -340,7 +351,7 @@ RecordObject* Atl08Dispatch::buildAncillaryRecord (Atl03Reader::extent_t* extent
}

/* Return Ancillary Record */
return AncillaryFields::createFieldArrayRecord(extent->extent_id, field_vec);
return AncillaryFields::createFieldArrayRecord(extent->extent_id | Icesat2Parms::EXTENT_ID_ELEVATION, field_vec);
}

/*----------------------------------------------------------------------------
Expand Down
Loading

0 comments on commit 389d4a6

Please sign in to comment.