Imagery.cc 20.5 KB
Newer Older
Lionel Heng's avatar
Lionel Heng committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
/*=====================================================================

QGroundControl Open Source Ground Control Station

(c) 2009, 2010 QGROUNDCONTROL PROJECT <http://www.qgroundcontrol.org>

This file is part of the QGROUNDCONTROL project

    QGROUNDCONTROL 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 3 of the License, or
    (at your option) any later version.

    QGROUNDCONTROL is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with QGROUNDCONTROL. If not, see <http://www.gnu.org/licenses/>.

======================================================================*/

/**
 * @file
 *   @brief Definition of the class Imagery.
 *
 *   @author Lionel Heng <hengli@student.ethz.ch>
 *
 */

#include "Imagery.h"

#include <cmath>
#include <iomanip>
#include <sstream>

const double WGS84_A = 6378137.0;
const double WGS84_ECCSQ = 0.00669437999013;

41
const int MAX_ZOOM_LEVEL = 20;
Lionel Heng's avatar
Lionel Heng committed
42
43
44
45
46
47
48

Imagery::Imagery()
    : textureCache(new TextureCache(1000))
{

}

49
50
51
52
53
54
Imagery::ImageryType
Imagery::getImageryType(void) const
{
    return currentImageryType;
}

Lionel Heng's avatar
Lionel Heng committed
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
void
Imagery::setImageryType(ImageryType type)
{
    currentImageryType = type;
}

void
Imagery::setOffset(double xOffset, double yOffset)
{
    this->xOffset = xOffset;
    this->yOffset = yOffset;
}

void
Imagery::prefetch2D(double windowWidth, double windowHeight,
                    double zoom, double xOrigin, double yOrigin,
                    const QString& utmZone)
{
73
    if (currentImageryType == BLANK_MAP) {
74
75
76
        return;
    }

77
    double tileResolution = 1.0;
Lionel Heng's avatar
Lionel Heng committed
78
    if (currentImageryType == GOOGLE_SATELLITE ||
79
            currentImageryType == GOOGLE_MAP) {
Lionel Heng's avatar
Lionel Heng committed
80
        tileResolution = 1.0;
81
        while (tileResolution * 3.0 / 2.0 < 1.0 / zoom) {
Lionel Heng's avatar
Lionel Heng committed
82
83
            tileResolution *= 2.0;
        }
84
        if (tileResolution > 512.0) {
Lionel Heng's avatar
Lionel Heng committed
85
86
            tileResolution = 512.0;
        }
87
    } else if (currentImageryType == SWISSTOPO_SATELLITE) {
Lionel Heng's avatar
Lionel Heng committed
88
89
90
        tileResolution = 0.25;
    }

91
92
    int minTileX, minTileY, maxTileX, maxTileY;
    int zoomLevel;
Lionel Heng's avatar
Lionel Heng committed
93
94

    tileBounds(tileResolution,
95
96
97
98
               xOrigin - windowWidth / 2.0 / zoom * 1.5,
               yOrigin - windowHeight / 2.0 / zoom * 1.5,
               xOrigin + windowWidth / 2.0 / zoom * 1.5,
               yOrigin + windowHeight / 2.0 / zoom * 1.5, utmZone,
Lionel Heng's avatar
Lionel Heng committed
99
100
               minTileX, minTileY, maxTileX, maxTileY, zoomLevel);

101
102
    for (int r = minTileY; r <= maxTileY; ++r) {
        for (int c = minTileX; c <= maxTileX; ++c) {
Lionel Heng's avatar
Lionel Heng committed
103
104
105
106
107
108
109
110
111
112
            QString url = getTileLocation(c, r, zoomLevel, tileResolution);

            TexturePtr t = textureCache->get(url);
        }
    }
}

