NASA Ocean Color

ocssw V2022

 #include "get_geospatial.hpp"
 #include "find_variable.hpp"
 #include "boost/algorithm/string.hpp"
 #include <boost/geometry.hpp>
 #include <boost/geometry/geometries/point_xy.hpp>
 #include <boost/geometry/geometries/polygon.hpp>
 #include <algorithm>
 #include <stack>
 namespace bg = boost::geometry;
 typedef bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree>> Point_t;
 typedef bg::model::linestring<Point_t> Linestring_t;
 typedef bg::model::polygon<Point_t> Polygon_t;
 typedef bg::model::multi_point<Point_t> MultiPoint_t;
  
 Geospatialbounds::Geospatialbounds() {
     scan_line_vars["elat"] = std::vector<float>();
     scan_line_vars["slat"] = std::vector<float>();
     scan_line_vars["clat"] = std::vector<float>();
     scan_line_vars["elon"] = std::vector<float>();
     scan_line_vars["slon"] = std::vector<float>();
 }
  
 std::string Geospatialbounds::get_day_night_flag(const netCDF::NcFile &nc_file) {
     netCDF::NcVar solzen_nc;
     for (const auto &name: possible_sozlen_names) {
         solzen_nc = find_nc_variable_cpp(name, nc_file);
         if (!solzen_nc.isNull()) break;
     }
     if (solzen_nc.isNull()) {
         std::cerr << "--Error--: Could not find solar zenith angle to compute Day Night Flag" << std::endl;
         exit(EXIT_FAILURE);
     }
     std::vector<float> solzen;
     allocate_var(solzen, solzen_nc);
     const float night_angle = 90.0f;
     std::stack<bool> state;
     for (const auto &angle: solzen) {
         if (angle > 180.0 || angle < -180 || std::isnan(angle) || !std::isfinite(angle))
             continue;
         bool flag = angle < night_angle;
         if (state.empty()) {
             state.push(flag);
         } else if (state.top() != flag) {
             return "Mixed";
         }
     }
     if (state.empty()) {
         return "Undefined";
     } else {
         if (state.top())
             return "Day";
         else
             return "Night";
     }
 }
  
 void Geospatialbounds::get_lat_lon(const netCDF::NcFile &nc_input) {
     netCDF::NcVar latnc = find_nc_variable_cpp("latitude", nc_input);
     netCDF::NcVar lonnc = find_nc_variable_cpp("longitude", nc_input);
     allocate_var(lat, latnc);
     allocate_var(lon, lonnc);
     std::vector<netCDF::NcDim> dims = latnc.getDims();
     if (dims.size() != 2) {
         std::cerr << "-Error- : Latitude array must be two-dimensional." << std::endl;
     } else {
         number_of_lines = dims.at(0).getSize();
         pixels_per_line = dims.at(1).getSize();
     }
 }
  
 void Geospatialbounds::allocate_var(std::vector<float> &data, const netCDF::NcVar &var) {
     if (!var.isNull()) {
         std::vector<netCDF::NcDim> dims = var.getDims();
         size_t size = 1;
         for (const auto &dim: dims) {
             size *= dim.getSize();
         }
         data.resize(size);
         var.getVar(data.data());
     }
 }
  
 void Geospatialbounds::allocate_attr(std::vector<float> &data, const netCDF::NcAtt &att) {
     if (!att.isNull()) {
         size_t size = att.getAttLength();
         if (size == 0)
             return;
         if (att.getType() == NC_STRING || att.getType() == NC_CHAR)
             return;
         data.resize(size);
         att.getValues(data.data());
     }
 }
  
 Geospatialbounds::Geospatialbounds(const std::string &path)
         : Geospatialbounds(netCDF::NcFile(path, netCDF::NcFile::read)) {
 }
  
 Geospatialbounds::Geospatialbounds(const netCDF::NcFile &nc_input) : Geospatialbounds() {
     netCDF::NcGroupAtt attr = nc_input.getAtt("instrument");
     if (attr.isNull()) {
         std::cerr << "Warning: Instrument not found. Returning an empty string" << std::endl;
     } else {
         attr.getValues(instrument);
     }
     attr = nc_input.getAtt("platform");
     if (attr.isNull()) {
         std::cerr << "Warning: Platform not found. Returning an empty string" << std::endl;
     } else {
         attr.getValues(platform);
     }
     // search for gring lat/lon or polygon in global attributes
  
     const auto attrs = nc_input.getAtts();
     for (const auto &_attr: attrs) {
         std::string attr_name = _attr.first;
         boost::algorithm::to_lower(attr_name);
         if (attr_name == "gringpointlatitude") {
             allocate_attr(gringpointlatitude, _attr.second);
         }
         if (attr_name == "gringpointlongitude") {
             allocate_attr(gringpointlongitude, _attr.second);
         }
         if (attr_name == "geospatial_bounds") {
             _attr.second.getValues(geospatial_bounds_wkt);
         }
     }
     // search for group attributes
     if (gringpointlatitude.empty()) {
         auto gringpointlatitude_attr = find_nc_grp_attribute_cpp("gringpointlatitude", nc_input);
         allocate_attr(gringpointlatitude, gringpointlatitude_attr);
     }
     if (gringpointlongitude.empty()) {
         auto gringpointlongitude_attr = find_nc_grp_attribute_cpp("gringpointlongitude", nc_input);
         allocate_attr(gringpointlongitude, gringpointlongitude_attr);
     }
     if (geospatial_bounds_wkt.empty()) {
         auto geospatial_bounds_wkt_attr = find_nc_grp_attribute_cpp("geospatial_bounds", nc_input);
         if (!geospatial_bounds_wkt_attr.isNull())
             geospatial_bounds_wkt_attr.getValues(geospatial_bounds_wkt);
     }
     auto crs = find_nc_grp_attribute_cpp("geospatial_bounds_crs", nc_input);
     if (!crs.isNull()) {
         Polygon_t poly;
         boost::geometry::read_wkt(geospatial_bounds_wkt, poly);
         geospatial_bounds_wkt = "POLYGON((";
         double first_ini = true;
         double x_first, y_first;
         // getting the vertices back
         for (auto it = boost::begin(boost::geometry::exterior_ring(poly));
              it != boost::end(boost::geometry::exterior_ring(poly)); ++it) {
             double x = bg::get<0>(*it);
             double y = bg::get<1>(*it);
             if (first_ini) {
                 x_first = x;
                 y_first = y;
                 first_ini = false;
             }
             geospatial_bounds_wkt += std::to_string(y);
             geospatial_bounds_wkt += " ";
             geospatial_bounds_wkt += std::to_string(x);
             geospatial_bounds_wkt += ", ";
             // use the coordinates...
         }
         // Polygon must be closed though it does not affect dateline crossing
         geospatial_bounds_wkt += std::to_string(y_first);
         geospatial_bounds_wkt += " ";
         geospatial_bounds_wkt += std::to_string(x_first);
         geospatial_bounds_wkt += "))";
     }
     if (std::isnan(max_lon)) {
         auto geospatial_lon_max_attr = find_nc_grp_attribute_cpp("geospatial_lon_max", nc_input);
         if (!geospatial_lon_max_attr.isNull())
             geospatial_lon_max_attr.getValues(&max_lon);
     }
     if (std::isnan(min_lon)) {
         auto geospatial_lon_min_attr = find_nc_grp_attribute_cpp("geospatial_lon_min", nc_input);
         if (!geospatial_lon_min_attr.isNull())
             geospatial_lon_min_attr.getValues(&min_lon);
     }
     if (time_coverage_start.empty()) {
         auto attr_time_coverage = find_nc_grp_attribute_cpp("time_coverage_start", nc_input);
         if (!attr_time_coverage.isNull())
             attr_time_coverage.getValues(time_coverage_start);
     }
     if (day_night_flag.empty()) {
         auto attr_day_night_flag = find_nc_grp_attribute_cpp("day_night_flag", nc_input);
         if (!attr_day_night_flag.isNull())
             attr_day_night_flag.getValues(day_night_flag);
         else
             day_night_flag = get_day_night_flag(nc_input);
     }
     for (auto &sc_ln_attr: scan_line_vars) {
         std::string var_name = sc_ln_attr.first;
         auto &var_val = sc_ln_attr.second;
         netCDF::NcVar var = find_nc_variable_cpp(var_name, nc_input);
         allocate_var(var_val, var);
     }
     for (const auto &sc_ln_attr: scan_line_vars) {
         auto &var_val = sc_ln_attr.second;
         if (var_val.empty()) {
             get_lat_lon(nc_input);
             break;
         }
         number_of_lines = var_val.size();
     }
 }
  
 void Geospatialbounds::calc_scan_line(const std::vector<float> &lat, const std::vector<float> &lon,
                                       size_t numScans, size_t numPixels,
                                       std::unordered_map<std::string, std::vector<float>> &scan_line_vars) {
     assert(lat.size() == numPixels * numScans);
     assert(lon.size() == numPixels * numScans);
     scan_line_vars.at("elat") = std::vector<float>(numScans, NAN);
     scan_line_vars.at("elon") = std::vector<float>(numScans, NAN);
     scan_line_vars.at("slon") = std::vector<float>(numScans, NAN);
     scan_line_vars.at("slat") = std::vector<float>(numScans, NAN);
     scan_line_vars.at("clat") = std::vector<float>(numScans, NAN);
     std::vector<int8_t> mask(lat.size(), 1);
     auto f_lat = [](const float &v) {
         if (v > 90.0 || v < -90 || std::isnan(v) || !std::isfinite(v))
             return 0;
         else
             return 1;
     };
     auto f_lon = [](const float &v, const int8_t &mask) {
         if (v > 180.0 || v < -180 || std::isnan(v) || !std::isfinite(v) || !mask)
             return 0;
         else
             return 1;
     };
     std::transform(lat.begin(), lat.end(), mask.begin(), f_lat);
     std::transform(lat.begin(), lat.end(), mask.begin(), mask.begin(), f_lon);
     auto get_index = [&](size_t is, size_t ip) { return numPixels * is + ip; };
     float last_lat = NAN;
     for (size_t is = 0; is < numScans; is++) {
         size_t spix = 0;
         size_t cpix = numPixels / 2;
         size_t epix = numPixels - 1;
         // find good geolocation for spix and epix
         size_t good_spix = spix;
         size_t good_epix = epix;
  
         for (size_t ip = spix; ip < epix; ip++) {
             size_t i = get_index(is, ip);
             if (mask.at(i)) {
                 good_spix = ip;
                 break;
             }
         }
         for (size_t ip = epix; ip >= spix; ip--) {
             size_t i = get_index(is, ip);
             if (mask.at(i)) {
                 good_epix = ip;
                 break;
             }
         }
         size_t i = get_index(is, good_epix);
         scan_line_vars.at("elat")[is] = lat.at(i);
         scan_line_vars.at("elon")[is] = lon.at(i);
         i = get_index(is, good_spix);
         scan_line_vars.at("slat")[is] = lat.at(i);
         scan_line_vars.at("slon")[is] = lon.at(i);
         if (std::isnan(last_lat))
             last_lat = (scan_line_vars.at("elat")[is] + scan_line_vars.at("slat")[is]) / 2;
         i = get_index(is, cpix);
         if (f_lat(lat.at(i))) {
             last_lat = lat.at(i);
         }
         scan_line_vars.at("clat")[is] = last_lat;
     }
 };
  
 std::pair<std::vector<float>, std::vector<float>> Geospatialbounds::calc_gring(
         const std::vector<float> &slat, const std::vector<float> &elat, const std::vector<float> &clat,
         const std::vector<float> &slon, const std::vector<float> &elon) {
     const float geobox_bounds{20.0f};
     std::vector<std::vector<float>> geobox{{},
                                            {},
                                            {},
                                            {}};
     // find last valid scan
     size_t nscan = clat.size();
     // find first valid scan
     size_t first_scan = nscan;
     size_t last_scan = nscan;
     for (size_t is = 0; is < nscan; is++) {
         if (!std::isnan(elat.at(is))) {
             if (first_scan == nscan)
                 first_scan = is;
             last_scan = is;
         }
     }
     if (first_scan == nscan) {
         std::cerr << "-Error- : not valid elat is supplied. Exiting";
         exit(EXIT_FAILURE);
     }
     geobox.at(0).push_back(slon.at(first_scan));
     geobox.at(1).push_back(slat.at(first_scan));
     geobox.at(2).push_back(elon.at(first_scan));
     geobox.at(3).push_back(elat.at(first_scan));
     //
     float last_lat = clat.at(first_scan);
     for (size_t is = first_scan + 1; is < last_scan; is++) {
         if (std::isnan(clat.at(is)))
             continue;
         if (std::isnan(elat.at(is)))
             continue;
         if (std::abs(last_lat - clat.at(is)) > geobox_bounds) {
             geobox.at(0).push_back(slon.at(is));
             geobox.at(1).push_back(slat.at(is));
             geobox.at(2).push_back(elon.at(is));
             geobox.at(3).push_back(elat.at(is));
             last_lat = clat.at(is);
         }
     }
     geobox.at(0).push_back(slon.at(last_scan));
     geobox.at(1).push_back(slat.at(last_scan));
     geobox.at(2).push_back(elon.at(last_scan));
     geobox.at(3).push_back(elat.at(last_scan));
     std::vector<float> lat_gring;
     std::vector<float> lon_gring;
     lat_gring.push_back(geobox.at(1).at(0));
     for (const auto &val: geobox.at(3)) {
         lat_gring.push_back(val);
     }
     while (geobox.at(1).size() > 1) {
         lat_gring.push_back(geobox.at(1).back());
         geobox.at(1).pop_back();
     }
  
     lon_gring.push_back(geobox.at(0).at(0));
     for (const auto &val: geobox.at(2)) {
         lon_gring.push_back(val);
     }
     while (geobox.at(0).size() > 1) {
         lon_gring.push_back(geobox.at(0).back());
         geobox.at(0).pop_back();
     }
     return {lat_gring, lon_gring};
 }
  
 std::string Geospatialbounds::calc_gring(const std::pair<std::vector<float>, std::vector<float>> &gring) {
     const std::vector<float> &lat_gring = gring.first;
     const std::vector<float> &lon_gring = gring.second;
     std::string wkt = "POLYGON((";
     size_t len = lat_gring.size();
     for (size_t i = 0; i < len; i++) {
         wkt += std::to_string(lon_gring.at(i));
         wkt += " ";
         wkt += std::to_string(lat_gring.at(i));
         wkt += ", ";
     }
     // Polygon must be closed though it does not affect dateline crossing
     wkt += std::to_string(lon_gring.at(0));
     wkt += " ";
     wkt += std::to_string(lat_gring.at(0));
     wkt += "))";
     return wkt;
 };
  
 const std::unordered_map<std::string, std::vector<float>> &Geospatialbounds::get_bounds() {
     if (scan_line_vars.at("elat").empty()) {
         calc_scan_line(lat, lon, number_of_lines, pixels_per_line, scan_line_vars);
     }
     return scan_line_vars;
 }
  
 std::pair<std::vector<float>, std::vector<float>> Geospatialbounds::get_gring() {
     if (gringpointlatitude.empty()) {
         scan_line_vars = get_bounds();
         const auto &slat = scan_line_vars.at("slat");
         const auto &elat = scan_line_vars.at("elat");
         const auto &clat = scan_line_vars.at("clat");
         const auto &slon = scan_line_vars.at("slon");
         const auto &elon = scan_line_vars.at("elon");
         auto res = calc_gring(slat, elat, clat, slon, elon);
         gringpointlatitude = res.first;
         gringpointlongitude = res.second;
     }
     return {gringpointlatitude, gringpointlongitude};
 }
  
 std::string Geospatialbounds::get_bounds_wkt() {
     if (geospatial_bounds_wkt.empty()) {
         geospatial_bounds_wkt = calc_gring(get_gring());
     }
     return geospatial_bounds_wkt;
 }
  
 float Geospatialbounds::get_max_lon() {
     if (std::isnan(max_lon)) {
         auto res = get_gring();
         const auto &longring = res.second;
         max_lon = *std::max_element(longring.begin(), longring.end());
         min_lon = *std::min_element(longring.begin(), longring.end());
     }
     return max_lon;
 }
  
 float Geospatialbounds::get_min_lon() {
     get_max_lon();
     return min_lon;
 }
  
 const float *Geospatialbounds::get_elat() {
     if (scan_line_vars.at("elat").empty())
         scan_line_vars = get_bounds();
     return scan_line_vars.at("elat").data();
 }
  
 const float *Geospatialbounds::get_slat() {
     if (scan_line_vars.at("slat").empty())
         scan_line_vars = get_bounds();
     return scan_line_vars.at("slat").data();
 }

