Program Listing for File GeoidUndulationSource.hpp

↰ Return to documentation for file (src/navtk/geospatial/sources/GeoidUndulationSource.hpp)

#pragma once

#include <fstream>
#include <memory>

#include <navtk/factory.hpp>
#include <navtk/geospatial/sources/SpatialMapDataSource.hpp>
#include <navtk/tensors.hpp>

namespace navtk {
namespace geospatial {

class GeoidUndulationSource : public SpatialMapDataSource {
public:
    static std::shared_ptr<GeoidUndulationSource> get_shared(
        const std::string& path = "WW15MGH.GRD");

    std::pair<bool, double> lookup_datum(double latitude, double longitude) const override;

    void set_chunk_size(Size size);

    Size get_chunk_size() const;

private:
    // private constructor because global instance is implemented via `get_shared()`.
    GeoidUndulationSource(const std::string& path);
    GeoidUndulationSource(const GeoidUndulationSource&)             = delete;
    GeoidUndulationSource& operator=(const GeoidUndulationSource&)  = delete;
    GeoidUndulationSource(const GeoidUndulationSource&&)            = delete;
    GeoidUndulationSource& operator=(const GeoidUndulationSource&&) = delete;

    // Check if either latitude/longitude pair are outside of the current chunk of undulation
    // values read into memory.
    bool is_outside_cache_bounds(double lat, double lon) const;

    // Moves file pointer to the start of a new chunk of data and reads the chunk into memory.
    void init_coverage(double lat, double lon) const;

    // Set the bounds for a new square chunk to be pulled into memory. This chunk will be centered
    // around (lat, lon). Will set internal variables lat_coverage_min, lat_coverage_max,
    // lon_coverage_min, and lat_coverage_max to reflect the range of latitudes and longitudes that
    // the Matrix available_undulations contains undulation values for.
    void set_coverage_bounds(double lat, double lon) const;

    // Read in a 'box' of undulation data into memory.
    void read_chunk() const;

    // Skip a section of the file equivalent to a given number of records. This is called after
    // resetting the file pointer to the start of the file to move the pointer to the start of the
    // first record within the coverage bounds.
    void skip_records(size_t num) const;

    // Reads one double value from the file and returns it.
    double read_value() const;

    // Reads the value from a line in the file with a single element (such as occurs at the end of
    // each record). Required that the file pointer at the end of block 9 (the last block) in a
    // given record. This value is the same as the first value in a record, so no need to store it.
    void ignore_value() const;

    // Read an entire record (1440 values) starting at the current location of the file pointer. A
    // record consists of 9 blocks containing 20 lines with 8 values per line (for a total of 160
    // values per block). Each block is preceded by an empty line. A record is followed by an empty
    // line and an additional single value on its own line.
    void read_record(Size min_lon_idx, Size max_lon_idx, size_t row_idx) const;

    // Get the record number of a given latitude.  Ex: latitude 90 is the first record in the file,
    // so return 0; latitude 80 is the 40th record, so return 40.
    size_t get_lat_idx(double lat) const;

    // Get the element in the available_undulations Matrix corresponding to the given lat and lon.
    // Assumes lat and lon will be multiples of 0.25 and that the corresponding point will be found
    // in memory.
    double get_value(double lat, double lon) const;

    mutable std::shared_ptr<std::ifstream> infile;

    /* 'WW15MGH.GRD' file attributes */
    const double min_lat, max_lat, min_lon, max_lon, lat_step, lon_step;

    /* Chunk attributes */
    // chunk_size is the length of a "side" of a chunk square (i.e. the number of latitude points or
    // the number of longitude points to cover). The latitude or longitude range that a
    // available_undulations spans in chunk_size/4. Ex: A chunk_size of 40 includes 41x41
    // latitude/longitude points and spans 10 degrees of latitude and 10 degrees of longitude.
    Size chunk_size = 40;  // total number of bytes = (chunk_size + 1)^2 * (# of bytes/value)

    // latitude coverage range is [lat_coverage_min, lat_coverage_max]
    mutable double lat_coverage_min = -100;
    mutable double lat_coverage_max = -100;
    // longitude coverage range is [lon_coverage_min, lon_coverage_max]
    mutable double lon_coverage_min = -1;
    mutable double lon_coverage_max = -1;

    // The chunk of read data stored in memory
    // NOTE: using xt::xtensor instead of navtk::Matrix as the latter will heap-allocate a
    // xt::pytensor object when using the python module. This is undesirable as the
    // GeoidUndulationSource singleton won't be cleaned up until after the python interpreter has
    // completely shut down, which will trigger a segfault.
    mutable xt::xtensor<navtk::Scalar, 2> available_undulations;
};
}  // namespace geospatial
}  // namespace navtk