void
Imagery::draw2D(double windowWidth, double windowHeight,
                double zoom, double xOrigin, double yOrigin,
113
                double xOffset, double yOffset, double zOffset,
Lionel Heng's avatar
Lionel Heng committed
114
115
                const QString& utmZone)
{
116
    if (getNumDrawables() > 0) {
117
118
119
        removeDrawables(0, getNumDrawables());
    }

120
    if (currentImageryType == BLANK_MAP) {
121
122
123
        return;
    }

124
    double tileResolution = 1.0;
Lionel Heng's avatar
Lionel Heng committed
125
    if (currentImageryType == GOOGLE_SATELLITE ||
126
            currentImageryType == GOOGLE_MAP) {
Lionel Heng's avatar
Lionel Heng committed
127
        tileResolution = 1.0;
128
        while (tileResolution * 3.0 / 2.0 < 1.0 / zoom) {
Lionel Heng's avatar
Lionel Heng committed
129
130
            tileResolution *= 2.0;
        }
131
        if (tileResolution > 512.0) {
Lionel Heng's avatar
Lionel Heng committed
132
133
            tileResolution = 512.0;
        }
134
    } else if (currentImageryType == SWISSTOPO_SATELLITE) {
Lionel Heng's avatar
Lionel Heng committed
135
136
137
        tileResolution = 0.25;
    }

138
139
    int minTileX, minTileY, maxTileX, maxTileY;
    int zoomLevel;
Lionel Heng's avatar
Lionel Heng committed
140
141

    tileBounds(tileResolution,
142
143
144
145
               xOrigin - windowWidth / 2.0 / zoom * 1.5,
               yOrigin - windowHeight / 2.0 / zoom * 1.5,
               xOrigin + windowWidth / 2.0 / zoom * 1.5,
               yOrigin + windowHeight / 2.0 / zoom * 1.5, utmZone,
Lionel Heng's avatar
Lionel Heng committed
146
147
               minTileX, minTileY, maxTileX, maxTileY, zoomLevel);

148
149
    for (int r = minTileY; r <= maxTileY; ++r) {
        for (int c = minTileX; c <= maxTileX; ++c) {
Lionel Heng's avatar
Lionel Heng committed
150
151
152
153
154
155
            QString tileURL = getTileLocation(c, r, zoomLevel, tileResolution);

            double x1, y1, x2, y2, x3, y3, x4, y4;
            imageBounds(c, r, tileResolution, x1, y1, x2, y2, x3, y3, x4, y4);

            TexturePtr t = textureCache->get(tileURL);
156
            if (!t.isNull()) {
157
158
159
160
                addDrawable(t->draw(y1 - yOffset, x1 - xOffset,
                                    y2 - yOffset, x2 - xOffset,
                                    y3 - yOffset, x3 - xOffset,
                                    y4 - yOffset, x4 - xOffset,
161
                                    zOffset,
162
                                    true));
Lionel Heng's avatar
Lionel Heng committed
163
164
165
166
167
168
169
170
            }
        }
    }
}

void
Imagery::prefetch3D(double radius, double tileResolution,
                    double xOrigin, double yOrigin,
171
                    const QString& utmZone)
Lionel Heng's avatar
Lionel Heng committed
172
{
173
    if (currentImageryType == BLANK_MAP) {
174
175
176
        return;
    }

177
178
    int minTileX, minTileY, maxTileX, maxTileY;
    int zoomLevel;
Lionel Heng's avatar
Lionel Heng committed
179
180

    tileBounds(tileResolution,
181
182
               xOrigin - radius, yOrigin - radius,
               xOrigin + radius, yOrigin + radius, utmZone,
Lionel Heng's avatar
Lionel Heng committed
183
184
               minTileX, minTileY, maxTileX, maxTileY, zoomLevel);

185
186
    for (int r = minTileY; r <= maxTileY; ++r) {
        for (int c = minTileX; c <= maxTileX; ++c) {
Lionel Heng's avatar
Lionel Heng committed
187
188
            QString url = getTileLocation(c, r, zoomLevel, tileResolution);

189
            TexturePtr t = textureCache->get(url);
Lionel Heng's avatar
Lionel Heng committed
190
191
192
193
194
195
196
        }
    }
}

void
Imagery::draw3D(double radius, double tileResolution,
                double xOrigin, double yOrigin,
197
                double xOffset, double yOffset, double zOffset,
198
                const QString& utmZone)