seadasutils.DictUtils.val

int val

Definition: DictUtils.py:140

convert2ocrvc.instrument

instrument

Definition: convert2ocrvc.py:29

nc_file

Definition: val_extract.h:241

Geospatialbounds::calc_scan_line

static void calc_scan_line(const std::vector< float > &, const std::vector< float > &, size_t, size_t, std::unordered_map< std::string, std::vector< float >> &)

Definition: get_geospatial.cpp:209

mapgen_overlay.transform

transform

Definition: mapgen_overlay.py:186

MDN.setup.name

name

Definition: setup.py:14

hico.value_locate.res

def res

Definition: value_locate.py:66

MultiPoint_t

bg::model::multi_point< Point_t > MultiPoint_t

Definition: get_geospatial.cpp:13

find_nc_grp_attribute_cpp

netCDF::NcGroupAtt find_nc_grp_attribute_cpp(const std::string &att_name, T &nc_file)

Definition: find_variable.hpp:22

color_dtdb.v

Definition: color_dtdb.py:430

Geospatialbounds::allocate_var

static void allocate_var(std::vector< float > &, const netCDF::NcVar &var)

Definition: get_geospatial.cpp:71

Geospatialbounds::get_bounds_wkt

std::string get_bounds_wkt()

Definition: get_geospatial.cpp:384

