]> git.mxchange.org Git - simgear.git/blob - simgear/bucket/newbucket.cxx
Ralf Gerlich: fix bucket numbering at extreme latitudes
[simgear.git] / simgear / bucket / newbucket.cxx
1 /**************************************************************************
2  * newbucket.hxx -- new bucket routines for better world modeling
3  *
4  * Written by Curtis L. Olson, started February 1999.
5  *
6  * Copyright (C) 1999  Curtis L. Olson - http://www.flightgear.org/~curt/
7  *
8  * This library is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU Library General Public
10  * License as published by the Free Software Foundation; either
11  * version 2 of the License, or (at your option) any later version.
12  *
13  * This library is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16  * Library General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this program; if not, write to the Free Software
20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
21  *
22  * $Id$
23  **************************************************************************/
24
25
26 #ifdef HAVE_CONFIG_H
27 #  include <simgear_config.h>
28 #endif
29
30 #include <math.h>
31
32 #include <simgear/misc/sg_path.hxx>
33
34 #include "newbucket.hxx"
35
36
37 // default constructor
38 SGBucket::SGBucket() {
39 }
40
41
42 // constructor for specified location
43 SGBucket::SGBucket(const double dlon, const double dlat) {
44     set_bucket(dlon, dlat);
45 }
46
47 SGBucket::SGBucket(const SGGeod& geod) {
48     set_bucket(geod);
49 }
50
51 // create an impossible bucket if false
52 SGBucket::SGBucket(const bool is_good) {
53     set_bucket(0.0, 0.0);
54     if ( !is_good ) {
55         lon = -1000;
56     }
57 }
58
59
60 // Parse a unique scenery tile index and find the lon, lat, x, and y
61 SGBucket::SGBucket(const long int bindex) {
62     long int index = bindex;
63         
64     lon = index >> 14;
65     index -= lon << 14;
66     lon -= 180;
67
68     lat = index >> 6;
69     index -= lat << 6;
70     lat -= 90;
71
72     y = index >> 3;
73     index -= y << 3;
74
75     x = index;
76 }
77
78
79 // Set the bucket params for the specified lat and lon
80 void SGBucket::set_bucket( double *lonlat ) {
81     set_bucket( lonlat[0], lonlat[1] );
82 }       
83
84
85 // Set the bucket params for the specified lat and lon
86 void SGBucket::set_bucket( double dlon, double dlat ) {
87     //
88     // latitude first
89     //
90     double span = sg_bucket_span( dlat );
91     double diff = dlon - (double)(int)dlon;
92
93     // cout << "diff = " << diff << "  span = " << span << endl;
94
95     /* Calculate the greatest integral longitude less than
96      * or equal to the given longitude (floor(dlon)),
97      * but attribute coordinates near the east border
98      * to the next tile.
99      */
100     if ( (dlon >= 0) || (fabs(diff) < SG_EPSILON) ) {
101         lon = (int)dlon;
102     } else {
103         lon = (int)dlon - 1;
104     }
105
106     // find subdivision or super lon if needed
107     if ( span < SG_EPSILON ) {
108         /* sg_bucket_span() never returns 0.0
109          * or anything near it, so this really
110          * should not occur at any time.
111          */
112         // polar cap
113         lon = 0;
114         x = 0;
115     } else if ( span <= 1.0 ) {
116         /* We have more than one tile per degree of
117          * longitude, so we need an x offset.
118          */
119         x = (int)((dlon - lon) / span);
120     } else {
121         /* We have one or more degrees per tile,
122          * so we need to find the base longitude
123          * of that tile.
124          *
125          * First we calculate the integral base longitude
126          * (e.g. -85.5 => -86) and then find the greatest
127          * multiple of span that is less than or equal to
128          * that longitude.
129          *
130          * That way, the Greenwich Meridian is always
131          * a tile border.
132          *
133          * This gets us into trouble with the polar caps,
134          * which have width 360 and thus either span
135          * the range from 0 to 360 or from -360 to 0
136          * degrees, depending on whether lon is positive
137          * or negative!
138          *
139          * We also get into trouble with the 8 degree tiles
140          * north of 88N and south of 88S, because the west-
141          * and east-most tiles in that range will cover 184W
142          * to 176W and 176E to 184E respectively, with their
143          * center at 180E/W!
144          */
145         lon=(int)floor(floor((lon+SG_EPSILON)/span)*span);
146         /* Correct the polar cap issue */
147         if ( lon < -180 ) {
148             lon = -180;
149         }
150         x = 0;
151     }
152
153     //
154     // then latitude
155     //
156     diff = dlat - (double)(int)dlat;
157
158     /* Again, a modified floor() function (see longitude) */
159     if ( (dlat >= 0) || (fabs(diff) < SG_EPSILON) ) {
160         lat = (int)dlat;
161     } else {
162         lat = (int)dlat - 1;
163     }
164     /* Latitude base and offset are easier, as
165      * tiles always are 1/8 degree of latitude wide.
166      */
167     y = (int)((dlat - lat) * 8);
168 }
169
170
171 void SGBucket::set_bucket(const SGGeod& geod)
172 {
173     set_bucket(geod.getLongitudeDeg(), geod.getLatitudeDeg());
174 }
175
176 // Build the path name for this bucket
177 std::string SGBucket::gen_base_path() const {
178     // long int index;
179     int top_lon, top_lat, main_lon, main_lat;
180     char hem, pole;
181     char raw_path[256];
182
183     top_lon = lon / 10;
184     main_lon = lon;
185     if ( (lon < 0) && (top_lon * 10 != lon) ) {
186         top_lon -= 1;
187     }
188     top_lon *= 10;
189     if ( top_lon >= 0 ) {
190         hem = 'e';
191     } else {
192         hem = 'w';
193         top_lon *= -1;
194     }
195     if ( main_lon < 0 ) {
196         main_lon *= -1;
197     }
198     
199     top_lat = lat / 10;
200     main_lat = lat;
201     if ( (lat < 0) && (top_lat * 10 != lat) ) {
202         top_lat -= 1;
203     }
204     top_lat *= 10;
205     if ( top_lat >= 0 ) {
206         pole = 'n';
207     } else {
208         pole = 's';
209         top_lat *= -1;
210     }
211     if ( main_lat < 0 ) {
212         main_lat *= -1;
213     }
214
215     sprintf(raw_path, "%c%03d%c%02d/%c%03d%c%02d", 
216             hem, top_lon, pole, top_lat, 
217             hem, main_lon, pole, main_lat);
218
219     SGPath path( raw_path );
220
221     return path.str();
222 }
223
224
225 // return width of the tile in degrees
226 double SGBucket::get_width() const {
227     if (lon==-180 && (lat==-89 || lat==88) ) {
228         /* Normally the tile at 180W in 88N and 89S
229          * would cover 184W to 176W and the next
230          * on the east side starts at 176W.
231          * To correct, make this a special tile
232          * from 180W to 176W with 4 degrees width
233          * instead of the normal 8 degrees at
234          * that latitude.
235          */
236          return 4.0;
237     }
238     return sg_bucket_span( get_center_lat() );
239 }
240
241
242 // return height of the tile in degrees
243 double SGBucket::get_height() const {
244     return SG_BUCKET_SPAN;
245 }
246
247
248 // return width of the tile in meters
249 double SGBucket::get_width_m() const {
250     double clat = (int)get_center_lat();
251     if ( clat > 0 ) {
252         clat = (int)clat + 0.5;
253     } else {
254         clat = (int)clat - 0.5;
255     }
256     double clat_rad = clat * SGD_DEGREES_TO_RADIANS;
257     double cos_lat = cos( clat_rad );
258     double local_radius = cos_lat * SG_EQUATORIAL_RADIUS_M;
259     double local_perimeter = local_radius * SGD_2PI;
260     double degree_width = local_perimeter / 360.0;
261
262     return get_width() * degree_width;
263 }
264
265
266 // return height of the tile in meters
267 double SGBucket::get_height_m() const {
268     double perimeter = SG_EQUATORIAL_RADIUS_M * SGD_2PI;
269     double degree_height = perimeter / 360.0;
270
271     return SG_BUCKET_SPAN * degree_height;
272 }
273
274
275 // find the bucket which is offset by the specified tile units in the
276 // X & Y direction.  We need the current lon and lat to resolve
277 // ambiguities when going from a wider tile to a narrower one above or
278 // below.  This assumes that we are feeding in
279 SGBucket sgBucketOffset( double dlon, double dlat, int dx, int dy ) {
280     SGBucket result( dlon, dlat );
281     double clat = result.get_center_lat() + dy * SG_BUCKET_SPAN;
282
283     // walk dy units in the lat direction
284     result.set_bucket( dlon, clat );
285
286     // find the lon span for the new latitude
287     double span = sg_bucket_span( clat );
288
289     // walk dx units in the lon direction
290     double tmp = dlon + dx * span;
291     while ( tmp < -180.0 ) {
292         tmp += 360.0;
293     }
294     while ( tmp >= 180.0 ) {
295         tmp -= 360.0;
296     }
297     result.set_bucket( tmp, clat );
298
299     return result;
300 }
301
302
303 // calculate the offset between two buckets
304 void sgBucketDiff( const SGBucket& b1, const SGBucket& b2, int *dx, int *dy ) {
305
306     // Latitude difference
307     double c1_lat = b1.get_center_lat();
308     double c2_lat = b2.get_center_lat();
309     double diff_lat = c2_lat - c1_lat;
310
311 #ifdef HAVE_RINT
312     *dy = (int)rint( diff_lat / SG_BUCKET_SPAN );
313 #else
314     if ( diff_lat > 0 ) {
315         *dy = (int)( diff_lat / SG_BUCKET_SPAN + 0.5 );
316     } else {
317         *dy = (int)( diff_lat / SG_BUCKET_SPAN - 0.5 );
318     }
319 #endif
320
321     // longitude difference
322     double diff_lon=0.0;
323     double span=0.0;
324
325     SGBucket tmp_bucket;
326     // To handle crossing the bucket size boundary
327     //  we need to account for different size buckets.
328
329     if ( sg_bucket_span(c1_lat) <= sg_bucket_span(c2_lat) )
330     {
331         span = sg_bucket_span(c1_lat);
332     } else {
333         span = sg_bucket_span(c2_lat);
334     }
335
336     diff_lon = b2.get_center_lon() - b1.get_center_lon();
337
338     if (diff_lon <0.0)
339     {
340        diff_lon -= b1.get_width()*0.5 + b2.get_width()*0.5 - span;
341     } 
342     else
343     {
344        diff_lon += b1.get_width()*0.5 + b2.get_width()*0.5 - span;
345     }
346
347
348 #ifdef HAVE_RINT
349     *dx = (int)rint( diff_lon / span );
350 #else
351     if ( diff_lon > 0 ) {
352         *dx = (int)( diff_lon / span + 0.5 );
353     } else {
354         *dx = (int)( diff_lon / span - 0.5 );
355     }
356 #endif
357 }
358
359