Lionel Heng's avatar
Lionel Heng committed
199
{
200
    if (getNumDrawables() > 0) {
201
202
203
        removeDrawables(0, getNumDrawables());
    }

204
    if (currentImageryType == BLANK_MAP) {
205
206
207
        return;
    }

208
209
    int minTileX, minTileY, maxTileX, maxTileY;
    int zoomLevel;
Lionel Heng's avatar
Lionel Heng committed
210
211

    tileBounds(tileResolution,
212
213
               xOrigin - radius, yOrigin - radius,
               xOrigin + radius, yOrigin + radius, utmZone,
Lionel Heng's avatar
Lionel Heng committed
214
215
               minTileX, minTileY, maxTileX, maxTileY, zoomLevel);

216
217
    for (int r = minTileY; r <= maxTileY; ++r) {
        for (int c = minTileX; c <= maxTileX; ++c) {
Lionel Heng's avatar
Lionel Heng committed
218
219
220
221
222
            QString tileURL = getTileLocation(c, r, zoomLevel, tileResolution);

            double x1, y1, x2, y2, x3, y3, x4, y4;
            imageBounds(c, r, tileResolution, x1, y1, x2, y2, x3, y3, x4, y4);

223
            TexturePtr t = textureCache->get(tileURL);
Lionel Heng's avatar
Lionel Heng committed
224

225
            if (!t.isNull()) {
226
227
228
229
                addDrawable(t->draw(y1 - yOffset, x1 - xOffset,
                                    y2 - yOffset, x2 - xOffset,
                                    y3 - yOffset, x3 - xOffset,
                                    y4 - yOffset, x4 - xOffset,
230
                                    zOffset,
231
                                    true));
Lionel Heng's avatar
Lionel Heng committed
232
233
234
235
236
237
238
239
240
241
242
243
244
245
            }
        }
    }
}

bool
Imagery::update(void)
{
    textureCache->sync();

    return true;
}

void
246
Imagery::imageBounds(int tileX, int tileY, double tileResolution,
Lionel Heng's avatar
Lionel Heng committed
247
248
249
250
                     double& x1, double& y1, double& x2, double& y2,
                     double& x3, double& y3, double& x4, double& y4) const
{
    if (currentImageryType == GOOGLE_MAP ||
251
            currentImageryType == GOOGLE_SATELLITE) {
252
253
        int zoomLevel = MAX_ZOOM_LEVEL - static_cast<int>(rint(log2(tileResolution)));
        int numTiles = static_cast<int>(exp2(static_cast<double>(zoomLevel)));
Lionel Heng's avatar
Lionel Heng committed
254
255
256
257
258
259
260
261
262
263
264
265

        double lon1 = tileXToLongitude(tileX, numTiles);
        double lon2 = tileXToLongitude(tileX + 1, numTiles);

        double lat1 = tileYToLatitude(tileY, numTiles);
        double lat2 = tileYToLatitude(tileY + 1, numTiles);

        QString utmZone;
        LLtoUTM(lat1, lon1, x1, y1, utmZone);
        LLtoUTM(lat1, lon2, x2, y2, utmZone);
        LLtoUTM(lat2, lon2, x3, y3, utmZone);
        LLtoUTM(lat2, lon1, x4, y4, utmZone);
266
    } else if (currentImageryType == SWISSTOPO_SATELLITE) {
Lionel Heng's avatar
Lionel Heng committed
267
268
269
270
271
272
        double utmMultiplier = tileResolution * 200.0;
        double minX = tileX * utmMultiplier;
        double maxX = minX + utmMultiplier;
        double minY = tileY * utmMultiplier;
        double maxY = minY + utmMultiplier;

273
274
275
276
277
278
279
280
        x1 = maxX;
        y1 = minY;
        x2 = maxX;
        y2 = maxY;
        x3 = minX;
        y3 = maxY;
        x4 = minX;
        y4 = minY;
Lionel Heng's avatar
Lionel Heng committed
281
282
283
284
285
286
287
    }
}

void
Imagery::tileBounds(double tileResolution,
                    double minUtmX, double minUtmY,
                    double maxUtmX, double maxUtmY, const QString& utmZone,
288
289
290
                    int& minTileX, int& minTileY,
                    int& maxTileX, int& maxTileY,
                    int& zoomLevel) const