Polygon_t

bg::model::polygon< Point_t > Polygon_t

Definition: get_geospatial.cpp:12

this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned resolving bug report Corrected to filter out Aqua attitude records with missing status helping resolve bug MOD_PR03 will still correctly write scan and pixel data that does not depend upon the start thereby resolving MODur00108 Changed header guard macro names to avoid reserved name resolving MODur00104 Maneuver list file for Terra satellite was updated to include data from Jecue DuChateu Maneuver list files for both Terra and Aqua were updated to include two maneuvers from recent IOT weekly reports The limits for Z component of angular momentum was and to set the fourth GEO scan quality flag for a scan depending upon whether or not it occurred during one of them Added _FillValue attributes to many and changed the fill value for the sector start times from resolving MODur00072 Writes boundingcoordinate metadata to L1A archived metadata For PERCENT *ECS change to treat rather resolving GSFcd01518 Added a LogReport Changed to create the Average Temperatures vdata even if the number of scans is

Definition: HISTORY.txt:301

nscan

int32 nscan

Definition: l1_czcs_hdf.c:19

string

@ string

Definition: GEO_read_param_file.c:64

find_nc_variable_cpp

netCDF::NcVar find_nc_variable_cpp(const std::string &var_name, T &nc_file)

Definition: find_variable.hpp:7

