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: |
16 | ProjectedMeters(double n_ = 0, double e_ = 0) |
17 | : _northing(n_), _easting(e_) { |
18 | if (std::isnan(x: _northing)) { |
19 | throw std::domain_error("northing must not be NaN" ); |
20 | } |
21 | if (std::isnan(x: _easting)) { |
22 | throw std::domain_error("easting must not be NaN" ); |
23 | } |
24 | } |
25 | |
26 | double northing() const { return _northing; } |
27 | double easting() const { return _easting; } |
28 | |
29 | friend bool operator==(const ProjectedMeters& a, const ProjectedMeters& b) { |
30 | return a._northing == b._northing && a._easting == b._easting; |
31 | } |
32 | |
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. |
43 | static double worldSize(double scale) { |
44 | return scale * util::tileSize; |
45 | } |
46 | |
47 | static double getMetersPerPixelAtLatitude(double lat, double zoom) { |
48 | const double constrainedZoom = util::clamp(value: zoom, min_: util::MIN_ZOOM, max_: util::MAX_ZOOM); |
49 | const double constrainedScale = std::pow(x: 2.0, y: constrainedZoom); |
50 | const double constrainedLatitude = util::clamp(value: lat, min_: -util::LATITUDE_MAX, max_: util::LATITUDE_MAX); |
51 | return std::cos(x: constrainedLatitude * util::DEG2RAD) * util::M2PI * util::EARTH_RADIUS_M / worldSize(scale: constrainedScale); |
52 | } |
53 | |
54 | static ProjectedMeters projectedMetersForLatLng(const LatLng& latLng) { |
55 | const double constrainedLatitude = util::clamp(value: latLng.latitude(), min_: -util::LATITUDE_MAX, max_: util::LATITUDE_MAX); |
56 | const double constrainedLongitude = util::clamp(value: latLng.longitude(), min_: -util::LONGITUDE_MAX, max_: util::LONGITUDE_MAX); |
57 | |
58 | const double m = 1 - 1e-15; |
59 | const double f = util::clamp(value: std::sin(x: util::DEG2RAD * constrainedLatitude), min_: -m, max_: 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(x: (1 + f) / (1 - f)); |
63 | |
64 | return ProjectedMeters(northing, easting); |
65 | } |
66 | |
67 | static LatLng latLngForProjectedMeters(const ProjectedMeters& projectedMeters) { |
68 | double latitude = (2 * std::atan(x: std::exp(x: 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(value: latitude, min_: -util::LATITUDE_MAX, max_: util::LATITUDE_MAX); |
72 | longitude = util::clamp(value: longitude, min_: -util::LONGITUDE_MAX, max_: util::LONGITUDE_MAX); |
73 | |
74 | return LatLng(latitude, longitude); |
75 | } |
76 | |
77 | static Point<double> project(const LatLng& latLng, double scale) { |
78 | return project_(latLng, worldSize: worldSize(scale)); |
79 | } |
80 | |
81 | //Returns point on tile |
82 | static Point<double> project(const LatLng& latLng, int32_t zoom) { |
83 | return project_(latLng, worldSize: 1 << zoom); |
84 | } |
85 | |
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(x: std::exp(x: (util::LONGITUDE_MAX - p2.y) * util::DEG2RAD)) - 90.0, |
90 | p2.x - util::LONGITUDE_MAX, |
91 | wrapMode |
92 | }; |
93 | } |
94 | |
95 | private: |
96 | static Point<double> project_(const LatLng& latLng, double worldSize) { |
97 | const double latitude = util::clamp(value: latLng.latitude(), min_: -util::LATITUDE_MAX, max_: util::LATITUDE_MAX); |
98 | return Point<double> { |
99 | util::LONGITUDE_MAX + latLng.longitude(), |
100 | util::LONGITUDE_MAX - util::RAD2DEG * std::log(x: std::tan(M_PI / 4 + latitude * M_PI / util::DEGREES_MAX)) |
101 | } * worldSize / util::DEGREES_MAX; |
102 | } |
103 | }; |
104 | |
105 | } // namespace mbgl |
106 | |