Lionel Heng's avatar
Lionel Heng committed
291
292
293
{
    double centerUtmX = (maxUtmX - minUtmX) / 2.0 + minUtmX;
    double centerUtmY = (maxUtmY - minUtmY) / 2.0 + minUtmY;
294
    int centerTileX, centerTileY;
Lionel Heng's avatar
Lionel Heng committed
295
296

    if (currentImageryType == GOOGLE_MAP ||
297
            currentImageryType == GOOGLE_SATELLITE) {
Lionel Heng's avatar
Lionel Heng committed
298
299
300
301
302
303
        UTMtoTile(minUtmX, minUtmY, utmZone, tileResolution,
                  minTileX, maxTileY, zoomLevel);
        UTMtoTile(centerUtmX, centerUtmY, utmZone, tileResolution,
                  centerTileX, centerTileY, zoomLevel);
        UTMtoTile(maxUtmX, maxUtmY, utmZone, tileResolution,
                  maxTileX, minTileY, zoomLevel);
304
    } else if (currentImageryType == SWISSTOPO_SATELLITE) {
Lionel Heng's avatar
Lionel Heng committed
305
306
        double utmMultiplier = tileResolution * 200;

307
308
309
310
311
312
        minTileX = static_cast<int>(rint(minUtmX / utmMultiplier));
        minTileY = static_cast<int>(rint(minUtmY / utmMultiplier));
        centerTileX = static_cast<int>(rint(centerUtmX / utmMultiplier));
        centerTileY = static_cast<int>(rint(centerUtmY / utmMultiplier));
        maxTileX = static_cast<int>(rint(maxUtmX / utmMultiplier));
        maxTileY = static_cast<int>(rint(maxUtmY / utmMultiplier));
Lionel Heng's avatar
Lionel Heng committed
313
314
    }

315
    if (maxTileX - minTileX + 1 > 14) {
Lionel Heng's avatar
Lionel Heng committed
316
317
318
        minTileX = centerTileX - 7;
        maxTileX = centerTileX + 6;
    }
319
    if (maxTileY - minTileY + 1 > 14) {
Lionel Heng's avatar
Lionel Heng committed
320
321
322
323
324
325
        minTileY = centerTileY - 7;
        maxTileY = centerTileY + 6;
    }
}

double
326
Imagery::tileXToLongitude(int tileX, int numTiles) const
Lionel Heng's avatar
Lionel Heng committed
327
328
329
330
331
332
{
    return 360.0 * (static_cast<double>(tileX)
                    / static_cast<double>(numTiles)) - 180.0;
}

double
333
Imagery::tileYToLatitude(int tileY, int numTiles) const
Lionel Heng's avatar
Lionel Heng committed
334
335
{
    double unnormalizedRad =
336
337
        (static_cast<double>(tileY) / static_cast<double>(numTiles))
        * 2.0 * M_PI - M_PI;
Lionel Heng's avatar
Lionel Heng committed
338
339
340
341
    double rad = 2.0 * atan(exp(unnormalizedRad)) - M_PI / 2.0;
    return -rad * 180.0 / M_PI;
}

342
343
int
Imagery::longitudeToTileX(double longitude, int numTiles) const
Lionel Heng's avatar
Lionel Heng committed
344
{
345
    return static_cast<int>((longitude / 180.0 + 1.0) / 2.0 * numTiles);
Lionel Heng's avatar
Lionel Heng committed
346
347
}

348
349
int
Imagery::latitudeToTileY(double latitude, int numTiles) const
Lionel Heng's avatar
Lionel Heng committed
350
351
352
{
    double rad = latitude * M_PI / 180.0;
    double normalizedRad = -log(tan(rad) + 1.0 / cos(rad));
353
    return static_cast<int>((normalizedRad + M_PI)
354
                            / (2.0 * M_PI) * numTiles);
Lionel Heng's avatar
Lionel Heng committed
355
356
357
358
}

void
Imagery::UTMtoTile(double northing, double easting, const QString& utmZone,
359
360
                   double tileResolution, int& tileX, int& tileY,
                   int& zoomLevel) const
Lionel Heng's avatar
Lionel Heng committed
361
362
363
364
365
{
    double latitude, longitude;

    UTMtoLL(northing, easting, utmZone, latitude, longitude);

366
367
    zoomLevel = MAX_ZOOM_LEVEL - static_cast<int>(rint(log2(tileResolution)));
    int numTiles = static_cast<int>(exp2(static_cast<double>(zoomLevel)));
Lionel Heng's avatar
Lionel Heng committed
368
369
370
371
372
373

    tileX = longitudeToTileX(longitude, numTiles);
    tileY = latitudeToTileY(latitude, numTiles);
}