Point_t

bg::model::point< double, 2, bg::cs::spherical_equatorial< bg::degree > > Point_t

Definition: get_geospatial.cpp:10

Geospatialbounds::Geospatialbounds

Geospatialbounds()

Definition: get_geospatial.cpp:15

find_variable.hpp

sort_gring.gringpointlatitude

gringpointlatitude

Definition: sort_gring.py:36

Geospatialbounds::get_max_lon

float get_max_lon()

Definition: get_geospatial.cpp:391

Geospatialbounds::get_elat

const float * get_elat()

Definition: get_geospatial.cpp:406

run_local.lat

lat

Definition: run_local.py:69

state

int state(double tjdTDB, JPLIntUtilType *util, double posvel[13][6], double *pnut)

run_local.lon

lon

Definition: run_local.py:69

Geospatialbounds

Definition: get_geospatial.hpp:11

sort_gring.gringpointlongitude

gringpointlongitude

Definition: sort_gring.py:35

color_dtdb.path

string path

Definition: color_dtdb.py:221

hico.value_locate.x

list x

Definition: value_locate.py:64

color_dtdb.y

Definition: color_dtdb.py:430

hktgen_pace.att

int att

Definition: hktgen_pace.py:28

Polygon_t

bg::model::polygon< Point_t > Polygon_t

