Branch data Line data Source code
1 : : #pragma once
2 : :
3 : : /*
4 : : MIT License
5 : :
6 : : Copyright (c) 2014-2024 Stephane Cuillerdier (aka aiekick)
7 : :
8 : : Permission is hereby granted, free of charge, to any person obtaining a copy
9 : : of this software and associated documentation files (the "Software"), to deal
10 : : in the Software without restriction, including without limitation the rights
11 : : to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 : : copies of the Software, and to permit persons to whom the Software is
13 : : furnished to do so, subject to the following conditions:
14 : :
15 : : The above copyright notice and this permission notice shall be included in all
16 : : copies or substantial portions of the Software.
17 : :
18 : : THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 : : IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 : : FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 : : AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 : : LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 : : OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24 : : SOFTWARE.
25 : : */
26 : :
27 : : // ezTile is part of the ezLibs project : https://github.com/aiekick/ezLibs.git
28 : :
29 : : #include <cmath>
30 : : #include <string>
31 : : #include <vector>
32 : : #include <cstdint>
33 : : #include <cassert>
34 : : #include <type_traits>
35 : :
36 : : #include "ezGeo.hpp"
37 : : #include "../ezMath.hpp"
38 : :
39 : : namespace ez {
40 : :
41 : : // --------------------------------------------------------------------------------------
42 : : // Convention used here (critical for consistency):
43 : : // - vec2<geo::dms> stores: x = latitude, y = longitude
44 : : // - Datas container is addressed as m_datas[latIndex][lonIndex]
45 : : // meaning "rows = latitudes, cols = longitudes"
46 : : // - Grid spacing is 1.0 degree in both latitude and longitude.
47 : : // --------------------------------------------------------------------------------------
48 : :
49 : : template <>
50 : : struct vec2<geo::dms> {
51 : : geo::dms x{}; // latitude
52 : : geo::dms y{}; // longitude
53 : :
54 : 10 : vec2() = default;
55 : :
56 : : // Constructor: (lat, lon) in decimal degrees
57 : 75 : explicit vec2(double vLat, double vLon) : x(geo::dms(vLat, true)), y(geo::dms(vLon, false)) {}
58 : :
59 : : // Convert to angles in decimal degrees (x = latDeg, y = lonDeg)
60 : 0 : dvec2 toAngle() const { return dvec2(x.toAngle(), y.toAngle()); }
61 : : };
62 : :
63 : : typedef vec2<geo::dms> dmsCoord;
64 : :
65 : : namespace geo {
66 : :
67 : : /*
68 : : tile is a container of geo datas built from a coord and a matrix buffer.
69 : : Coordinates are (lat, lon) mapped to (x, y) as stated above:
70 : : - x = latitude, y = longitude
71 : : - m_datas is indexed as m_datas[lat][lon]
72 : : */
73 : : template <typename TDATAS>
74 : : class tile {
75 : : public:
76 : : typedef std::vector<std::vector<TDATAS>> DatasContainer;
77 : :
78 : : private:
79 : : DatasContainer m_datas;
80 : : dmsCoord m_min; // min latitude (x) and min longitude (y) in DMS
81 : : dmsCoord m_max; // max latitude (x) and max longitude (y) in DMS
82 : : uint32_t m_nLats{}; // number of rows (latitude samples)
83 : : uint32_t m_nLons{}; // number of cols (longitude samples)
84 : : ez::range<TDATAS> m_range;
85 : : bool m_valid{false};
86 : :
87 : : public:
88 : : tile() = default;
89 : :
90 : : // vMin is the (lat, lon) of the minimal corner (in degrees)
91 : : // m_datas must be a rectangular matrix: rows = lats, cols = lons
92 : 10 : tile(const DatasContainer& vDatas, const dmsCoord& vMin) : m_datas(vDatas), m_min(vMin) { m_valid = check(); }
93 : :
94 : 10 : bool isValid() const { return m_valid; }
95 : :
96 : : // Counts
97 : 1 : uint32_t getNLats() const { return m_nLats; }
98 : 1 : uint32_t getNLons() const { return m_nLons; }
99 : :
100 : : // Min/Max DMS coords
101 : 1 : const dmsCoord& getMinDms() const { return m_min; }
102 : 1 : const dmsCoord& getMaxDms() const { return m_max; }
103 : :
104 : : // range
105 : 3 : ez::range<TDATAS> getRange() const { return m_range; }
106 : :
107 : : DatasContainer& getDatasRef() { return m_datas; }
108 : : const DatasContainer& getDatas() const { return m_datas; }
109 : :
110 : : // -----------------------------------------------------------------------------
111 : : // Discrete access by indices
112 : : // vCoord = {lonIndex, latIndex}
113 : : // Returns the raw TDATAS sample (no int32_terpolation).
114 : : // -----------------------------------------------------------------------------
115 : 88 : bool getValue(const uvec2& vCoord, TDATAS& voValue) const {
116 [ - + ]: 88 : if (!m_valid) {
117 : 0 : return false;
118 : 0 : }
119 : :
120 : 88 : const uint32_t lon = vCoord.x;
121 : 88 : const uint32_t lat = vCoord.y;
122 : :
123 [ + + ][ + + ]: 88 : if (lat >= m_nLats || lon >= m_nLons) {
124 : 7 : return false;
125 : 7 : }
126 : :
127 : 81 : voValue = m_datas.at(static_cast<size_t>(lat)).at(static_cast<size_t>(lon));
128 : 81 : return true;
129 : 88 : }
130 : :
131 : :
132 : : // -----------------------------------------------------------------------------
133 : : // Continuous access by geo coordinate (lat, lon in DMS)
134 : : // Bilinear int32_terpolation with optional wrapping on lat/lon.
135 : : // Always outputs a double in voValue (no int32_teger truncation).
136 : : // -----------------------------------------------------------------------------
137 : 30 : bool getValue(const dmsCoord& vCoord, double& voValue, bool vWrapLat = false, bool vWrapLon = false) const {
138 [ - + ]: 30 : if (!m_valid) {
139 : 0 : return false;
140 : 0 : }
141 : :
142 : : // Convert DMS to decimal degrees
143 : 30 : const double latMin = m_min.x.toAngle();
144 : 30 : const double lonMin = m_min.y.toAngle();
145 : 30 : const double latMax = m_max.x.toAngle();
146 : 30 : const double lonMax = m_max.y.toAngle();
147 : :
148 : 30 : double lat = vCoord.x.toAngle(); // latitude
149 : 30 : double lon = vCoord.y.toAngle(); // longitude
150 : :
151 : 30 : const double latSpan = static_cast<double>(m_nLats);
152 : 30 : const double lonSpan = static_cast<double>(m_nLons);
153 : :
154 : : // Optional wrapping int32_to tile range [min, min+span)
155 [ + + ]: 30 : if (vWrapLat) {
156 [ - + ]: 6 : if (latSpan <= 0.0) {
157 : 0 : return false;
158 : 0 : }
159 : 6 : double d = std::fmod(lat - latMin, latSpan);
160 [ - + ]: 6 : if (d < 0.0) {
161 : 0 : d += latSpan;
162 : 0 : }
163 : 6 : lat = latMin + d;
164 : 6 : }
165 [ + + ]: 30 : if (vWrapLon) {
166 [ - + ]: 6 : if (lonSpan <= 0.0) {
167 : 0 : return false;
168 : 0 : }
169 : 6 : double d = std::fmod(lon - lonMin, lonSpan);
170 [ - + ]: 6 : if (d < 0.0) {
171 : 0 : d += lonSpan;
172 : 0 : }
173 : 6 : lon = lonMin + d;
174 : 6 : }
175 : :
176 : : // Strict bounds if no wrap
177 [ + + ][ + + ]: 30 : if (!vWrapLat && (lat < latMin || lat > latMax)) {
[ + + ]
178 : 3 : return false;
179 : 3 : }
180 [ + + ][ + + ]: 27 : if (!vWrapLon && (lon < lonMin || lon > lonMax)) {
[ + + ]
181 : 3 : return false;
182 : 3 : }
183 : :
184 : : // Fractional indices in the grid (1 degree per step)
185 : 24 : double fLat = (lat - latMin);
186 : 24 : double fLon = (lon - lonMin);
187 : :
188 : : // Clamp int32_to [0, n-1] to keep indices valid
189 [ - + ]: 24 : if (fLat < 0.0) {
190 : 0 : fLat = 0.0;
191 : 0 : }
192 [ - + ]: 24 : if (fLon < 0.0) {
193 : 0 : fLon = 0.0;
194 : 0 : }
195 [ - + ]: 24 : if (fLat > static_cast<double>(m_nLats - 1)) {
196 : 0 : fLat = static_cast<double>(m_nLats - 1);
197 : 0 : }
198 [ + + ]: 24 : if (fLon > static_cast<double>(m_nLons - 1)) {
199 : 6 : fLon = static_cast<double>(m_nLons - 1);
200 : 6 : }
201 : :
202 : 24 : const int32_t iLat = static_cast<int32_t>(std::floor(fLat));
203 : 24 : const int32_t iLon = static_cast<int32_t>(std::floor(fLon));
204 : 24 : const double tLat = fLat - static_cast<double>(iLat);
205 : 24 : const double tLon = fLon - static_cast<double>(iLon);
206 : :
207 : : // Helper to wrap/clamp indices according to flags and sizes
208 : 96 : auto wrapOrClamp = [](int32_t idx, int32_t size, bool doWrap) -> uint32_t {
209 [ + + ]: 96 : if (doWrap) {
210 [ - + ]: 24 : if (size <= 0) {
211 : 0 : return 0u;
212 : 0 : }
213 : 24 : int32_t m = idx % size;
214 [ - + ]: 24 : if (m < 0) {
215 : 0 : m += size;
216 : 0 : }
217 : 24 : return static_cast<uint32_t>(m);
218 : 72 : } else {
219 [ - + ]: 72 : if (idx < 0) {
220 : 0 : idx = 0;
221 : 0 : }
222 [ + + ]: 72 : if (idx >= size) {
223 : 12 : idx = size - 1;
224 : 12 : }
225 : 72 : return static_cast<uint32_t>(idx);
226 : 72 : }
227 : 96 : };
228 : :
229 : : // Compute neighbor indices (with wrap if requested)
230 : 24 : const uint32_t iLat0 = wrapOrClamp(iLat, static_cast<int32_t>(m_nLats), vWrapLat);
231 : 24 : const uint32_t iLat1 = wrapOrClamp(iLat + 1, static_cast<int32_t>(m_nLats), vWrapLat);
232 : 24 : const uint32_t iLon0 = wrapOrClamp(iLon, static_cast<int32_t>(m_nLons), vWrapLon);
233 : 24 : const uint32_t iLon1 = wrapOrClamp(iLon + 1, static_cast<int32_t>(m_nLons), vWrapLon);
234 : :
235 : : // If no neighbor and no wrap requested, fallback to nearest neighbor
236 [ + + ][ + + ]: 24 : if ((!vWrapLat && (iLat + 1 >= static_cast<int32_t>(m_nLats))) || (!vWrapLon && (iLon + 1 >= static_cast<int32_t>(m_nLons)))) {
[ + + ][ + + ]
237 : 9 : TDATAS nearest{};
238 [ - + ]: 9 : if (!getValue(uvec2(iLon0, iLat0), nearest)) {
239 : 0 : return false;
240 : 0 : }
241 : 9 : voValue = static_cast<double>(nearest);
242 : 9 : return true;
243 : 9 : }
244 : :
245 : : // Fetch the 4 corners (order: v00=(lat0,lon0), v10=(lat1,lon0), v01=(lat0,lon1), v11=(lat1,lon1))
246 : 15 : TDATAS v00{};
247 : 15 : TDATAS v10{};
248 : 15 : TDATAS v01{};
249 : 15 : TDATAS v11{};
250 : 15 : bool ok = true;
251 : 15 : ok &= getValue(uvec2(iLon0, iLat0), v00);
252 : 15 : ok &= getValue(uvec2(iLon0, iLat1), v10);
253 : 15 : ok &= getValue(uvec2(iLon1, iLat0), v01);
254 : 15 : ok &= getValue(uvec2(iLon1, iLat1), v11);
255 [ - + ]: 15 : if (!ok) {
256 : 0 : return false;
257 : 0 : }
258 : :
259 : : // Bilinear int32_terpolation (lat vertically, lon horizontally)
260 : 15 : const double d00 = static_cast<double>(v00);
261 : 15 : const double d10 = static_cast<double>(v10);
262 : 15 : const double d01 = static_cast<double>(v01);
263 : 15 : const double d11 = static_cast<double>(v11);
264 : :
265 : 15 : const double a = d00 * (1.0 - tLat) + d10 * tLat; // along latitude
266 : 15 : const double b = d01 * (1.0 - tLat) + d11 * tLat; // along latitude
267 : 15 : const double d = a * (1.0 - tLon) + b * tLon; // along longitude
268 : :
269 : 15 : voValue = d;
270 : 15 : return true;
271 : 15 : }
272 : :
273 : : // Validate the matrix, set sizes and deduce m_max from m_min and dimensions.
274 : 10 : bool check() {
275 [ - + ]: 10 : if (m_datas.empty()) {
276 : 0 : #ifdef EZ_TOOLS_LOG
277 : 0 : LogVarError(u8R"(tile: data matrix is empty)");
278 : 0 : #endif
279 : 0 : return false;
280 : 0 : }
281 : :
282 : 10 : m_range = {};
283 : :
284 : : // rows = latitudes
285 : 10 : m_nLats = static_cast<uint32_t>(m_datas.size());
286 : :
287 : : // ensure rectangular matrix and compute m_nLons
288 : 10 : m_nLons = 0u;
289 [ + + ]: 20 : for (const auto& row : m_datas) {
290 [ - + ]: 20 : if (row.empty()) {
291 : 0 : #ifdef EZ_TOOLS_LOG
292 : 0 : LogVarError(u8R"(tile: a row is empty)");
293 : 0 : #endif
294 : 0 : return false;
295 : 0 : }
296 [ + + ]: 20 : if (m_nLons == 0u) {
297 : 10 : m_nLons = static_cast<uint32_t>(row.size());
298 [ - + ]: 10 : } else if (m_nLons != static_cast<uint32_t>(row.size())) {
299 : 0 : #ifdef EZ_TOOLS_LOG
300 : 0 : LogVarError(u8R"(tile: non-rectangular matrix (row sizes differ))");
301 : 0 : #endif
302 : 0 : return false;
303 : 0 : }
304 [ + + ]: 40 : for (const auto& v : row) {
305 : 40 : m_range.combine(v);
306 : 40 : }
307 : 20 : }
308 : :
309 : : // Compute max bounds from min + extents (1 degree per step)
310 : 10 : m_max = m_min;
311 : : // x = latitude, y = longitude
312 : 10 : m_max.x.offsetDeg(static_cast<int32_t>(m_nLats));
313 : 10 : m_max.y.offsetDeg(static_cast<int32_t>(m_nLons));
314 : :
315 : 10 : return true;
316 : 10 : }
317 : : };
318 : :
319 : : } // namespace geo
320 : : } // namespace ez
|
