1 #pragma once
2 
3 #include <mbgl/util/constants.hpp>
4 #include <mbgl/util/geo.hpp>
5 #include <mbgl/util/geometry.hpp>
6 #include <mbgl/math/clamp.hpp>
7 
8 namespace mbgl {
9 
10 class ProjectedMeters {
11 private:
12     double _northing; // Distance measured northwards.
13     double _easting;  // Distance measured eastwards.
14 
15 public:
ProjectedMeters(double n_=0,double e_=0)16     ProjectedMeters(double n_ = 0, double e_ = 0)
17         : _northing(n_), _easting(e_) {
18         if (std::isnan(_northing)) {
19             throw std::domain_error("northing must not be NaN");
20         }
21         if (std::isnan(_easting)) {
22             throw std::domain_error("easting must not be NaN");
23         }
24     }
25 
northing() const26     double northing() const { return _northing; }
easting() const27     double easting() const { return _easting; }
28 
operator ==(const ProjectedMeters & a,const ProjectedMeters & b)29     friend bool operator==(const ProjectedMeters& a, const ProjectedMeters& b) {
30         return a._northing == b._northing && a._easting == b._easting;
31     }
32 
operator !=(const ProjectedMeters & a,const ProjectedMeters & b)33     friend bool operator!=(const ProjectedMeters& a, const ProjectedMeters& b) {
34         return !(a == b);
35     }
36 };
37 
38 // Spherical Mercator projection
39 // http://docs.openlayers.org/library/spherical_mercator.html
40 class Projection {
41 public:
42     // Map pixel width at given scale.
worldSize(double scale)43     static double worldSize(double scale) {
44         return scale * util::tileSize;
45     }
46 
getMetersPerPixelAtLatitude(double lat,double zoom)47     static double getMetersPerPixelAtLatitude(double lat, double zoom) {
48         const double constrainedZoom = util::clamp(zoom, util::MIN_ZOOM, util::MAX_ZOOM);
49         const double constrainedScale = std::pow(2.0, constrainedZoom);
50         const double constrainedLatitude = util::clamp(lat, -util::LATITUDE_MAX, util::LATITUDE_MAX);
51         return std::cos(constrainedLatitude * util::DEG2RAD) * util::M2PI * util::EARTH_RADIUS_M / worldSize(constrainedScale);
52     }
53 
projectedMetersForLatLng(const LatLng & latLng)54     static ProjectedMeters projectedMetersForLatLng(const LatLng& latLng) {
55         const double constrainedLatitude = util::clamp(latLng.latitude(), -util::LATITUDE_MAX, util::LATITUDE_MAX);
56         const double constrainedLongitude = util::clamp(latLng.longitude(), -util::LONGITUDE_MAX, util::LONGITUDE_MAX);
57 
58         const double m = 1 - 1e-15;
59         const double f = util::clamp(std::sin(util::DEG2RAD * constrainedLatitude), -m, m);
60 
61         const double easting  = util::EARTH_RADIUS_M * constrainedLongitude * util::DEG2RAD;
62         const double northing = 0.5 * util::EARTH_RADIUS_M * std::log((1 + f) / (1 - f));
63 
64         return ProjectedMeters(northing, easting);
65     }
66 
latLngForProjectedMeters(const ProjectedMeters & projectedMeters)67     static LatLng latLngForProjectedMeters(const ProjectedMeters& projectedMeters) {
68         double latitude = (2 * std::atan(std::exp(projectedMeters.northing() / util::EARTH_RADIUS_M)) - (M_PI / 2.0)) * util::RAD2DEG;
69         double longitude = projectedMeters.easting() * util::RAD2DEG / util::EARTH_RADIUS_M;
70 
71         latitude = util::clamp(latitude, -util::LATITUDE_MAX, util::LATITUDE_MAX);
72         longitude = util::clamp(longitude, -util::LONGITUDE_MAX, util::LONGITUDE_MAX);
73 
74         return LatLng(latitude, longitude);
75     }
76 
project(const LatLng & latLng,double scale)77     static Point<double> project(const LatLng& latLng, double scale) {
78         return project_(latLng, worldSize(scale));
79     }
80 
81     //Returns point on tile
project(const LatLng & latLng,int32_t zoom)82     static Point<double> project(const LatLng& latLng, int32_t zoom) {
83         return project_(latLng, 1 << zoom);
84     }
85 
unproject(const Point<double> & p,double scale,LatLng::WrapMode wrapMode=LatLng::Unwrapped)86     static LatLng unproject(const Point<double>& p, double scale, LatLng::WrapMode wrapMode = LatLng::Unwrapped) {
87         auto p2 = p * util::DEGREES_MAX / worldSize(scale);
88         return LatLng {
89             util::DEGREES_MAX / M_PI * std::atan(std::exp((util::LONGITUDE_MAX - p2.y) * util::DEG2RAD)) - 90.0,
90             p2.x - util::LONGITUDE_MAX,
91             wrapMode
92         };
93     }
94 
95 private:
project_(const LatLng & latLng,double worldSize)96     static Point<double> project_(const LatLng& latLng, double worldSize) {
97         const double latitude = util::clamp(latLng.latitude(), -util::LATITUDE_MAX, util::LATITUDE_MAX);
98         return Point<double> {
99             util::LONGITUDE_MAX + latLng.longitude(),
100             util::LONGITUDE_MAX - util::RAD2DEG * std::log(std::tan(M_PI / 4 + latitude * M_PI / util::DEGREES_MAX))
101         } * worldSize / util::DEGREES_MAX;
102     }
103 };
104 
105 } // namespace mbgl
106