Definition: get_dataday.cpp:25

mask

a context in which it is NOT documented to do so subscript which cannot be easily calculated when extracting TONS attitude data from the Terra L0 files Corrected several defects in extraction of entrained ephemeris and and as HDF file for both the L1A and Geolocation enabling retrieval of South Polar DEM data Resolved Bug by changing to opent the geolocation file only after a successful read of the L1A and also by checking for fatal errors from not restoring C5 and to report how many of those high resolution values were water in the new WaterPresent SDS Added valid_range attribute to Land SeaMask Changed to bilinearly interpolate the geoid_height to remove artifacts at one degree lines Made corrections to const qualification of pointers allowed by new version of M API library Removed casts that are no longer for same not the geoid Corrected off by one error in calculation of high resolution offsets Corrected parsing of maneuver list configuration parameter Corrected to set Height SDS to fill values when geolocation when for elevation and land water mask

Definition: HISTORY.txt:114

Geospatialbounds::get_day_night_flag

const std::string & get_day_night_flag() const

Definition: get_geospatial.hpp:62

get_geospatial.hpp

data

no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data

Definition: HISTORY.txt:356

convert2ocrvc.platform

platform

Definition: convert2ocrvc.py:30

Linestring_t

bg::model::linestring< Point_t > Linestring_t