QChar
374
Imagery::UTMLetterDesignator(double latitude)
Lionel Heng's avatar
Lionel Heng committed
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
{
    // This routine determines the correct UTM letter designator for the given latitude
    // returns 'Z' if latitude is outside the UTM limits of 84N to 80S
    // Written by Chuck Gantz- chuck.gantz@globalstar.com
    char letterDesignator;

    if ((84.0 >= latitude) && (latitude >= 72.0)) letterDesignator = 'X';
    else if ((72.0 > latitude) && (latitude >= 64.0)) letterDesignator = 'W';
    else if ((64.0 > latitude) && (latitude >= 56.0)) letterDesignator = 'V';
    else if ((56.0 > latitude) && (latitude >= 48.0)) letterDesignator = 'U';
    else if ((48.0 > latitude) && (latitude >= 40.0)) letterDesignator = 'T';
    else if ((40.0 > latitude) && (latitude >= 32.0)) letterDesignator = 'S';
    else if ((32.0 > latitude) && (latitude >= 24.0)) letterDesignator = 'R';
    else if ((24.0 > latitude) && (latitude >= 16.0)) letterDesignator = 'Q';
    else if ((16.0 > latitude) && (latitude >= 8.0)) letterDesignator = 'P';
    else if (( 8.0 > latitude) && (latitude >= 0.0)) letterDesignator = 'N';
    else if (( 0.0 > latitude) && (latitude >= -8.0)) letterDesignator = 'M';
    else if ((-8.0 > latitude) && (latitude >= -16.0)) letterDesignator = 'L';
    else if ((-16.0 > latitude) && (latitude >= -24.0)) letterDesignator = 'K';
    else if ((-24.0 > latitude) && (latitude >= -32.0)) letterDesignator = 'J';
    else if ((-32.0 > latitude) && (latitude >= -40.0)) letterDesignator = 'H';
    else if ((-40.0 > latitude) && (latitude >= -48.0)) letterDesignator = 'G';
    else if ((-48.0 > latitude) && (latitude >= -56.0)) letterDesignator = 'F';
    else if ((-56.0 > latitude) && (latitude >= -64.0)) letterDesignator = 'E';
    else if ((-64.0 > latitude) && (latitude >= -72.0)) letterDesignator = 'D';
    else if ((-72.0 > latitude) && (latitude >= -80.0)) letterDesignator = 'C';
    else letterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits

    return letterDesignator;
}

void
Imagery::LLtoUTM(double latitude, double longitude,
                 double& utmNorthing, double& utmEasting,
409
                 QString& utmZone)
