Imagery.cc 20.8 KB
Newer Older
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
/*=====================================================================

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>
35
#include <cstdio>
36 37 38 39 40 41
#include <iomanip>
#include <sstream>

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

42
const int MAX_ZOOM_LEVEL = 20;
43 44

Imagery::Imagery()
45
 : mTextureCache(new TextureCache(1000))
46 47 48 49
 , mImageryType(Imagery::BLANK_MAP)
 , mXOffset(0.0)
 , mYOffset(0.0)
 , mZOffset(0.0)
50
{
51
    setCullingActive(false);
52 53
}

54
Imagery::Type
55 56
Imagery::getImageryType(void) const
{
57
    return mImageryType;
58 59
}

60
void
61
Imagery::setImageryType(Imagery::Type type)
62
{
63
    mImageryType = type;
64 65 66
}

void
67
Imagery::setOffset(double xOffset, double yOffset, double zOffset)
68
{
69 70 71 72 73 74 75 76 77
    mXOffset = xOffset;
    mYOffset = yOffset;
    mZOffset = zOffset;
}

void
Imagery::setPath(const QString &path)
{
    mImageryPath = path.toStdString();
78 79 80 81 82 83 84
}

void
Imagery::prefetch2D(double windowWidth, double windowHeight,
                    double zoom, double xOrigin, double yOrigin,
                    const QString& utmZone)
{
85 86
    if (mImageryType == BLANK_MAP)
    {
87 88 89
        return;
    }

90
    double tileResolution = 1.0;
91 92 93
    if (mImageryType == GOOGLE_SATELLITE ||
        mImageryType == GOOGLE_MAP)
    {
94
        tileResolution = 1.0;
95 96
        while (tileResolution * 3.0 / 2.0 < 1.0 / zoom)
        {
97 98
            tileResolution *= 2.0;
        }
99 100
        if (tileResolution > 512.0)
        {
101 102
            tileResolution = 512.0;
        }
103 104 105
    }
    else if (mImageryType == OFFLINE_SATELLITE)
    {
106 107 108
        tileResolution = 0.25;
    }

109 110
    int minTileX, minTileY, maxTileX, maxTileY;
    int zoomLevel;
111 112

    tileBounds(tileResolution,
113 114 115 116
               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,
117 118
               minTileX, minTileY, maxTileX, maxTileY, zoomLevel);

119 120 121 122
    for (int r = minTileY; r <= maxTileY; ++r)
    {
        for (int c = minTileX; c <= maxTileX; ++c)
        {
123 124
            QString url = getTileLocation(c, r, zoomLevel, tileResolution);

125
            TexturePtr t = mTextureCache->get(url);
126 127 128 129 130 131 132 133 134
        }
    }
}

void
Imagery::draw2D(double windowWidth, double windowHeight,
                double zoom, double xOrigin, double yOrigin,
                const QString& utmZone)
{
135 136
    if (getNumDrawables() > 0)
    {
137 138 139
        removeDrawables(0, getNumDrawables());
    }

140 141
    if (mImageryType == BLANK_MAP)
    {
142 143 144
        return;
    }

145
    double tileResolution = 1.0;
146 147 148
    if (mImageryType == GOOGLE_SATELLITE ||
        mImageryType == GOOGLE_MAP)
    {
149
        tileResolution = 1.0;
150 151
        while (tileResolution * 3.0 / 2.0 < 1.0 / zoom)
        {
152 153
            tileResolution *= 2.0;
        }
154 155
        if (tileResolution > 512.0)
        {
156 157
            tileResolution = 512.0;
        }
158 159 160
    }
    else if (mImageryType == OFFLINE_SATELLITE)
    {
161 162 163
        tileResolution = 0.25;
    }

164 165
    int minTileX, minTileY, maxTileX, maxTileY;
    int zoomLevel;
166 167

    tileBounds(tileResolution,
168 169 170 171
               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,
172 173
               minTileX, minTileY, maxTileX, maxTileY, zoomLevel);

174 175 176 177
    for (int r = minTileY; r <= maxTileY; ++r)
    {
        for (int c = minTileX; c <= maxTileX; ++c)
        {
178 179 180 181 182
            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);

183 184 185 186 187 188 189 190
            TexturePtr t = mTextureCache->get(tileURL);
            if (!t.isNull())
            {
                addDrawable(t->draw(y1, x1,
                                    y2, x2,
                                    y3, x3,
                                    y4, x4,
                                    - mZOffset,
191
                                    true));
192 193 194 195 196 197 198 199
            }
        }
    }
}

void
Imagery::prefetch3D(double radius, double tileResolution,
                    double xOrigin, double yOrigin,
200
                    const QString& utmZone)
201
{
202 203
    if (mImageryType == BLANK_MAP)
    {
204 205 206
        return;
    }

207 208
    int minTileX, minTileY, maxTileX, maxTileY;
    int zoomLevel;
209 210

    tileBounds(tileResolution,
211 212
               xOrigin + mXOffset - radius, yOrigin + mYOffset - radius,
               xOrigin + mXOffset + radius, yOrigin + mYOffset + radius, utmZone,
213 214
               minTileX, minTileY, maxTileX, maxTileY, zoomLevel);

215 216 217 218
    for (int r = minTileY; r <= maxTileY; ++r)
    {
        for (int c = minTileX; c <= maxTileX; ++c)
        {
219 220
            QString url = getTileLocation(c, r, zoomLevel, tileResolution);

221
            TexturePtr t = mTextureCache->get(url);
222 223 224 225 226 227 228
        }
    }
}

void
Imagery::draw3D(double radius, double tileResolution,
                double xOrigin, double yOrigin,
229
                double xOffset, double yOffset,
230
                const QString& utmZone)
231
{
232 233
    if (getNumDrawables() > 0)
    {
234 235 236
        removeDrawables(0, getNumDrawables());
    }

237 238
    if (mImageryType == BLANK_MAP)
    {
239 240 241
        return;
    }

242 243
    int minTileX, minTileY, maxTileX, maxTileY;
    int zoomLevel;
244 245

    tileBounds(tileResolution,
246 247
               xOrigin + mXOffset - radius, yOrigin + mYOffset - radius,
               xOrigin + mXOffset + radius, yOrigin + mYOffset + radius, utmZone,
248 249
               minTileX, minTileY, maxTileX, maxTileY, zoomLevel);

250 251 252 253
    for (int r = minTileY; r <= maxTileY; ++r)
    {
        for (int c = minTileX; c <= maxTileX; ++c)
        {
254 255 256 257 258
            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);

259
            TexturePtr t = mTextureCache->get(tileURL);
260

261 262 263 264 265 266 267
            if (!t.isNull())
            {
                addDrawable(t->draw(y1 - mYOffset + yOffset, x1 - mXOffset + xOffset,
                                    y2 - mYOffset + yOffset, x2 - mXOffset + xOffset,
                                    y3 - mYOffset + yOffset, x3 - mXOffset + xOffset,
                                    y4 - mYOffset + yOffset, x4 - mXOffset + xOffset,
                                    - mZOffset,
268
                                    true));
269 270 271 272 273 274 275 276
            }
        }
    }
}

bool
Imagery::update(void)
{
277
    mTextureCache->sync();
278 279 280 281 282

    return true;
}

void
283
Imagery::imageBounds(int tileX, int tileY, double tileResolution,
284 285 286
                     double& x1, double& y1, double& x2, double& y2,
                     double& x3, double& y3, double& x4, double& y4) const
{
287 288 289
    if (mImageryType == GOOGLE_MAP ||
        mImageryType == GOOGLE_SATELLITE)
    {
290 291
        int zoomLevel = MAX_ZOOM_LEVEL - static_cast<int>(rint(log2(tileResolution)));
        int numTiles = static_cast<int>(exp2(static_cast<double>(zoomLevel)));
292 293 294 295 296 297 298 299 300 301 302 303

        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);
304 305 306
    }
    else if (mImageryType == OFFLINE_SATELLITE)
    {
307 308 309 310 311 312
        double utmMultiplier = tileResolution * 200.0;
        double minX = tileX * utmMultiplier;
        double maxX = minX + utmMultiplier;
        double minY = tileY * utmMultiplier;
        double maxY = minY + utmMultiplier;

313 314 315 316 317 318 319 320
        x1 = maxX;
        y1 = minY;
        x2 = maxX;
        y2 = maxY;
        x3 = minX;
        y3 = maxY;
        x4 = minX;
        y4 = minY;
321 322 323 324 325 326 327
    }
}

void
Imagery::tileBounds(double tileResolution,
                    double minUtmX, double minUtmY,
                    double maxUtmX, double maxUtmY, const QString& utmZone,
328 329 330
                    int& minTileX, int& minTileY,
                    int& maxTileX, int& maxTileY,
                    int& zoomLevel) const
331 332 333
{
    double centerUtmX = (maxUtmX - minUtmX) / 2.0 + minUtmX;
    double centerUtmY = (maxUtmY - minUtmY) / 2.0 + minUtmY;
334
    int centerTileX, centerTileY;
335

336 337 338
    if (mImageryType == GOOGLE_MAP ||
        mImageryType == GOOGLE_SATELLITE)
    {
339 340 341 342 343 344
        UTMtoTile(minUtmX, minUtmY, utmZone, tileResolution,
                  minTileX, maxTileY, zoomLevel);
        UTMtoTile(centerUtmX, centerUtmY, utmZone, tileResolution,
                  centerTileX, centerTileY, zoomLevel);
        UTMtoTile(maxUtmX, maxUtmY, utmZone, tileResolution,
                  maxTileX, minTileY, zoomLevel);
345 346 347
    }
    else if (mImageryType == OFFLINE_SATELLITE)
    {
348 349
        double utmMultiplier = tileResolution * 200;

350 351 352 353 354 355
        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));
356 357
    }

358 359
    if (maxTileX - minTileX + 1 > 14)
    {
360 361 362
        minTileX = centerTileX - 7;
        maxTileX = centerTileX + 6;
    }
363 364
    if (maxTileY - minTileY + 1 > 14)
    {
365 366 367 368 369 370
        minTileY = centerTileY - 7;
        maxTileY = centerTileY + 6;
    }
}

double
371
Imagery::tileXToLongitude(int tileX, int numTiles) const
372 373 374 375 376 377
{
    return 360.0 * (static_cast<double>(tileX)
                    / static_cast<double>(numTiles)) - 180.0;
}

double
378
Imagery::tileYToLatitude(int tileY, int numTiles) const
379 380
{
    double unnormalizedRad =
381 382
        (static_cast<double>(tileY) / static_cast<double>(numTiles))
        * 2.0 * M_PI - M_PI;
383 384 385 386
    double rad = 2.0 * atan(exp(unnormalizedRad)) - M_PI / 2.0;
    return -rad * 180.0 / M_PI;
}

387 388
int
Imagery::longitudeToTileX(double longitude, int numTiles) const
389
{
390
    return static_cast<int>((longitude / 180.0 + 1.0) / 2.0 * numTiles);
391 392
}

393 394
int
Imagery::latitudeToTileY(double latitude, int numTiles) const
395 396 397
{
    double rad = latitude * M_PI / 180.0;
    double normalizedRad = -log(tan(rad) + 1.0 / cos(rad));
398
    return static_cast<int>((normalizedRad + M_PI)
399
                            / (2.0 * M_PI) * numTiles);
400 401 402 403
}

void
Imagery::UTMtoTile(double northing, double easting, const QString& utmZone,
404 405
                   double tileResolution, int& tileX, int& tileY,
                   int& zoomLevel) const
406 407 408 409 410
{
    double latitude, longitude;

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

411 412
    zoomLevel = MAX_ZOOM_LEVEL - static_cast<int>(rint(log2(tileResolution)));
    int numTiles = static_cast<int>(exp2(static_cast<double>(zoomLevel)));
413 414 415 416 417 418

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

QChar
419
Imagery::UTMLetterDesignator(double latitude)
420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453
{
    // 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,
454
                 QString& utmZone)
455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471
{
    // 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;

472
    int ZoneNumber = static_cast<int>((longitude + 180.0) / 6.0) + 1;
473 474

    if (latitude >= 56.0 && latitude < 64.0 &&
475
            longitude >= 3.0 && longitude < 12.0) {
476 477 478 479
        ZoneNumber = 32;
    }

    // Special zones for Svalbard
480
    if (latitude >= 72.0 && latitude < 84.0) {
481 482 483 484
        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;
485
    }
486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524
    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));
525
    if (latitude < 0.0) {
526 527 528 529 530 531
        utmNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
    }
}

void
Imagery::UTMtoLL(double utmNorthing, double utmEasting, const QString& utmZone,
532
                 double& latitude, double& longitude)
533 534 535 536 537 538 539 540 541 542 543 544 545 546
{
    // 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;
547
    int ZoneNumber;
548 549 550 551 552 553 554 555
    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;
556
    if ((ZoneLetter - 'N') >= 0) {
557
        NorthernHemisphere = true;//point is in northern hemisphere
558
    } else {
559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598
        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;
}
599

600
QString
601
Imagery::getTileLocation(int tileX, int tileY, int zoomLevel,
602 603 604 605
                         double tileResolution) const
{
    std::ostringstream oss;

606 607
    switch (mImageryType)
    {
608 609 610 611 612 613 614 615
    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;
616 617
    case OFFLINE_SATELLITE:
        oss << mImageryPath << "/200/color/" << tileY
618
            << "/tile-";
619 620
        if (tileResolution < 1.0)
        {
621
            oss << std::fixed << std::setprecision(2) << tileResolution;
622 623 624
        }
        else
        {
625
            oss << static_cast<int>(rint(tileResolution));
626 627 628
        }
        oss << "-" << tileY << "-" << tileX << ".jpg";
    default:
629
    {};
630 631 632 633 634 635
    }

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

    return url;
}