Definition: get_geospatial.cpp:11

color_dtdb.bg

Definition: color_dtdb.py:451

geometry

subroutine geometry

Definition: atrem_app_refl_f90_cubeio.f:1329

Geospatialbounds::get_bounds

const std::unordered_map< std::string, std::vector< float > > & get_bounds()

Definition: get_geospatial.cpp:362

ncattredit.attrs

attrs

Definition: ncattredit.py:34

get_index

int16 get_index(char *Vdata_name)

Definition: get_index.c:6

Geospatialbounds::get_slat

const float * get_slat()

Definition: get_geospatial.cpp:412

spix

int32 spix

Definition: l1_czcs_hdf.c:21

georegion_gen.wkt

wkt

Definition: georegion_gen.py:149

max_lon

double max_lon

Definition: GEO_DEM.c:64

mapgen_overlay.size

size

Definition: mapgen_overlay.py:242

Geospatialbounds::allocate_attr

static void allocate_attr(std::vector< float > &, const netCDF::NcAtt &att)

Definition: get_geospatial.cpp:83

epix

int32 epix

Definition: l1_czcs_hdf.c:23

Geospatialbounds::get_gring

std::pair< std::vector< float >, std::vector< float > > get_gring()

Definition: get_geospatial.cpp:369

check_products.end

end

Definition: check_products.py:189

abs

#define abs(a)

Definition: misc.h:90

int i

Definition: decode_rs.h:71

Geospatialbounds::get_min_lon

float get_min_lon()

Definition: get_geospatial.cpp:401

Geospatialbounds::calc_gring

static std::pair< std::vector< float >, std::vector< float > > calc_gring(const std::vector< float > &, const std::vector< float > &, const std::vector< float > &, const std::vector< float > &, const std::vector< float > &)

Definition: get_geospatial.cpp:274

mapgen_overlay.crs

def crs

Definition: mapgen_overlay.py:176