Lionel Heng's avatar
Lionel Heng committed
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
{
    // converts lat/long to UTM coords.  Equations from USGS Bulletin 1532
    // East Longitudes are positive, West longitudes are negative.
    // North latitudes are positive, South latitudes are negative
    // Lat and Long are in decimal degrees
    // Written by Chuck Gantz- chuck.gantz@globalstar.com

    double k0 = 0.9996;

    double LongOrigin;
    double eccPrimeSquared;
    double N, T, C, A, M;

    double LatRad = latitude * M_PI / 180.0;
    double LongRad = longitude * M_PI / 180.0;
    double LongOriginRad;

427
    int ZoneNumber = static_cast<int>((longitude + 180.0) / 6.0) + 1;
Lionel Heng's avatar
Lionel Heng committed
428
429

    if (latitude >= 56.0 && latitude < 64.0 &&
430
            longitude >= 3.0 && longitude < 12.0) {
Lionel Heng's avatar
Lionel Heng committed
431
432
433
434
        ZoneNumber = 32;
    }

    // Special zones for Svalbard
435
    if (latitude >= 72.0 && latitude < 84.0) {
Lionel Heng's avatar
Lionel Heng committed
436
437
438
439
        if (     longitude >= 0.0  && longitude <  9.0) ZoneNumber = 31;
        else if (longitude >= 9.0  && longitude < 21.0) ZoneNumber = 33;
        else if (longitude >= 21.0 && longitude < 33.0) ZoneNumber = 35;
        else if (longitude >= 33.0 && longitude < 42.0) ZoneNumber = 37;
440
    }
Lionel Heng's avatar
Lionel Heng committed
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
    LongOrigin = static_cast<double>((ZoneNumber - 1) * 6 - 180 + 3);  //+3 puts origin in middle of zone
    LongOriginRad = LongOrigin * M_PI / 180.0;

    // compute the UTM Zone from the latitude and longitude
    utmZone = QString("%1%2").arg(ZoneNumber).arg(UTMLetterDesignator(latitude));

    eccPrimeSquared = WGS84_ECCSQ / (1.0 - WGS84_ECCSQ);

    N = WGS84_A / sqrt(1.0f - WGS84_ECCSQ * sin(LatRad) * sin(LatRad));
    T = tan(LatRad) * tan(LatRad);
    C = eccPrimeSquared * cos(LatRad) * cos(LatRad);
    A = cos(LatRad) * (LongRad - LongOriginRad);

    M = WGS84_A * ((1.0 - WGS84_ECCSQ / 4.0
                    - 3.0 * WGS84_ECCSQ * WGS84_ECCSQ / 64.0
                    - 5.0 * WGS84_ECCSQ * WGS84_ECCSQ * WGS84_ECCSQ / 256.0)
                   * LatRad
                   - (3.0 * WGS84_ECCSQ / 8.0
                      + 3.0 * WGS84_ECCSQ * WGS84_ECCSQ / 32.0
                      + 45.0 * WGS84_ECCSQ * WGS84_ECCSQ * WGS84_ECCSQ / 1024.0)
                   * sin(2.0 * LatRad)
                   + (15.0 * WGS84_ECCSQ * WGS84_ECCSQ / 256.0
                      + 45.0 * WGS84_ECCSQ * WGS84_ECCSQ * WGS84_ECCSQ / 1024.0)
                   * sin(4.0 * LatRad)
                   - (35.0 * WGS84_ECCSQ * WGS84_ECCSQ * WGS84_ECCSQ / 3072.0)
                   * sin(6.0 * LatRad));

    utmEasting = k0 * N * (A + (1.0 - T + C) * A * A * A / 6.0
                           + (5.0 - 18.0 * T + T * T + 72.0 * C
                              - 58.0 * eccPrimeSquared)
                           * A * A * A * A * A / 120.0)
                 + 500000.0;

    utmNorthing = k0 * (M + N * tan(LatRad) *
                        (A * A / 2.0 +
                         (5.0 - T + 9.0 * C + 4.0 * C * C) * A * A * A * A / 24.0
                         + (61.0 - 58.0 * T + T * T + 600.0 * C
                            - 330.0 * eccPrimeSquared)
                         * A * A * A * A * A * A / 720.0));
480
    if (latitude < 0.0) {
Lionel Heng's avatar
Lionel Heng committed
481
482
483
484
485
486
        utmNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
    }
}

void
Imagery::UTMtoLL(double utmNorthing, double utmEasting, const QString& utmZone,
487
                 double& latitude, double& longitude)
