coordinates.hpp

// Copyright 2021 Arnaud Becheler    <abechele@umich.edu>
 
/*********************************************************************** * This program is free software; you can
 *redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free
 *Software Foundation; either version 2 of the License, or    * (at your option) any later version. *
 *                                                                      *
 ***************************************************************************/
 
#pragma once
 
#include <assert.h>
#include <cmath> // trigonometry
#include <ostream>
#include <tuple>
#include <numbers>
 
#include <mp-units/systems/si.h>
#include <mp-units/ostream.h>
 
namespace quetzal::geography
{
 
using std::numbers::pi;
 
class latlon;
class lonlat;
class rowcol;
class colrow;
 
template <typename T>
concept GridCoordinate = requires(T t) {
                             {
                                 t.col <=> t.row
                             };
                         };
 
template <GridCoordinate C1, GridCoordinate C2> auto operator==(const C1 &c1, const C2 &c2)
{
    return std::tie(c1.col, c1.row) == std::tie(c2.col, c2.row);
}
 
template <GridCoordinate C1, GridCoordinate C2> auto operator<=>(const C1 &c1, const C2 &c2)
{
    return std::tie(c1.col, c1.row) <=> std::tie(c2.col, c2.row);
}
 
template <GridCoordinate C> std::ostream &operator<<(std::ostream &stream, const C &c)
{
    stream << "(Column: " << c.col << ", Row: " << c.row << ")";
    return stream;
}
 
struct rowcol
{
    using value_type = int;
 
    value_type row;
 
    value_type col;
 
    constexpr explicit rowcol(value_type row, value_type column) noexcept : row(row), col(column)
    {
    }
};
static_assert(GridCoordinate<rowcol>);
 
struct colrow
{
    using value_type = int;
 
    value_type col;
 
    value_type row;
 
    constexpr explicit colrow(value_type column, value_type row) noexcept : col(column), row(row)
    {
    }
};
static_assert(GridCoordinate<colrow>);
 
template <typename T>
concept Coordinate = requires(T t) {
                         {
                             t.lon <=> t.lat
                         };
                     };
 
template <Coordinate C1, Coordinate C2> auto operator==(const C1 &c1, const C2 &c2)
{
    return std::tie(c1.lat, c1.lon) == std::tie(c2.lat, c2.lon);
}
 
template <Coordinate C1, Coordinate C2> auto operator<=>(const C1 &c1, const C2 &c2)
{
    return std::tie(c1.lat, c1.lon) <=> std::tie(c2.lat, c2.lon);
}
 
template <Coordinate C> std::ostream &operator<<(std::ostream &stream, const C &c)
{
    stream << "(Lat: " << c.lat << ", Lon: " << c.lon << ")";
    return stream;
}
 
namespace details
{
 
void check(double lat, double lon)
{
    if ((lat > 90.) || (lat < -90.) || (lon < -180.) || (lon > 180.))
        throw std::invalid_argument("invalid coordinates");
}
 
constexpr double deg2rad(double deg) noexcept
{
    return (deg * pi / 180.0);
}
 
constexpr double rad2deg(double rad) noexcept
{
    return (rad * 180.0 / pi);
}
 
constexpr double distanceEarth(double lat1d, double lon1d, double lat2d, double lon2d)
{
 
    const double earthRadiusKm = 6371.0;
 
    const double lat1r = deg2rad(lat1d);
    const double lon1r = deg2rad(lon1d);
    const double lat2r = deg2rad(lat2d);
    const double lon2r = deg2rad(lon2d);
 
    const double u = sin((lat2r - lat1r) / 2);
    const double v = sin((lon2r - lon1r) / 2);
 
    const double d = 2.0 * earthRadiusKm * asin(sqrt(u * u + cos(lat1r) * cos(lat2r) * v * v));
    assert(d >= 0);
    return d;
}
 
} // namespace details
 
struct latlon
{
    using decimal_degree = double;
 
    decimal_degree lat;
 
    decimal_degree lon;
 
    constexpr explicit latlon(decimal_degree latitude, decimal_degree longitude) noexcept
        : lat(latitude), lon(longitude)
    {
    }
 
    template <class T> constexpr 
    mp_units::quantity<mp_units::si::metre> great_circle_distance_to(const T &other) const noexcept
    {
        const double d = details::distanceEarth(this->lat, this->lon, other.lat, other.lon);
        assert(d >= 0.);
        return d * 1000.0 * mp_units::si::metre;
    }
};
 
static_assert(Coordinate<latlon>);
 
struct lonlat
{
    using decimal_degree = double;
 
    decimal_degree lon;
 
    decimal_degree lat;
 
    constexpr explicit lonlat(decimal_degree longitude, decimal_degree latitude) noexcept
        : lon(longitude), lat(latitude)
    {
    }
 
    template <class T> 
    constexpr 
    mp_units::quantity<mp_units::si::metre> great_circle_distance_to(const T &other) const noexcept
    {
        const double d = details::distanceEarth(this->lat, this->lon, other.lat, other.lon);
        assert(d >= 0.0);
        return d * 1000.0 * mp_units::si::metre;
    }
 
}; // end struct latlon
static_assert(Coordinate<lonlat>);
 
} // namespace quetzal::geography