src/map/geotiff.cpp

   1 #include "pcs.h"
   2 #include "tifffile.h"
   3 #include "geotiff.h"
   4
   5
   6 #define TIFF_DOUBLE 12
   7 #define TIFF_SHORT  3
   8 #define TIFF_LONG   4
   9
  10 #define ModelPixelScaleTag            33550
  11 #define ModelTiepointTag              33922
  12 #define ModelTransformationTag        34264
  13 #define GeoKeyDirectoryTag            34735
  14 #define GeoDoubleParamsTag            34736
  15
  16 #define GTModelTypeGeoKey             1024
  17 #define GTRasterTypeGeoKey            1025
  18 #define GeographicTypeGeoKey          2048
  19 #define GeogGeodeticDatumGeoKey       2050
  20 #define GeogPrimeMeridianGeoKey       2051
  21 #define GeogAngularUnitsGeoKey        2054
  22 #define GeogEllipsoidGeoKey           2056
  23 #define GeogAzimuthUnitsGeoKey        2060
  24 #define ProjectedCSTypeGeoKey         3072
  25 #define ProjectionGeoKey              3074
  26 #define ProjCoordTransGeoKey          3075
  27 #define ProjLinearUnitsGeoKey         3076
  28 #define ProjStdParallel1GeoKey        3078
  29 #define ProjStdParallel2GeoKey        3079
  30 #define ProjNatOriginLongGeoKey       3080
  31 #define ProjNatOriginLatGeoKey        3081
  32 #define ProjFalseEastingGeoKey        3082
  33 #define ProjFalseNorthingGeoKey       3083
  34 #define ProjFalseOriginLongGeoKey     3084
  35 #define ProjFalseOriginLatGeoKey      3085
  36 #define ProjFalseOriginEastingGeoKey  3086
  37 #define ProjFalseOriginNorthingGeoKey 3087
  38 #define ProjCenterLongGeoKey          3088
  39 #define ProjCenterLatGeoKey           3089
  40 #define ProjCenterEastingGeoKey       3090
  41 #define ProjCenterNorthingGeoKey      3091
  42 #define ProjScaleAtNatOriginGeoKey    3092
  43 #define ProjScaleAtCenterGeoKey       3093
  44 #define ProjAzimuthAngleGeoKey        3094
  45 #define ProjRectifiedGridAngleGeoKey  3096
  46
  47 #define ModelTypeProjected            1
  48 #define ModelTypeGeographic           2
  49 #define ModelTypeGeocentric           3
  50
  51 #define CT_TransverseMercator         1
  52 #define CT_ObliqueMercator            3
  53 #define CT_Mercator                   7
  54 #define CT_LambertConfConic_2SP       8
  55 #define CT_LambertConfConic_1SP       9
  56 #define CT_LambertAzimEqualArea       10
  57 #define CT_AlbersEqualArea            11
  58
  59
  60 #define IS_SET(map, key) \
  61         ((map).contains(key) && (map).value(key).SHORT != 32767)
  62
  63
  64 typedef struct {
  65         quint16 KeyDirectoryVersion;
  66         quint16 KeyRevision;
  67         quint16 MinorRevision;
  68         quint16 NumberOfKeys;
  69 } Header;
  70
  71 typedef struct {
  72         quint16 KeyID;
  73         quint16 TIFFTagLocation;
  74         quint16 Count;
  75         quint16 ValueOffset;
  76 } KeyEntry;
  77
  78
  79 bool GeoTIFF::readEntry(TIFFFile &file, Ctx &ctx) const
  80 {
  81         quint16 tag;
  82         quint16 type;
  83         quint32 count, offset;
  84
  85         if (!file.readValue(tag))
  86                 return false;
  87         if (!file.readValue(type))
  88                 return false;
  89         if (!file.readValue(count))
  90                 return false;
  91         if (!file.readValue(offset))
  92                 return false;
  93
  94         switch (tag) {
  95                 case ModelPixelScaleTag:
  96                         if (type != TIFF_DOUBLE || count != 3)
  97                                 return false;
  98                         ctx.scale = offset;
  99                         break;
 100                 case ModelTiepointTag:
 101                         if (type != TIFF_DOUBLE || count < 6 || count % 6)
 102                                 return false;
 103                         ctx.tiepoints = offset;
 104                         ctx.tiepointCount = count / 6;
 105                         break;
 106                 case GeoKeyDirectoryTag:
 107                         if (type != TIFF_SHORT)
 108                                 return false;
 109                         ctx.keys = offset;
 110                         break;
 111                 case GeoDoubleParamsTag:
 112                         if (type != TIFF_DOUBLE)
 113                                 return false;
 114                         ctx.values = offset;
 115                         break;
 116                 case ModelTransformationTag:
 117                         if (type != TIFF_DOUBLE || count != 16)
 118                                 return false;
 119                         ctx.matrix = offset;
 120                         break;
 121         }
 122
 123         return true;
 124 }
 125
 126 bool GeoTIFF::readIFD(TIFFFile &file, quint32 &offset, Ctx &ctx) const
 127 {
 128         quint16 count;
 129
 130         if (!file.seek(offset))
 131                 return false;
 132         if (!file.readValue(count))
 133                 return false;
 134
 135         for (quint16 i = 0; i < count; i++)
 136                 if (!readEntry(file, ctx))
 137                         return false;
 138
 139         if (!file.readValue(offset))
 140                 return false;
 141
 142         return true;
 143 }
 144
 145 bool GeoTIFF::readScale(TIFFFile &file, quint32 offset, PointD &scale) const
 146 {
 147         if (!file.seek(offset))
 148                 return false;
 149
 150         if (!file.readValue(scale.rx()))
 151                 return false;
 152         if (!file.readValue(scale.ry()))
 153                 return false;
 154
 155         return true;
 156 }
 157
 158 bool GeoTIFF::readTiepoints(TIFFFile &file, quint32 offset, quint32 count,
 159   QList<ReferencePoint> &points) const
 160 {
 161         double z, pz;
 162         PointD xy, pp;
 163
 164         if (!file.seek(offset))
 165                 return false;
 166
 167         for (quint32 i = 0; i < count; i++) {
 168                 if (!file.readValue(xy.rx()))
 169                         return false;
 170                 if (!file.readValue(xy.ry()))
 171                         return false;
 172                 if (!file.readValue(z))
 173                         return false;
 174
 175                 if (!file.readValue(pp.rx()))
 176                         return false;
 177                 if (!file.readValue(pp.ry()))
 178                         return false;
 179                 if (!file.readValue(pz))
 180                         return false;
 181
 182                 points.append(ReferencePoint(xy, pp));
 183         }
 184
 185         return true;
 186 }
 187
 188 bool GeoTIFF::readMatrix(TIFFFile &file, quint32 offset, double matrix[16]) const
 189 {
 190         if (!file.seek(offset))
 191                 return false;
 192
 193         for (int i = 0; i < 16; i++)
 194                 if (!file.readValue(matrix[i]))
 195                         return false;
 196
 197         return true;
 198 }
 199
 200 bool GeoTIFF::readKeys(TIFFFile &file, Ctx &ctx, QMap<quint16, Value> &kv) const
 201 {
 202         Header header;
 203         KeyEntry entry;
 204         Value value;
 205
 206         if (!file.seek(ctx.keys))
 207                 return false;
 208
 209         if (!file.readValue(header.KeyDirectoryVersion))
 210                 return false;
 211         if (!file.readValue(header.KeyRevision))
 212                 return false;
 213         if (!file.readValue(header.MinorRevision))
 214                 return false;
 215         if (!file.readValue(header.NumberOfKeys))
 216                 return false;
 217
 218         for (int i = 0; i < header.NumberOfKeys; i++) {
 219                 if (!file.readValue(entry.KeyID))
 220                         return false;
 221                 if (!file.readValue(entry.TIFFTagLocation))
 222                         return false;
 223                 if (!file.readValue(entry.Count))
 224                         return false;
 225                 if (!file.readValue(entry.ValueOffset))
 226                         return false;
 227
 228                 switch (entry.KeyID) {
 229                         case GTModelTypeGeoKey:
 230                         case GTRasterTypeGeoKey:
 231                         case GeographicTypeGeoKey:
 232                         case GeogGeodeticDatumGeoKey:
 233                         case GeogPrimeMeridianGeoKey:
 234                         case GeogAngularUnitsGeoKey:
 235                         case GeogAzimuthUnitsGeoKey:
 236                         case ProjectedCSTypeGeoKey:
 237                         case ProjectionGeoKey:
 238                         case ProjCoordTransGeoKey:
 239                         case ProjLinearUnitsGeoKey:
 240                         case GeogEllipsoidGeoKey:
 241                                 if (entry.TIFFTagLocation != 0 || entry.Count != 1)
 242                                         return false;
 243                                 value.SHORT = entry.ValueOffset;
 244                                 kv.insert(entry.KeyID, value);
 245                                 break;
 246                         case ProjScaleAtNatOriginGeoKey:
 247                         case ProjNatOriginLongGeoKey:
 248                         case ProjNatOriginLatGeoKey:
 249                         case ProjFalseEastingGeoKey:
 250                         case ProjFalseNorthingGeoKey:
 251                         case ProjStdParallel1GeoKey:
 252                         case ProjStdParallel2GeoKey:
 253                         case ProjCenterLongGeoKey:
 254                         case ProjCenterLatGeoKey:
 255                         case ProjScaleAtCenterGeoKey:
 256                         case ProjAzimuthAngleGeoKey:
 257                         case ProjRectifiedGridAngleGeoKey:
 258                         case ProjFalseOriginLongGeoKey:
 259                         case ProjFalseOriginLatGeoKey:
 260                         case ProjCenterEastingGeoKey:
 261                         case ProjCenterNorthingGeoKey:
 262                         case ProjFalseOriginEastingGeoKey:
 263                         case ProjFalseOriginNorthingGeoKey:
 264                                 if (!readGeoValue(file, ctx.values, entry.ValueOffset,
 265                                   value.DOUBLE))
 266                                         return false;
 267                                 kv.insert(entry.KeyID, value);
 268                                 break;
 269                         default:
 270                                 break;
 271                 }
 272         }
 273
 274         return true;
 275 }
 276
 277 bool GeoTIFF::readGeoValue(TIFFFile &file, quint32 offset, quint16 index,
 278   double &val) const
 279 {
 280         qint64 pos = file.pos();
 281
 282         if (!file.seek(offset + index * sizeof(double)))
 283                 return false;
 284         if (!file.readValue(val))
 285                 return false;
 286
 287         if (!file.seek(pos))
 288                 return false;
 289
 290         return true;
 291 }
 292
 293 const GCS *GeoTIFF::gcs(QMap<quint16, Value> &kv)
 294 {
 295         const GCS *gcs = 0;
 296
 297         if (IS_SET(kv, GeographicTypeGeoKey)) {
 298                 if (!(gcs = GCS::gcs(kv.value(GeographicTypeGeoKey).SHORT)))
 299                         _errorString = QString("%1: unknown GCS")
 300                           .arg(kv.value(GeographicTypeGeoKey).SHORT);
 301         } else if (IS_SET(kv, GeogGeodeticDatumGeoKey)
 302           || kv.value(GeogEllipsoidGeoKey).SHORT == 7019
 303           || kv.value(GeogEllipsoidGeoKey).SHORT == 7030) {
 304                 int pm = IS_SET(kv, GeogPrimeMeridianGeoKey)
 305                   ? kv.value(GeogPrimeMeridianGeoKey).SHORT : 8901;
 306                 int au = IS_SET(kv, GeogAngularUnitsGeoKey)
 307                   ? kv.value(GeogAngularUnitsGeoKey).SHORT : 9102;
 308
 309                 // If only the ellipsoid is defined and it is GRS80 or WGS84, handle
 310                 // such definition as a WGS84 geodetic datum.
 311                 int gd = IS_SET(kv, GeogGeodeticDatumGeoKey)
 312                   ? kv.value(GeogGeodeticDatumGeoKey).SHORT : 6326;
 313
 314                 if (!(gcs = GCS::gcs(gd, pm, au)))
 315                         _errorString = QString("%1+%2: unknown geodetic datum + prime"
 316                           " meridian combination").arg(gd).arg(pm);
 317         } else
 318                 _errorString = "Can not determine GCS";
 319
 320         return gcs;
 321 }
 322
 323 Projection::Method GeoTIFF::method(QMap<quint16, Value> &kv)
 324 {
 325         if (!IS_SET(kv, ProjCoordTransGeoKey)) {
 326                 _errorString = "Missing coordinate transformation method";
 327                 return Projection::Method();
 328         }
 329
 330         switch (kv.value(ProjCoordTransGeoKey).SHORT) {
 331                 case CT_TransverseMercator:
 332                         return Projection::Method(9807);
 333                 case CT_ObliqueMercator:
 334                         return Projection::Method(9815);
 335                 case CT_Mercator:
 336                         return Projection::Method(9804);
 337                 case CT_LambertConfConic_2SP:
 338                         return Projection::Method(9802);
 339                 case CT_LambertConfConic_1SP:
 340                         return Projection::Method(9801);
 341                 case CT_LambertAzimEqualArea:
 342                         return Projection::Method(9820);
 343                 case CT_AlbersEqualArea:
 344                         return Projection::Method(9822);
 345                 default:
 346                         _errorString = QString("%1: unknown coordinate transformation method")
 347                           .arg(kv.value(ProjCoordTransGeoKey).SHORT);
 348                         return Projection::Method();
 349         }
 350 }
 351
 352 bool GeoTIFF::projectedModel(QMap<quint16, Value> &kv)
 353 {
 354         if (IS_SET(kv, ProjectedCSTypeGeoKey)) {
 355                 const PCS *pcs;
 356                 if (!(pcs = PCS::pcs(kv.value(ProjectedCSTypeGeoKey).SHORT))) {
 357                         _errorString = QString("%1: unknown PCS")
 358                           .arg(kv.value(ProjectedCSTypeGeoKey).SHORT);
 359                         return false;
 360                 }
 361                 _projection = Projection(pcs);
 362         } else if (IS_SET(kv, ProjectionGeoKey)) {
 363                 const GCS *g = gcs(kv);
 364                 if (!g)
 365                         return false;
 366                 const PCS *pcs = PCS::pcs(g, kv.value(ProjectionGeoKey).SHORT);
 367                 if (!pcs) {
 368                         _errorString = QString("%1: unknown projection code")
 369                           .arg(kv.value(GeographicTypeGeoKey).SHORT)
 370                           .arg(kv.value(ProjectionGeoKey).SHORT);
 371                         return false;
 372                 }
 373                 _projection = Projection(pcs);
 374         } else {
 375                 double lat0, lon0, scale, fe, fn, sp1, sp2;
 376
 377                 const GCS *g = gcs(kv);
 378                 if (!g)
 379                         return false;
 380                 Projection::Method m(method(kv));
 381                 if (m.isNull())
 382                         return false;
 383
 384                 AngularUnits au(IS_SET(kv, GeogAngularUnitsGeoKey)
 385                   ? kv.value(GeogAngularUnitsGeoKey).SHORT : 9102);
 386                 AngularUnits zu(IS_SET(kv, GeogAzimuthUnitsGeoKey)
 387                   ? kv.value(GeogAzimuthUnitsGeoKey).SHORT : 9102);
 388                 LinearUnits lu(IS_SET(kv, ProjLinearUnitsGeoKey)
 389                   ? kv.value(ProjLinearUnitsGeoKey).SHORT : 9001);
 390                 if (lu.isNull()) {
 391                         _errorString = QString("%1: unknown projection linear units code")
 392                           .arg(kv.value(ProjLinearUnitsGeoKey).SHORT);
 393                         return false;
 394                 }
 395
 396                 if (kv.contains(ProjNatOriginLatGeoKey))
 397                         lat0 = au.toDegrees(kv.value(ProjNatOriginLatGeoKey).DOUBLE);
 398                 else if (kv.contains(ProjCenterLatGeoKey))
 399                         lat0 = au.toDegrees(kv.value(ProjCenterLatGeoKey).DOUBLE);
 400                 else if (kv.contains(ProjFalseOriginLatGeoKey))
 401                         lat0 = au.toDegrees(kv.value(ProjFalseOriginLatGeoKey).DOUBLE);
 402                 else
 403                         lat0 = NAN;
 404                 if (kv.contains(ProjNatOriginLongGeoKey))
 405                         lon0 = au.toDegrees(kv.value(ProjNatOriginLongGeoKey).DOUBLE);
 406                 else if (kv.contains(ProjCenterLongGeoKey))
 407                         lon0 = au.toDegrees(kv.value(ProjCenterLongGeoKey).DOUBLE);
 408                 else if (kv.contains(ProjFalseOriginLongGeoKey))
 409                         lon0 = au.toDegrees(kv.value(ProjFalseOriginLongGeoKey).DOUBLE);
 410                 else
 411                         lon0 = NAN;
 412                 if (kv.contains(ProjScaleAtNatOriginGeoKey))
 413                         scale = kv.value(ProjScaleAtNatOriginGeoKey).DOUBLE;
 414                 else if (kv.contains(ProjScaleAtCenterGeoKey))
 415                         scale = kv.value(ProjScaleAtCenterGeoKey).DOUBLE;
 416                 else
 417                         scale = NAN;
 418                 if (kv.contains(ProjStdParallel1GeoKey))
 419                         sp1 = au.toDegrees(kv.value(ProjStdParallel1GeoKey).DOUBLE);
 420                 else if (kv.contains(ProjAzimuthAngleGeoKey))
 421                         sp1 = zu.toDegrees(kv.value(ProjAzimuthAngleGeoKey).DOUBLE);
 422                 else
 423                         sp1 = NAN;
 424                 if (kv.contains(ProjStdParallel2GeoKey))
 425                         sp2 = au.toDegrees(kv.value(ProjStdParallel2GeoKey).DOUBLE);
 426                 else if (kv.contains(ProjRectifiedGridAngleGeoKey))
 427                         sp2 = au.toDegrees(kv.value(ProjRectifiedGridAngleGeoKey).DOUBLE);
 428                 else
 429                         sp2 = NAN;
 430                 if (kv.contains(ProjFalseNorthingGeoKey))
 431                         fn = lu.toMeters(kv.value(ProjFalseNorthingGeoKey).DOUBLE);
 432                 else if (kv.contains(ProjCenterNorthingGeoKey))
 433                         fn = lu.toMeters(kv.value(ProjCenterNorthingGeoKey).DOUBLE);
 434                 else if (kv.contains(ProjFalseOriginNorthingGeoKey))
 435                         fn = lu.toMeters(kv.value(ProjFalseOriginNorthingGeoKey).DOUBLE);
 436                 else
 437                         fn = NAN;
 438                 if (kv.contains(ProjFalseEastingGeoKey))
 439                         fe = lu.toMeters(kv.value(ProjFalseEastingGeoKey).DOUBLE);
 440                 else if (kv.contains(ProjCenterEastingGeoKey))
 441                         fe = lu.toMeters(kv.value(ProjCenterEastingGeoKey).DOUBLE);
 442                 else if (kv.contains(ProjFalseOriginEastingGeoKey))
 443                         fe = lu.toMeters(kv.value(ProjFalseOriginEastingGeoKey).DOUBLE);
 444                 else
 445                         fe = NAN;
 446
 447                 Projection::Setup setup(lat0, lon0, scale, fe, fn, sp1, sp2);
 448                 PCS pcs(g, m, setup, lu, CoordinateSystem());
 449                 _projection = Projection(&pcs);
 450         }
 451
 452         return true;
 453 }
 454
 455 bool GeoTIFF::geographicModel(QMap<quint16, Value> &kv)
 456 {
 457         const GCS *g = gcs(kv);
 458         if (!g)
 459                 return false;
 460
 461         _projection = Projection(g);
 462
 463         return true;
 464 }
 465
 466 GeoTIFF::GeoTIFF(const QString &path)
 467 {
 468         quint32 ifd;
 469         QList<ReferencePoint> points;
 470         PointD scale;
 471         QMap<quint16, Value> kv;
 472         Ctx ctx;
 473         TIFFFile file;
 474
 475
 476         file.setFileName(path);
 477         if (!file.open(QIODevice::ReadOnly)) {
 478                 _errorString = QString("Error opening TIFF file: %1")
 479                   .arg(file.errorString());
 480                 return;
 481         }
 482         if (!file.readHeader(ifd)) {
 483                 _errorString = "Invalid TIFF header";
 484                 return;
 485         }
 486
 487         while (ifd) {
 488                 if (!readIFD(file, ifd, ctx)) {
 489                         _errorString = "Invalid IFD";
 490                         return;
 491                 }
 492         }
 493
 494         if (!ctx.keys) {
 495                 _errorString = "Not a GeoTIFF file";
 496                 return;
 497         }
 498
 499         if (ctx.scale) {
 500                 if (!readScale(file, ctx.scale, scale)) {
 501                         _errorString = "Error reading model pixel scale";
 502                         return;
 503                 }
 504         }
 505         if (ctx.tiepoints) {
 506                 if (!readTiepoints(file, ctx.tiepoints, ctx.tiepointCount, points)) {
 507                         _errorString = "Error reading raster->model tiepoint pairs";
 508                         return;
 509                 }
 510         }
 511
 512         if (!readKeys(file, ctx, kv)) {
 513                 _errorString = "Error reading Geo key/value";
 514                 return;
 515         }
 516
 517         switch (kv.value(GTModelTypeGeoKey).SHORT) {
 518                 case ModelTypeProjected:
 519                         if (!projectedModel(kv))
 520                                 return;
 521                         break;
 522                 case ModelTypeGeographic:
 523                         if (!geographicModel(kv))
 524                                 return;
 525                         break;
 526                 case ModelTypeGeocentric:
 527                         _errorString = "Geocentric models are not supported";
 528                         return;
 529                 default:
 530                         _errorString = "Unknown/missing model type";
 531                         return;
 532         }
 533
 534         if (ctx.scale && ctx.tiepoints)
 535                 _transform = Transform(points.first(), scale);
 536         else if (ctx.tiepointCount > 1)
 537                 _transform = Transform(points);
 538         else if (ctx.matrix) {
 539                 double matrix[16];
 540                 if (!readMatrix(file, ctx.matrix, matrix)) {
 541                         _errorString = "Error reading transformation matrix";
 542                         return;
 543                 }
 544                 _transform = Transform(matrix);
 545         } else {
 546                 _errorString = "Incomplete/missing model transformation info";
 547                 return;
 548         }
 549
 550         if (!_transform.isValid())
 551                 _errorString = _transform.errorString();
 552 }