Lionel Heng's avatar
Lionel Heng committed
488
489
490
491
492
493
494
495
496
497
498
499
500
501
{
    // converts UTM coords to lat/long.  Equations from USGS Bulletin 1532
    // East Longitudes are positive, West longitudes are negative.
    // North latitudes are positive, South latitudes are negative
    // Lat and Long are in decimal degrees.
    // Written by Chuck Gantz- chuck.gantz@globalstar.com

    double k0 = 0.9996;
    double eccPrimeSquared;
    double e1 = (1.0 - sqrt(1.0 - WGS84_ECCSQ)) / (1.0 + sqrt(1.0 - WGS84_ECCSQ));
    double N1, T1, C1, R1, D, M;
    double LongOrigin;
    double mu, phi1, phi1Rad;
    double x, y;
502
    int ZoneNumber;
Lionel Heng's avatar
Lionel Heng committed
503
504
505
506
507
508
509
510
    char ZoneLetter;
    bool NorthernHemisphere;

    x = utmEasting - 500000.0; //remove 500,000 meter offset for longitude
    y = utmNorthing;

    std::istringstream iss(utmZone.toStdString());
    iss >> ZoneNumber >> ZoneLetter;
511
    if ((ZoneLetter - 'N') >= 0) {
Lionel Heng's avatar
Lionel Heng committed
512
        NorthernHemisphere = true;//point is in northern hemisphere
513
    } else {
Lionel Heng's avatar
Lionel Heng committed
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
        NorthernHemisphere = false;//point is in southern hemisphere
        y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
    }

    LongOrigin = (ZoneNumber - 1.0) * 6.0 - 180.0 + 3.0;  //+3 puts origin in middle of zone

    eccPrimeSquared = WGS84_ECCSQ / (1.0 - WGS84_ECCSQ);

    M = y / k0;
    mu = M / (WGS84_A * (1.0 - WGS84_ECCSQ / 4.0
                         - 3.0 * WGS84_ECCSQ * WGS84_ECCSQ / 64.0
                         - 5.0 * WGS84_ECCSQ * WGS84_ECCSQ * WGS84_ECCSQ / 256.0));

    phi1Rad = mu + (3.0 * e1 / 2.0 - 27.0 * e1 * e1 * e1 / 32.0) * sin(2.0 * mu)
              + (21.0 * e1 * e1 / 16.0 - 55.0 * e1 * e1 * e1 * e1 / 32.0)
              * sin(4.0 * mu)
              + (151.0 * e1 * e1 * e1 / 96.0) * sin(6.0 * mu);
    phi1 = phi1Rad / M_PI * 180.0;

    N1 = WGS84_A / sqrt(1.0 - WGS84_ECCSQ * sin(phi1Rad) * sin(phi1Rad));
    T1 = tan(phi1Rad) * tan(phi1Rad);
    C1 = eccPrimeSquared * cos(phi1Rad) * cos(phi1Rad);
    R1 = WGS84_A * (1.0 - WGS84_ECCSQ) /
         pow(1.0 - WGS84_ECCSQ * sin(phi1Rad) * sin(phi1Rad), 1.5);
    D = x / (N1 * k0);

    latitude = phi1Rad - (N1 * tan(phi1Rad) / R1)
               * (D * D / 2.0 - (5.0 + 3.0 * T1 + 10.0 * C1 - 4.0 * C1 * C1
                                 - 9.0 * eccPrimeSquared) * D * D * D * D / 24.0
                  + (61.0 + 90.0 * T1 + 298.0 * C1 + 45.0 * T1 * T1
                     - 252.0 * eccPrimeSquared - 3.0 * C1 * C1)
                  * D * D * D * D * D * D / 720.0);
    latitude *= 180.0 / M_PI;

    longitude = (D - (1.0 + 2.0 * T1 + C1) * D * D * D / 6.0
                 + (5.0 - 2.0 * C1 + 28.0 * T1 - 3.0 * C1 * C1
                    + 8.0 * eccPrimeSquared + 24.0 * T1 * T1)
                 * D * D * D * D * D / 120.0) / cos(phi1Rad);
    longitude = LongOrigin + longitude / M_PI * 180.0;
}
554

Lionel Heng's avatar
Lionel Heng committed
555
QString
556
Imagery::getTileLocation(int tileX, int tileY, int zoomLevel,
Lionel Heng's avatar
Lionel Heng committed
557
558
559
560
                         double tileResolution) const
{
    std::ostringstream oss;

561
    switch (currentImageryType) {
Lionel Heng's avatar
Lionel Heng committed
562
563
564
565
566
567
568
569
    case GOOGLE_MAP:
        oss << "http://mt0.google.com/vt/lyrs=m@120&x=" << tileX
            << "&y=" << tileY << "&z=" << zoomLevel;
        break;
    case GOOGLE_SATELLITE:
        oss << "http://khm.google.com/vt/lbw/lyrs=y&x=" << tileX
            << "&y=" << tileY << "&z=" << zoomLevel;
        break;
570
    case SWISSTOPO_SATELLITE:
Lionel Heng's avatar
Lionel Heng committed
571
572
        oss << "../map/eth_zurich_swissimage_025/200/color/" << tileY
            << "/tile-";
573
        if (tileResolution < 1.0) {
Lionel Heng's avatar
Lionel Heng committed
574
            oss << std::fixed << std::setprecision(2) << tileResolution;
575
        } else {
576
            oss << static_cast<int>(rint(tileResolution));
Lionel Heng's avatar
Lionel Heng committed
577
578
579
        }
        oss << "-" << tileY << "-" << tileX << ".jpg";
    default:
580
    {};
Lionel Heng's avatar
Lionel Heng committed
581
582
583
584
585
586
    }

    QString url(oss.str().c_str());

    return url;
}