OSSIM - Open Source Software Image Map  Version 1.9.0 (20180803)
ossimGeoTiff.cpp
Go to the documentation of this file.
1 //***************************************************************************
2 // FILE: ossimGeoTiff.cpp
3 //
4 // License: See top level LICENSE.txt file.
5 //
6 // Description:
7 //
8 // Class definition for ossimGeoTiff which is designed to read and hold tag
9 // information.
10 //
11 //***************************************************************************
12 // $Id: ossimGeoTiff.cpp 21024 2012-05-30 08:45:13Z dburken $
13 
15 #include <ossim/base/ossimTrace.h>
16 #include <ossim/base/ossimCommon.h>
39 #include <tiff.h>
40 #include <tiffio.h>
41 #include <xtiffio.h>
42 #include <geotiff.h>
43 #include <geo_normalize.h>
44 #include <geovalues.h>
45 #include <string.h>
46 #include <iomanip>
47 #include <iterator>
48 #include <sstream>
49 #include <cstdlib>
50 
51 static const ossimGeoTiffCoordTransformsLut COORD_TRANS_LUT;
52 static const ossimGeoTiffDatumLut DATUM_LUT;
53 std::mutex ossimGeoTiff::theMutex;
54 
55 #ifdef OSSIM_ID_ENABLED
56 static const char OSSIM_ID[] = "$Id: ossimGeoTiff.cpp 21024 2012-05-30 08:45:13Z dburken $";
57 #endif
58 
59 //---
60 // Static trace for debugging
61 //---
62 static ossimTrace traceDebug("ossimGeoTiff:debug");
63 
64 // Prototype, defined at bottom of this file. ArcMAP 9.2 bug workaround.
66 
67 //---
68 // This was created to remove/hide "libgetiff/geo_normalize.h" in ossimGeoTiff.h.
69 //---
71 {
72 public:
74  : m_defs(0)
75  {}
77  {
78  if ( m_defs )
79  {
80  delete m_defs;
81  m_defs = 0;
82  }
83  }
84  GTIFDefn* m_defs;
85 };
86 
87 //*************************************************************************************************
88 // CONSTRUCTOR
89 //*************************************************************************************************
91  :
92  theTiffPtr(0),
93  theGeoKeyOffset(0),
94  theGeoKeyLength(0),
95  theGeoKeysPresentFlag(false),
96  theZone(0),
97  theHemisphere("N"),
98  theDoubleParamLength(0),
99  theAsciiParamLength(0),
100  theProjectionName("unknown"),
101  theDatumName("unknown"),
102  theWidth(0),
103  theLength(0),
104  theBitsPerSample(0),
105  theModelType(0),
106  theRasterType(ossimGeoTiff::UNDEFINED),
107  theGcsCode(0),
108  theDatumCode(0),
109  theAngularUnits(0),
110  thePcsCode(0),
111  theCoorTransGeoCode(0),
112  theLinearUnitsCode(ossimGeoTiff::UNDEFINED),
113  theStdPar1(0.0),
114  theStdPar2(0.0),
115  theOriginLon(0.0),
116  theOriginLat(0.0),
117  theFalseEasting(0.0),
118  theFalseNorthing(0.0),
119  theScaleFactor(0.0),
120  thePrivateDefinitions(new ossimPrivateGtifDef())
121 {
122 }
123 
124 //*************************************************************************************************
125 // CONSTRUCTOR
126 //*************************************************************************************************
128  :
129  theTiffPtr(0),
130  theGeoKeyOffset(0),
131  theGeoKeyLength(0),
132  theGeoKeysPresentFlag(false),
133  theZone(0),
134  theHemisphere("N"),
135  theDoubleParamLength(0),
136  theAsciiParamLength(0),
137  theProjectionName("unknown"),
138  theDatumName("unknown"),
139  theWidth(0),
140  theLength(0),
141  theBitsPerSample(0),
142  theModelType(0),
143  theRasterType(ossimGeoTiff::UNDEFINED),
144  theGcsCode(0),
145  theDatumCode(0),
146  theAngularUnits(0),
147  thePcsCode(0),
148  theCoorTransGeoCode(0),
149  theLinearUnitsCode(ossimGeoTiff::UNDEFINED),
150  theStdPar1(0.0),
151  theStdPar2(0.0),
152  theOriginLon(0.0),
153  theOriginLat(0.0),
154  theFalseEasting(0.0),
155  theFalseNorthing(0.0),
156  theScaleFactor(0.0),
157  thePrivateDefinitions(new ossimPrivateGtifDef())
158 {
159  if (traceDebug())
160  {
162  << "DEBUG ossimGeoTiff::ossimGeoTiff: Entered..." << std::endl;
163 #ifdef OSSIM_ID_ENABLED
165  << "DEBUG ossimGeoTiff::ossimGeoTiff: OSSIM_ID = "
166  << OSSIM_ID << endl;
167 #endif
168  }
169 
170  if(readTags(file, entryIdx) == false)
171  {
173  if (traceDebug())
174  {
176  << "DEBUG ossimGeoTiff::ossimGeoTiff: "
177  << "Unable to reade tags."
178  << std::endl;
179  }
181  << "FATAL ossimGeoTiff::ossimGeoTiff: "
182  << "Unable to reade tags."
183  << std::endl;
184  }
185  if (traceDebug())
186  {
188  }
189 
190 }
191 
193 {
195  {
196  delete thePrivateDefinitions;
198  }
199  if(theTiffPtr)
200  {
201  XTIFFClose(theTiffPtr);
202  theTiffPtr = 0;
203  }
204 }
205 
207 {
208  int pcsUnits = ossimGeoTiff::UNDEFINED;
209  ossimUnitType units = OSSIM_UNIT_UNKNOWN; // default
210 
211  // Need to instantiate a projection given the pcs code:
213  ossimEpsgProjectionDatabase::instance()->findProjection(pcsCode));
214  if (proj)
215  units = proj->getProjectionUnits();
216  else
218 
219  switch (units)
220  {
221  case OSSIM_METERS:
222  pcsUnits = ossimGeoTiff::LINEAR_METER;
223  break;
224  case OSSIM_FEET:
225  pcsUnits = ossimGeoTiff::LINEAR_FOOT;
226  break;
229  break;
230  default:
231  break;
232  }
233 
234  return pcsUnits;
235 }
236 
237 #define EPSG_CODE_MAX 32767
239  const ossimRefPtr<ossimMapProjectionInfo> projectionInfo,
240  bool imagineNad27Flag)
241 {
242  std::lock_guard<std::mutex> lock(theMutex);
243 
244  const ossimMapProjection* mapProj = projectionInfo->getProjection();
245 
246  if(!mapProj) return false;
247 
248  GTIF* gtif = GTIFNew(tifPtr);
249 
250  // Get some things we need throughout.
251  ossimGpt origin = mapProj->origin();
252  double falseEasting = mapProj->getFalseEasting();
253  double falseNorthing = mapProj->getFalseNorthing();
254 
255  ossimKeywordlist kwl;
256  mapProj->saveState(kwl);
257  const char* stdParallel1 = kwl.find(ossimKeywordNames::STD_PARALLEL_1_KW);
258  const char* stdParallel2 = kwl.find(ossimKeywordNames::STD_PARALLEL_2_KW);
259  const char* scaleFactor = kwl.find(ossimKeywordNames::SCALE_FACTOR_KW);
260 
261  bool gcsTypeSet = false;
262 
263  //---
264  // Since using a pcs code is the easiest way to go, look for that first.
265  //---
266  ossim_int16 pcsCode = mapProj->getPcsCode();
267 
268  //---
269  // Get the units now. If user has overriden pcs units then go user defined
270  // projection by setting pcs code to 0.
271  //---
272  ossimString projName = mapProj->getClassName();
273 
274  int units = ossimGeoTiff::UNDEFINED;
275  if(mapProj->isGeographic())
277  else
278  units = getPcsUnitType(pcsCode);
279  if (units == UNDEFINED)
280  units = LINEAR_METER;
281 
282  if (pcsCode)
283  {
284  if ((units==LINEAR_FOOT_US_SURVEY) || (units==LINEAR_FOOT))
285  {
286  // ArcMap 9.2 bug workaround (originally implemented by ESH 2008, reworked by OLK 04/2010
287  ossim_uint16 meter_pcs = getMetersEquivalentHarnCode(pcsCode);
288  if (meter_pcs)
289  pcsCode = meter_pcs;
290  }
291 
292  //int gcs_code = mapProj->getGcsCode();
293  int datum_code = USER_DEFINED;
294  int ellipsoid_code = USER_DEFINED;
295  const ossimDatum* datum = mapProj->getDatum();
296  if (datum)
297  {
298  datum_code = (int) datum->epsgCode();
299  const ossimEllipsoid* ellipsoid = datum->ellipsoid();
300  if (ellipsoid)
301  ellipsoid_code = ellipsoid->getEpsgCode();
302  }
303 
304  if(mapProj->isGeographic())
305  GTIFKeySet(gtif, GeographicTypeGeoKey, TYPE_SHORT, 1, pcsCode);
306 
307  GTIFKeySet(gtif, GeogGeodeticDatumGeoKey, TYPE_SHORT, 1, datum_code);
308  GTIFKeySet(gtif, ProjectionGeoKey , TYPE_SHORT, 1, pcsCode);
309  GTIFKeySet(gtif, GeogEllipsoidGeoKey, TYPE_SHORT, 1, ellipsoid_code);
310  GTIFKeySet(gtif, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, pcsCode);
311  }
312  else
313  {
314  ossimString datumCode = "WGE";
315  ossimString datumName = "WGE";
316  // Attemp to get the datum code
317  const ossimDatum* datum = mapProj->getDatum();
318  if(datum)
319  {
320  datumCode = datum->code();
321  datumName = datum->name();
322  }
323 
324  short gcs = USER_DEFINED;
325 
326  if (datumCode == "WGE") gcs = GCS_WGS_84;
327  else if (datumCode == "WGD") gcs = GCS_WGS_72;
328  else if (datumCode == "NAR-C") gcs = GCS_NAD83;
329  else if (datumCode == "NAR") gcs = GCS_NAD83;
330  else if (datumCode == "NAS-C") gcs = GCS_NAD27;
331  else if (datumCode == "NAS") gcs = GCS_NAD27;
332  else if (datumCode == "ADI-M") gcs = GCS_Adindan;
333  else if (datumCode == "ARF-M") gcs = GCS_Arc_1950;
334  else if (datumCode == "ARS-M") gcs = GCS_Arc_1960;
335  else if (datumCode == "EUR-7" || datumCode == "EUR-M") gcs = GCS_ED50;
336  else if ((datumCode == "OGB-7") ||
337  (datumCode == "OGB-M") ||
338  (datumCode == "OGB-A") ||
339  (datumCode == "OGB-B") ||
340  (datumCode == "OGB-C") ||
341  (datumCode == "OGB-D")) gcs = GCS_OSGB_1936;
342  else if (datumCode == "TOY-M") gcs = GCS_Tokyo;
343  else if(datumCode == "6055") gcs = 3785; // this is a special google datum so we will give it a gcs code of 3785
344  else
345  {
346  if(traceDebug())
347  {
349  << "DATUM = " << datumCode << " tag not written " << std::endl
350  << "Please let us know so we can add it" << std::endl;
351  }
352  }
353 
354  // ***
355  // ERDAS Imagine < v8.7 has a NAD27 Conus Bug. They are not using the
356  // proper GCS code. They use user-defined fields and Geog citation tag to
357  // define. Sucks! It is an open issue at Leica. This is a work around
358  // flag for this issue.
359  // ***
360  if((datumCode == "NAS-C") && imagineNad27Flag)
361  {
362  gcs = USER_DEFINED;
363 
365  os << "IMAGINE GeoTIFF Support\nCopyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved\n@(#)$RCSfile$ $Revision: 21024 $ $Date: 2012-05-30 04:45:13 -0400 (Wed, 30 May 2012) $\nUnable to match Ellipsoid (Datum) to a GeographicTypeGeoKey value\nEllipsoid = Clarke 1866\nDatum = NAD27 (CONUS)";
366 
367  GTIFKeySet(gtif,
368  GeogCitationGeoKey,
369  TYPE_ASCII,
370  1,
371  os.str().c_str());
372 
373  // User-Defined
374  GTIFKeySet(gtif, GeogGeodeticDatumGeoKey, TYPE_SHORT, 1,
375  KvUserDefined );
376  // User-Defined
377  GTIFKeySet(gtif, GeogEllipsoidGeoKey, TYPE_SHORT, 1,
378  KvUserDefined );
379  }
380  else
381  {
382  GTIFKeySet( gtif, GeographicTypeGeoKey, TYPE_SHORT, 1, gcs );
383  gcsTypeSet = true;
384  }
385 
386  // Write the projection parameters.
387 
388  bool setFalseEastingNorthingFlag = false;
389 
390  if ( projName == "ossimUtmProjection" )
391  {
392  //---
393  // UTM tags needed example from the geo tiff spec page:
394  // ModelTiepointTag = (0, 0, 0, 350807.4, 5316081.3, 0.0)
395  // ModelPixelScaleTag = (100.0, 100.0, 0.0)
396  // GeoKeyDirectoryTag:
397  // GTModelTypeGeoKey = 1 (ModelTypeProjected)
398  // GTRasterTypeGeoKey = 1 (RasterPixelIsArea)
399  // ProjectedCSTypeGeoKey = 32660 (PCS_WGS84_UTM_zone_60N)
400  // PCSCitationGeoKey = "UTM Zone 60 N with WGS84"
401  //
402  // NOTE:
403  // The "ProjectedCSTypeGeoKey" can be constructed using the map zone
404  // and the datum.
405  //---
406  const ossimUtmProjection* utmProjection = dynamic_cast<const ossimUtmProjection*>(mapProj);
407 
408  // Attempt to get the pcs key.
409  int mapZone = utmProjection->getZone();
410  ossimString hemisphere = utmProjection->getHemisphere();
411  short projSysCode=0;
412 
413  //---
414  // Use a projection code that does not imply a datum.
415  // See section "6.3.3.2 Projection Codes" for definition.
416  //
417  // NOTE: The ossim code does not use negative zones to discriminate between
418  // hemispheres. (drb 26 Oct. 2016)
419  //---
420  if ( hemisphere == "N" ) // Northern hemisphere.
421  {
422  projSysCode = 16000 + mapZone;
423  }
424  else // Southern hemisphere.
425  {
426  projSysCode = 16100 + mapZone;
427  }
428 
429  // Set the Projected Coordinate System Type to be user defined.
430  GTIFKeySet(gtif,
431  ProjectedCSTypeGeoKey,
432  TYPE_SHORT,
433  1,
434  USER_DEFINED);
435 
436  if ( !gcsTypeSet )
437  {
438  // Set the geographic type to be user defined.
439  GTIFKeySet(gtif,
440  GeographicTypeGeoKey,
441  TYPE_SHORT,
442  1,
443  USER_DEFINED);
444  }
445 
446  // Set the ProjectionGeoKey in place of the ProjectedCSTypeGeoKey.
447  GTIFKeySet(gtif,
448  ProjectionGeoKey,
449  TYPE_SHORT,
450  1,
451  projSysCode);
452 
454  os << "UTM Zone " << dec << mapZone << hemisphere.c_str()
455  << " with " << datumName << " datum";
456 
457  GTIFKeySet(gtif,
458  PCSCitationGeoKey,
459  TYPE_ASCII,
460  1,
461  os.str().c_str());
462 
463  } // End of "if ( projName == "ossimUtmProjection" )
464 
465  else if(projName == "ossimBngProjection")
466  {
467  // User-Defined
468  GTIFKeySet(gtif, ProjectedCSTypeGeoKey, TYPE_SHORT, 1,
469  PCS_BRITISH_NATIONAL_GRID);//KvUserDefined );
470 
471  // User-Defined
472  GTIFKeySet(gtif, ProjectionGeoKey, TYPE_SHORT, 1,
473  KvUserDefined );
474 
475  GTIFKeySet(gtif,
476  PCSCitationGeoKey,
477  TYPE_ASCII,
478  26,
479  "PCS_British_National_Grid");
480 
481  GTIFKeySet(gtif,
482  ProjCoordTransGeoKey,
483  TYPE_SHORT,
484  1,
485  (uint16)CT_TransverseMercator);
486 
487  GTIFKeySet(gtif,
488  ProjNatOriginLongGeoKey,
489  TYPE_DOUBLE,
490  1,
491  origin.lond());
492 
493  GTIFKeySet(gtif,
494  ProjNatOriginLatGeoKey,
495  TYPE_DOUBLE,
496  1,
497  origin.latd());
498 
499  setFalseEastingNorthingFlag = true;
500 
501  double scale = ossimString(scaleFactor).toDouble();
502 
503  GTIFKeySet(gtif,
504  ProjScaleAtNatOriginGeoKey,
505  TYPE_DOUBLE,
506  1,
507  scale);
508  }
509  else if( projName == "ossimSinusoidalProjection")
510  {
511  GTIFKeySet(gtif,
512  ProjCoordTransGeoKey,
513  TYPE_SHORT,
514  1,
515  (uint16)CT_Sinusoidal);
516 
517  GTIFKeySet(gtif,
518  ProjNatOriginLongGeoKey,
519  TYPE_DOUBLE,
520  1,
521  origin.lond());
522 
523  GTIFKeySet(gtif,
524  ProjNatOriginLatGeoKey,
525  TYPE_DOUBLE,
526  1,
527  origin.latd());
528 
529  setFalseEastingNorthingFlag = true;
530  }
531  else if( (projName == "ossimEquDistCylProjection")||
532  (projName == "ossimLlxyProjection"))
533  {
534  GTIFKeySet(gtif,
535  ProjNatOriginLongGeoKey,
536  TYPE_DOUBLE,
537  1,
538  origin.lond());
539 
540  GTIFKeySet(gtif,
541  ProjNatOriginLatGeoKey,
542  TYPE_DOUBLE,
543  1,
544  origin.latd());
545  }
546  else if ( (projName == "ossimLambertConformalConicProjection") ||
547  (projName == "ossimAlbersProjection") )
548  {
549  //---
550  // Lambert Conformal Conic:
551  // tags needed example from the geo tiff spec page:
552  // ModelTiepointTag = ( 80, 100, 0, 200000, 1500000, 0)
553  // ModelPixelScaleTag = (1000, 1000, 0)
554  // GeoKeyDirectoryTag:
555  // GTModelTypeGeoKey = 1 (ModelTypeProjected)
556  // GTRasterTypeGeoKey = 1 (RasterPixelIsArea)
557  // GeographicTypeGeoKey = 4267 (GCS_NAD27)
558  // ProjectedCSTypeGeoKey = 32767 (user-defined)
559  // ProjectionGeoKey = 32767 (user-defined)
560  // ProjLinearUnitsGeoKey = 9001 (Linear_Meter)
561  // ProjCoordTransGeoKey = 8 (CT_LambertConfConic_2SP)
562  // ProjStdParallel1GeoKey = 41.333
563  // ProjStdParallel2GeoKey = 48.666
564  // ProjCenterLongGeoKey =-120.0
565  // ProjNatOriginLatGeoKey = 45.0
566  // ProjFalseEastingGeoKey, = 200000.0
567  // ProjFalseNorthingGeoKey, = 1500000.0
568  //
569  // NOTE: Albers Same as Lambert with the exception of the
570  // ProjCoordTransGeoKey which is: CT_AlbersEqualArea.
571  //---
572 
573  if (projName == "ossimLambertConformalConicProjection")
574  {
575  GTIFKeySet(gtif,
576  ProjCoordTransGeoKey,
577  TYPE_SHORT,
578  1,
579  (uint16)CT_LambertConfConic_2SP );
580  }
581  else // Albers
582  {
583  GTIFKeySet(gtif,
584  ProjCoordTransGeoKey,
585  TYPE_SHORT,
586  1,
587  (uint16)CT_AlbersEqualArea);
588  }
589 
590  // User-Defined
591  GTIFKeySet(gtif, ProjectedCSTypeGeoKey, TYPE_SHORT, 1,
592  KvUserDefined );
593 
594  // User-Defined
595  GTIFKeySet(gtif, ProjectionGeoKey, TYPE_SHORT, 1,
596  KvUserDefined );
597 
598  double phi1 = ossimString(stdParallel1).toDouble();
599 
600  GTIFKeySet(gtif,
601  ProjStdParallel1GeoKey,
602  TYPE_DOUBLE,
603  1,
604  phi1); // 1st parallel
605 
606  double phi2 = ossimString(stdParallel2).toDouble();
607 
608  GTIFKeySet(gtif,
609  ProjStdParallel2GeoKey,
610  TYPE_DOUBLE,
611  1,
612  phi2); // 2nd parallel
613 
614  GTIFKeySet(gtif,
615  ProjCenterLongGeoKey,
616  TYPE_DOUBLE,
617  1,
618  origin.lond()); // Longitude at the origin.
619 
620  GTIFKeySet(gtif,
621  ProjNatOriginLatGeoKey,
622  TYPE_DOUBLE,
623  1,
624  origin.latd()); // Origin
625 
626  setFalseEastingNorthingFlag = true;
627 
628  } // End of Lambert.
629 
630 
631  else if ( projName == "ossimMercatorProjection" )
632  {
633  GTIFKeySet(gtif,
634  ProjCoordTransGeoKey,
635  TYPE_SHORT,
636  1,
637  (uint16)CT_Mercator);
638 
639  GTIFKeySet(gtif,
640  ProjNatOriginLongGeoKey,
641  TYPE_DOUBLE,
642  1,
643  origin.lond());
644 
645  GTIFKeySet(gtif,
646  ProjNatOriginLatGeoKey,
647  TYPE_DOUBLE,
648  1,
649  origin.latd());
650 
651  setFalseEastingNorthingFlag = true;
652 
653  double scale = ossimString(scaleFactor).toDouble();
654 
655  GTIFKeySet(gtif,
656  ProjScaleAtNatOriginGeoKey,
657  TYPE_DOUBLE,
658  1,
659  scale);
660  }
661  else if ( projName == "ossimTransMercatorProjection" )
662  {
663  //---
664  // Transverse Mercator ( no example in the geo tiff spec.
665  // Requires:
666  // - latitude/longitude of the origin
667  // - easting/northing of some tie point(line/sample 0,0)
668  // - false easting/northing
669  // - The scale factor.
670  //---
671 
672  //---
673  // The easting/northing is the distance from the origin plus the
674  // false easting/northing. In other words if line 0 is 5,000
675  // meters from the origin and the false northing is 5,000 meters,
676  // then the northing would be 10,000. The same goes for the easting.
677  //---
678  GTIFKeySet(gtif,
679  ProjCoordTransGeoKey,
680  TYPE_SHORT,
681  1,
682  (uint16)CT_TransverseMercator);
683 
684  // User-Defined
685  GTIFKeySet(gtif, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined );
686 
687  // User-Defined
688  GTIFKeySet(gtif, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined );
689 
690  GTIFKeySet(gtif,
691  ProjNatOriginLongGeoKey,
692  TYPE_DOUBLE,
693  1,
694  origin.lond());
695 
696  GTIFKeySet(gtif,
697  ProjNatOriginLatGeoKey,
698  TYPE_DOUBLE,
699  1,
700  origin.latd());
701 
702  setFalseEastingNorthingFlag = true;
703 
704  double scale = ossimString(scaleFactor).toDouble();
705 
706  GTIFKeySet(gtif,
707  ProjScaleAtNatOriginGeoKey,
708  TYPE_DOUBLE,
709  1,
710  scale);
711  } // End of TM
712 
713  if (setFalseEastingNorthingFlag == true)
714  {
715 
716  GTIFKeySet(gtif,
717  ProjFalseEastingGeoKey,
718  TYPE_DOUBLE,
719  1,
720  falseEasting);
721 
722  GTIFKeySet(gtif,
723  ProjFalseNorthingGeoKey,
724  TYPE_DOUBLE,
725  1,
726  falseNorthing);
727  }
728  }
729 
730  //---
731  // Set the model type and units.
732  //---
733  if (units == ossimGeoTiff::ANGULAR_DEGREE)
734  {
735  GTIFKeySet(gtif,
736  GTModelTypeGeoKey,
737  TYPE_SHORT,
738  1,
739  ModelTypeGeographic);
740 
741  // Set the units key.
742  GTIFKeySet(gtif,
743  GeogAngularUnitsGeoKey,
744  TYPE_SHORT,
745  1,
746  units);
747  }
748  else
749  {
750  GTIFKeySet(gtif,
751  GTModelTypeGeoKey,
752  TYPE_SHORT,
753  1,
754  ModelTypeProjected);
755 
756  // Set the units key.
757  GTIFKeySet(gtif,
758  ProjLinearUnitsGeoKey,
759  TYPE_SHORT,
760  1,
761  units);
762  }
763 
764  // Set the ellipsoid major/minor axis.
765  GTIFKeySet(gtif,
766  GeogSemiMajorAxisGeoKey,
767  TYPE_DOUBLE,
768  1,
769  mapProj->getA());
770 
771  GTIFKeySet(gtif,
772  GeogSemiMinorAxisGeoKey,
773  TYPE_DOUBLE,
774  1,
775  mapProj->getB());
776 
777  // Set the pixel type.
778  if (projectionInfo->getPixelType() == OSSIM_PIXEL_IS_POINT)
779  {
780  // Tie point relative to center of pixel.
781  GTIFKeySet(gtif, GTRasterTypeGeoKey, TYPE_SHORT, 1, RasterPixelIsPoint);
782  }
783  else
784  {
785  // Tie point relative to upper left corner of pixel
786  GTIFKeySet(gtif, GTRasterTypeGeoKey, TYPE_SHORT, 1, RasterPixelIsArea);
787  }
788 
789  //---
790  // Set the tie point and scale.
791  //---
792  double tiePoints[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
793  double pixScale[3] = { 0.0, 0.0, 0.0 };
794  switch (units)
795  {
796  case LINEAR_FOOT:
797  {
798  tiePoints[3] = ossim::mtrs2ft(projectionInfo->ulEastingNorthingPt().x);
799  tiePoints[4] = ossim::mtrs2ft(projectionInfo->ulEastingNorthingPt().y);
800  pixScale[0] = ossim::mtrs2ft(projectionInfo->getMetersPerPixel().x);
801  pixScale[1] = ossim::mtrs2ft(projectionInfo->getMetersPerPixel().y);
802  falseEasting = ossim::mtrs2ft(falseEasting);
803  falseNorthing = ossim::mtrs2ft(falseNorthing);
804 
805  break;
806  }
808  {
809  tiePoints[3] = ossim::mtrs2usft(projectionInfo->ulEastingNorthingPt().x);
810  tiePoints[4] = ossim::mtrs2usft(projectionInfo->ulEastingNorthingPt().y);
811  pixScale[0] = ossim::mtrs2usft(projectionInfo->getMetersPerPixel().x);
812  pixScale[1] = ossim::mtrs2usft(projectionInfo->getMetersPerPixel().y);
813  falseEasting = ossim::mtrs2usft(falseEasting);
814  falseNorthing = ossim::mtrs2usft(falseNorthing);
815  break;
816  }
817  case ANGULAR_DEGREE:
818  {
819  tiePoints[3] = projectionInfo->ulGroundPt().lond();
820  tiePoints[4] = projectionInfo->ulGroundPt().latd();
821  pixScale[0] = projectionInfo->getDecimalDegreesPerPixel().x;
822  pixScale[1] = projectionInfo->getDecimalDegreesPerPixel().y;
823  break;
824  }
825  case LINEAR_METER:
826  default:
827  {
828  tiePoints[3] = projectionInfo->ulEastingNorthingPt().x;
829  tiePoints[4] = projectionInfo->ulEastingNorthingPt().y;
830  pixScale[0] = projectionInfo->getMetersPerPixel().x;
831  pixScale[1] = projectionInfo->getMetersPerPixel().y;
832  break;
833  }
834 
835  } // End of "switch (units)"
836 
837  TIFFSetField( tifPtr, TIFFTAG_GEOTIEPOINTS, 6, tiePoints );
838  TIFFSetField( tifPtr, TIFFTAG_GEOPIXELSCALE, 3, pixScale );
839 
840 
841  GTIFWriteKeys(gtif); // Write out geotiff tags.
842  GTIFFree(gtif);
843 
844  return true;
845 }
846 
848  const ossimIrect& rect,
849  const ossimProjection* proj,
850  std::vector<ossim_uint8>& buf,
851  ossimPixelType pixelType )
852 {
853  //---
854  // Snip from The "GeoTIFF Box" Specification for JPEG 2000 Metadata:
855  // This box contains a valid GeoTIFF image. The image is "degenerate", in
856  // that it represents a very simple image with specific constraints:
857  // . the image height and width are both 1
858  // . the datatype is 8-bit
859  // . the colorspace is grayscale
860  // . the (single) pixel must have a value of 0 for its (single) sample
861  //
862  // NOTE: It also states little endian but I think libtiff writes whatever
863  // endianesss the host is.
864  //
865  // Also assuming class tiff for now. Not big tiff.
866  //---
867  bool result = true;
868 
869  TIFF* tiff = XTIFFOpen(tmpFile.c_str(), "w");
870  if (tiff)
871  {
872  // Write the projection info out.
874  if(mapProj)
875  {
877  = new ossimMapProjectionInfo(mapProj, ossimDrect(rect));
878 
879  // Set the pixel type to point of area.
880  projectionInfo->setPixelType(pixelType);
881 
882  // Write the geotiff keys.
883  ossimGeoTiff::writeTags(tiff, projectionInfo, false);
884  }
885 
886  // Basic tiff tags.
887  TIFFSetField( tiff, TIFFTAG_IMAGEWIDTH, 1 );
888  TIFFSetField( tiff, TIFFTAG_IMAGELENGTH, 1 );
889  TIFFSetField( tiff, TIFFTAG_BITSPERSAMPLE, 8 );
890  TIFFSetField( tiff, TIFFTAG_SAMPLESPERPIXEL, 1 );
891  TIFFSetField( tiff, TIFFTAG_ROWSPERSTRIP, 1 );
892  TIFFSetField( tiff, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG );
893  TIFFSetField( tiff, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK );
894 
895  // One pixel image:
896  ossim_uint8 pixel = 0;
897  TIFFWriteEncodedStrip( tiff, 0, (char *) &pixel, 1 );
898 
899  TIFFWriteDirectory( tiff );
900  XTIFFClose( tiff );
901 
902  // Get the size. Note 16 bytes added for the JP2 UUID:
903  const std::vector<ossim_uint8>::size_type UUID_SIZE = 16;
904  const std::vector<ossim_uint8>::size_type BOX_SIZE = UUID_SIZE +
905  static_cast<std::vector<ossim_uint8>::size_type>(tmpFile.fileSize());
906 
907  // Create the buffer.
908  buf.resize( BOX_SIZE );
909 
910  if ( BOX_SIZE == buf.size() )
911  {
912  const ossim_uint8 GEOTIFF_UUID[UUID_SIZE] =
913  {
914  0xb1, 0x4b, 0xf8, 0xbd,
915  0x08, 0x3d, 0x4b, 0x43,
916  0xa5, 0xae, 0x8c, 0xd7,
917  0xd5, 0xa6, 0xce, 0x03
918  };
919 
920  // Copy the UUID.
921  std::vector<ossim_uint8>::size_type i;
922  for (i = 0; i < UUID_SIZE; ++i)
923  {
924  buf[i] = GEOTIFF_UUID[i];
925  }
926 
927  // Copy the tiff.
928  std::ifstream str;
929  str.open(tmpFile.c_str(), ios::in | ios::binary);
930  if (str.is_open())
931  {
932  char ch;
933  for (; i < BOX_SIZE; ++i)
934  {
935  str.get(ch);
936  buf[i] = static_cast<ossim_uint8>(ch);
937  }
938  }
939  }
940  else
941  {
942  result = false;
943  }
944 
945  // Remove the temp file.
946  tmpFile.remove();
947 
948  }
949  else
950  {
951  result = false;
952 
953  if (traceDebug())
954  {
956  << "ossimGeoTiff::writeJp2GeotiffBox ERROR:\n"
957  << "Could not open " << tmpFile << std::endl;
958  }
959  }
960  return result;
961 }
962 
964 {
965  bool result = false;
966 
967  TIFF* tiff = XTIFFOpen(file.c_str(), "r");
968  if(tiff)
969  {
970  result = readTags(tiff, entryIdx, true);
971  }
972 
973  return result;
974 }
975 
976 bool ossimGeoTiff::readTags(std::shared_ptr<ossim::TiffHandlerState> state, ossim_uint32 entryIdx)
977 {
978  std::lock_guard<std::mutex> lock(theMutex);
979 
980  if(!state) return false;
981 
982  if(!state->checkBool(entryIdx, "is_geotiff"))
983  {
984  return false;
985  }
986  ossimString value;
987  theLength = state->getImageLength(entryIdx);
988  theWidth = state->getImageWidth(entryIdx);
989  theScaleFactor = 0.0;
990  theModelType = state->getModelType(entryIdx);
991  theGcsCode = state->getGcsCode(entryIdx);
992  theDatumCode = state->getDatumCode(entryIdx);
993  theAngularUnits = state->getAngularUnits(entryIdx);
994  theLinearUnitsCode = state->getLinearUnits(entryIdx);
995  theGeoKeysPresentFlag = true;
996 
998  {
999  //---
1000  // Hack for bug, where the libgeotiff funtion GTIFGetDefn sets the angular units
1001  // incorrectly to ANGULAR_DMS_HEMISPHERE:
1002  if ( traceDebug() )
1003  {
1005  << " WARNING ossimGeoTiff::addImageGeometry:"
1006  << "The angular units (key 2054) is set to ANGULAR_DMS_HEMISPHERE!"
1007  << "\nAssuming \"Angular_Degree\"..." << std::endl;
1008  }
1010  }
1011 
1012  if (state->getValue(value, entryIdx, "tifftag.pcs_code"))
1013  {
1014  thePcsCode = value.toInt32();
1015  parsePcsCode();
1016  }
1017  //---
1018  // ESH 2/2008 -- Handle geotiff's with state plane coordinate systems produced by ERDAS.
1019  // They use the citation filed to specify the geometry (complete HACK by Erdas)
1020  //---
1021  else
1022  {
1023  ossimString gTCitation;
1024  if(state->getValue(gTCitation, entryIdx, "tifftag.citation"))
1025  {
1026  // Extract state plane string from the citation key
1027  ossimString projStrTemp =
1028  gTCitation.afterRegExp( "Projection Name = " );
1029 
1030  ossimString projStr = projStrTemp.beforeRegExp( "\n" );
1031  if ( projStr.empty() == false )
1032  {
1034  ossimProjection* proj = f->createProjection(projStr);
1035  ossimMapProjection* map_proj = PTR_CAST(ossimMapProjection, proj);
1036  parseProjection(map_proj);
1037  }
1038  } // End of "if(GTIFKeyGet(gtif, GTCitationGeoKey..."
1039  }
1040 
1043  state->getValue(thePcsCitation, entryIdx, "tifftag.pcs_citation");
1044  if(state->getValue(value, entryIdx, "tifftag.coord_trans_code"))
1045  {
1046  theCoorTransGeoCode = value.toInt32();
1047  }
1048  if(state->getValue(value, entryIdx, "tifftag.std_parallel_1"))
1049  {
1050  theStdPar1 = value.toDouble();
1051  }
1052  if(state->getValue(value, entryIdx, "tifftag.std_parallel_2"))
1053  {
1054  theStdPar2 = value.toDouble();
1055  }
1056  if(state->getValue(value, entryIdx, "tifftag.origin_lon"))
1057  {
1058  theOriginLon = value.toDouble();
1059  }
1060  else if(state->getValue(value, entryIdx, "tifftag.center_lon"))
1061  {
1062  theOriginLon = value.toDouble();
1063  }
1064  else if(state->getValue(value, entryIdx, "tifftag.false_origin_lon"))
1065  {
1066  theOriginLon = value.toDouble();
1067  }
1068  else if(state->getValue(value, entryIdx, "tifftag.straight_vert_pole_lon"))
1069  {
1070  theOriginLon = value.toDouble();
1071  }
1072 
1073  if(state->getValue(value, entryIdx, "tifftag.origin_lat"))
1074  {
1075  theOriginLat = value.toDouble();
1076  }
1077  else if(state->getValue(value, entryIdx, "tifftag.center_lat"))
1078  {
1079  theOriginLat = value.toDouble();
1080  }
1081  else if(state->getValue(value, entryIdx, "tifftag.false_origin_lat"))
1082  {
1083  theOriginLat = value.toDouble();
1084  }
1085 
1086  if(state->getValue(value, entryIdx, "tifftag.false_easting"))
1087  {
1088  theFalseEasting = value.toDouble();
1089  }
1090  if(state->getValue(value, entryIdx, "tifftag.false_northing"))
1091  {
1092  theFalseNorthing = value.toDouble();
1093  }
1094  if(state->getValue(value, entryIdx, "tifftag.scale_factor"))
1095  {
1096  theScaleFactor = value.toDouble();
1097  }
1098  theScale.clear();
1099  if(state->getGeoPixelScale(theScale, entryIdx))
1100  {
1101  if ( theModelType == ModelTypeGeographic )
1102  {
1103  // The origin latitude must be computed so as to achieve the proper horizontal scaling:
1105  }
1106  }
1107  theTiePoint.clear();
1108  state->getGeoTiePoints(theTiePoint, entryIdx);
1109  theModelTransformation.clear();
1110  state->getGeoTransMatrix(theModelTransformation, entryIdx);
1111  theDoubleParam.clear();
1112  state->getGeoDoubleParams(theDoubleParam, entryIdx);
1113 
1114  state->getValue(theAsciiParam, entryIdx, "tifftag.geo_ascii_params");
1115 
1116  setOssimProjectionName(state, entryIdx);
1117  setOssimDatumName(state, entryIdx);
1118 
1119  return true;
1120 }
1121 
1123  TIFF* tiff, ossim_uint32 entryIdx, bool ownTiffPtrFlag)
1124 {
1125  std::lock_guard<std::mutex> lock(theMutex);
1126 
1127  if ( !tiff )
1128  {
1129  return false;
1130  }
1131 
1132  if(!TIFFSetDirectory(tiff, (ossim_uint16)entryIdx))
1133  {
1134  return false;
1135  }
1136 
1137  GTIF* gtif = GTIFNew(tiff);
1138  if(!gtif)
1139  {
1140  return false;
1141  }
1142 
1143  if(theTiffPtr)
1144  {
1145  XTIFFClose(theTiffPtr);
1146  }
1147  theTiffPtr = tiff;
1148 
1150  {
1151  delete thePrivateDefinitions->m_defs;;
1152  }
1153  thePrivateDefinitions->m_defs = new GTIFDefn;
1154  GTIFGetDefn(gtif, thePrivateDefinitions->m_defs);
1155  ossim_uint32 idx = 0;
1156  theGeoKeysPresentFlag = true;
1157  if(traceDebug())
1158  {
1160  << "ossimGeoTiff::readTags: Raw Geotiff Tags are\n";
1161  GTIFPrint(gtif,0,0);
1162  }
1163  if(TIFFGetField( theTiffPtr,
1164  TIFFTAG_IMAGELENGTH,
1165  &theLength ))
1166  {
1167  }
1168  if(TIFFGetField( theTiffPtr,
1169  TIFFTAG_IMAGEWIDTH,
1170  &theWidth ))
1171  {
1172  }
1173  theScaleFactor = 0.0;
1180 
1182  {
1183  //---
1184  // Hack for bug, where the libgeotiff funtion GTIFGetDefn sets the angular units
1185  // incorrectly to ANGULAR_DMS_HEMISPHERE:
1186  if ( traceDebug() )
1187  {
1189  << " WARNING ossimGeoTiff::addImageGeometry:"
1190  << "The angular units (key 2054) is set to ANGULAR_DMS_HEMISPHERE!"
1191  << "\nAssuming \"Angular_Degree\"..." << std::endl;
1192  }
1194  }
1195 
1196 #if 0
1197  ossim_uint16 modelType;
1198  if(GTIFKeyGet(gtif, GTModelTypeGeoKey, &modelType, 0, 1))
1199  {
1200  theModelType = (ModelType)modelType;
1201  }
1202 #endif
1203  if(GTIFKeyGet(gtif, GTRasterTypeGeoKey, &theRasterType, 0, 1))
1204  {
1205  }
1206 
1207 #if 0
1208  if(GTIFKeyGet(gtif, GeographicTypeGeoKey, &theGcsCode, 0, 1))
1209  {
1210  }
1211  if(GTIFKeyGet(gtif, GeogGeodeticDatumGeoKey, &theDatumCode, 0, 1))
1212  {
1213  }
1214  if(GTIFKeyGet(gtif, GeogAngularUnitsGeoKey, &theAngularUnits, 0, 1))
1215  {
1216  }
1217 #endif
1218 
1219  if (GTIFKeyGet(gtif, ProjectedCSTypeGeoKey, &thePcsCode, 0, 1))
1220  parsePcsCode();
1221 
1222  //---
1223  // ESH 2/2008 -- Handle geotiff's with state plane coordinate systems produced by ERDAS.
1224  // They use the citation filed to specify the geometry (complete HACK by Erdas)
1225  //---
1226  else
1227  {
1228  const int CITATION_STRING_SIZE = 512;
1229  char citationStr[CITATION_STRING_SIZE];
1230  if ( GTIFKeyGet(gtif, GTCitationGeoKey, &citationStr,
1231  0, CITATION_STRING_SIZE))
1232  {
1233  ossimString gTCitation = citationStr; // key 1026
1234 
1235  // Extract state plane string from the citation key
1236  ossimString projStrTemp =
1237  gTCitation.afterRegExp( "Projection Name = " );
1238 
1239  ossimString projStr = projStrTemp.beforeRegExp( "\n" );
1240  if ( projStr.empty() == false )
1241  {
1243  ossimProjection* proj = f->createProjection(projStr);
1244  ossimMapProjection* map_proj = PTR_CAST(ossimMapProjection, proj);
1245  parseProjection(map_proj);
1246  }
1247  } // End of "if(GTIFKeyGet(gtif, GTCitationGeoKey..."
1248  }
1249 
1250  char* buf = 0;
1253  if(GTIFKeyGet(gtif, PCSCitationGeoKey , &buf, 0, 1))
1254  {
1255  thePcsCitation = ossimString(buf);
1256  }
1257  GTIFKeyGet(gtif, ProjCoordTransGeoKey , &theCoorTransGeoCode, 0, 1);
1258  for(idx = 0; idx < (ossim_uint32)thePrivateDefinitions->m_defs->nParms; ++idx)
1259  {
1260  switch(thePrivateDefinitions->m_defs->ProjParmId[idx])
1261  {
1262  case ProjStdParallel1GeoKey:
1263  {
1264  theStdPar1 = thePrivateDefinitions->m_defs->ProjParm[idx];
1265  break;
1266  }
1267  case ProjStdParallel2GeoKey:
1268  {
1269  theStdPar2 = thePrivateDefinitions->m_defs->ProjParm[idx];
1270  break;
1271  }
1272  case ProjOriginLongGeoKey:
1273  {
1274  theOriginLon = thePrivateDefinitions->m_defs->ProjParm[idx];
1275  break;
1276  }
1277  case ProjOriginLatGeoKey:
1278  {
1279  theOriginLat = thePrivateDefinitions->m_defs->ProjParm[idx];
1280  break;
1281  }
1282  case ProjFalseEastingGeoKey:
1283  {
1284  theFalseEasting = thePrivateDefinitions->m_defs->ProjParm[idx];
1285  break;
1286  }
1287  case ProjFalseNorthingGeoKey:
1288  {
1289  theFalseNorthing = thePrivateDefinitions->m_defs->ProjParm[idx];
1290  break;
1291  }
1292  case ProjCenterLongGeoKey:
1293  {
1294  theOriginLon = thePrivateDefinitions->m_defs->ProjParm[idx];
1295  break;
1296  }
1297  case ProjCenterLatGeoKey:
1298  {
1299  theOriginLat = thePrivateDefinitions->m_defs->ProjParm[idx];
1300  break;
1301  }
1302  case ProjFalseOriginLatGeoKey:
1303  {
1304  theOriginLat = thePrivateDefinitions->m_defs->ProjParm[idx];
1305  break;
1306  }
1307  case ProjFalseOriginLongGeoKey:
1308  case ProjStraightVertPoleLongGeoKey:
1309  {
1310  theOriginLon = thePrivateDefinitions->m_defs->ProjParm[idx];
1311  break;
1312  }
1313  case ProjScaleAtNatOriginGeoKey:
1314  {
1315  theScaleFactor = thePrivateDefinitions->m_defs->ProjParm[idx];
1316  break;
1317  }
1318  }
1319  }
1320 
1321 #if 0
1322  if(GTIFKeyGet(gtif, ProjStdParallel1GeoKey, &theStdPar1, 0, 1))
1323  {
1324  }
1325  if(GTIFKeyGet(gtif, ProjStdParallel2GeoKey, &theStdPar2, 0, 1))
1326  {
1327  }
1328  if(GTIFKeyGet(gtif, ProjNatOriginLongGeoKey, &tempDouble, 0, 1))
1329  {
1330  theOriginLon = tempDouble;
1331  }
1332  else if(GTIFKeyGet(gtif, ProjOriginLongGeoKey, &tempDouble, 0, 1))
1333  {
1334  theOriginLon = tempDouble;
1335  }
1336  if(GTIFKeyGet(gtif, ProjNatOriginLatGeoKey, &tempDouble, 0, 1))
1337  {
1338  theOriginLat = tempDouble;
1339  }
1340  else if(GTIFKeyGet(gtif, ProjOriginLatGeoKey, &tempDouble, 0, 1))
1341  {
1342  theOriginLat = tempDouble;
1343  }
1344  if(GTIFKeyGet(gtif, ProjFalseEastingGeoKey, &theFalseEasting, 0, 1))
1345  {
1346  }
1347  if(GTIFKeyGet(gtif, ProjFalseNorthingGeoKey, &theFalseNorthing, 0, 1))
1348  {
1349  }
1350  if(GTIFKeyGet(gtif, ProjCenterLongGeoKey, &theCenterLon, 0, 1))
1351  {
1352  }
1353  if(GTIFKeyGet(gtif, ProjCenterLatGeoKey, &theCenterLat, 0, 1))
1354  {
1355  }
1356  if(GTIFKeyGet(gtif, ProjScaleAtNatOriginGeoKey, &theScaleFactor, 0, 1))
1357  {
1358  }
1359 #endif
1360  theScale.clear();
1361  ossim_uint16 pixScaleSize = 0;
1362  double* pixScale=0;
1363  if(TIFFGetField(theTiffPtr, TIFFTAG_GEOPIXELSCALE, &pixScaleSize, &pixScale))
1364  {
1365  theScale.insert(theScale.begin(), pixScale, pixScale+pixScaleSize);
1366  if ( theModelType == ModelTypeGeographic )
1367  {
1368  // The origin latitude must be computed so as to achieve the proper horizontal scaling:
1370  }
1371  }
1372  theTiePoint.clear();
1373  ossim_uint16 tiePointSize = 0;
1374  double* tiepoints=0;
1375  if(TIFFGetField(theTiffPtr, TIFFTAG_GEOTIEPOINTS, &tiePointSize, &tiepoints))
1376  {
1377  theTiePoint.insert(theTiePoint.begin(), tiepoints, tiepoints+tiePointSize);
1378 
1379  // ESH 05/2009 -- If the image is in a projected coordinate system, the
1380  // tiepoints will be projected coordinates not lat/lon. Let's avoid setting
1381  // the origin lon/lat to projected x/y. Fix for ticket #711.
1382  //if ( theModelType == ModelTypeGeographic )
1383  //{
1384  // if(ossim::isnan(theOriginLon) &&
1385  // (pixScaleSize > 1) &&
1386  // (tiePointSize > 3))
1387  // {
1388  // theOriginLon = tiepoints[3] - tiepoints[0] * pixScale[0];
1389  // }
1390  //
1391  // if(ossim::isnan(theOriginLat) && (pixScaleSize > 1) && (tiePointSize > 3))
1392  // {
1393  // theOriginLat = tiepoints[4] + tiepoints[1] * fabs(pixScale[1]);
1394  // }
1395  //}
1396  }
1397  theModelTransformation.clear();
1398  ossim_uint16 transSize = 0;
1399  double* trans = 0;
1400 
1401  if(TIFFGetField(theTiffPtr, TIFFTAG_GEOTRANSMATRIX, &transSize, &trans))
1402  {
1404  trans, trans+transSize);
1405  }
1406 // if(!theTiePoint.size()&&(theModelTransform.size()==16))
1407 // {
1408 // // for now we will stuff the tie point with the model transform tie points.
1409 // theTiePoint.push_back(0.0);
1410 // theTiePoint.push_back(0.0);
1411 // theTiePoint.push_back(0.0);
1412 // theTiePoint.push_back(theModelTransformation[3]);
1413 // theTiePoint.push_back(theModelTransformation[7]);
1414 // theTiePoint.push_back(0.0);
1415 // }
1416  ossim_uint16 doubleParamSize = 0;
1417  double* tempDoubleParam = 0;
1418  theDoubleParam.clear();
1419  if(TIFFGetField(theTiffPtr, TIFFTAG_GEODOUBLEPARAMS, &doubleParamSize, &tempDoubleParam))
1420  {
1421  theDoubleParam.insert(theDoubleParam.begin(),
1422  tempDoubleParam,
1423  tempDoubleParam+doubleParamSize);
1424  }
1425 
1426  char* tempAsciiParam=0;
1427  theAsciiParam = "";
1428 
1429  // Note: this tag does not have the setting set to return the size
1430  // so this call is only a 3 argument call without a size parameter
1431  if(TIFFGetField(theTiffPtr, TIFFTAG_GEOASCIIPARAMS, &tempAsciiParam))
1432  {
1433  theAsciiParam = tempAsciiParam;
1434  }
1435 
1437  {
1438  setOssimProjectionName(); // Initialize the ossim projection name.
1439  setOssimDatumName(); // Initialize the ossim datum name.
1440  }
1441 
1442  // commenting this out. Frank mentioned the GTIFFGetDefn which in geo_normalize
1443  // this should be all we need.
1444  //
1445 #if 0
1446  /*
1447  ESH 05/2009: Replacing badly broken code for making
1448  use of TIFFTAG_GEODOUBLEPARAMS.
1449 
1450  Read the geokey directory tag in order to see how
1451  the TIFFTAG_GEODOUBLEPARAMS are defined.
1452 
1453  For structure of geokey directory, see:
1454  http://www.remotesensing.org/geotiff/spec/geotiff2.4.html
1455  */
1456  ossim_uint16 gkdParamSize = 0;
1457  ossim_uint16* gkdParams = 0;
1458  if(TIFFGetField(theTiffPtr, TIFFTAG_GEOKEYDIRECTORY, &gkdParamSize, &gkdParams))
1459  {
1460  ossim_uint16 numKeys = gkdParams ? gkdParams[3] : 0;
1461  ossim_uint16 key = 0;
1462  for( key=0; key<numKeys; ++key )
1463  {
1464  ossim_uint16 loc = (key+1) * 4;
1465  ossim_uint16 ind = gkdParams[loc+3];
1466 
1467  if ( gkdParams[loc+1] == TIFFTAG_GEODOUBLEPARAMS &&
1468  gkdParams[loc+2] == 1 &&
1469  ind >= 0 && ind < doubleParamSize )
1470  {
1471  double dval = theDoubleParam[ind];
1472  switch( gkdParams[loc] )
1473  {
1474  case ProjStdParallel1GeoKey: theStdPar1 = dval; break;
1475  case ProjStdParallel2GeoKey: theStdPar2 = dval; break;
1476  case ProjNatOriginLongGeoKey: theOriginLon = dval; break;
1477  /* case ProjOriginLongGeoKey: theOriginLon = dval; break; (alias) */
1478  case ProjNatOriginLatGeoKey: theOriginLat = dval; break;
1479  /* case ProjOriginLatGeoKey: theOriginLat = dval; break; (alias) */
1480  case ProjFalseEastingGeoKey: theFalseEasting = dval; break;
1481  case ProjFalseNorthingGeoKey: theFalseNorthing = dval; break;
1482  case ProjCenterLongGeoKey: theCenterLon = dval; break;
1483  case ProjCenterLatGeoKey: theCenterLat = dval; break;
1484  case ProjScaleAtNatOriginGeoKey: theScaleFactor = dval; break;
1485  default:
1486  if(traceDebug())
1487  {
1489  << "ossimGeoTiff::readTags: Unrecognized geokey directory entry."
1490  << "\ngeokey directory index = " << loc
1491  << "\ngeokey = " << gkdParams[loc]
1492  << "\ndouble array index = " << ind
1493  << "\ndval = " << dval
1494  << std::endl;
1495  }
1496  break;
1497  }
1498  }
1499  }
1500  }
1501 #endif
1502 
1503  GTIFFree(gtif);
1504 
1505  if (ownTiffPtrFlag == false)
1506  {
1507  //---
1508  // Zero out the pointer so the destructor doesn't close it on some
1509  // external code.
1510  //---
1511  theTiffPtr = 0;
1512  }
1513 
1514  return true;
1515 }
1516 
1517 bool ossimGeoTiff::addImageGeometry(ossimKeywordlist& kwl, const char* prefix) const
1518 {
1519  if(traceDebug())
1520  {
1522  << "ossimGeoTiff::addImageGeometry: Entered............."
1523  << std::endl;
1524  }
1525 
1526  // NOT SURE THIS IS A GOOD IDEA HERE. KEPT FOR LEGACY SAKE. (OLK 5/10)
1527  if(theGcsCode == 3785)
1528  {
1531  if(proj.valid())
1532  {
1533  proj->saveState(kwl, prefix);
1534  }
1535  }
1536 
1537  //---
1538  // Sanity check...
1539  // NOTE: It takes six doubles to make one tie point ie:
1540  // x,y,z,longitude,latitude,height or x,y,z,easting,northing,height
1541  //---
1542  if (theErrorStatus ||
1543  (!usingModelTransform() &&
1544  ( (theScale.size() < 2) && // no scale
1545  ( theTiePoint.size() < 24) ) ) )//need at least 3 ties if no scale.
1546  {
1547  if(traceDebug())
1548  {
1550  << "ossimGeoTiff::addImageGeometry: Failed sanity check "
1551  << std::endl;
1552  if(theErrorStatus)
1553  {
1555  << "for error status" << std::endl;
1556  }
1557  else if( theTiePoint.size()<5)
1558  {
1560  << "for tie points, size = " << theTiePoint.size()
1561  << std::endl;
1562  }
1563  else
1564  {
1565  ossimNotify(ossimNotifyLevel_DEBUG) << "for scale" << std::endl;
1566  }
1567  }
1568  return false;
1569  }
1570 
1571  double x_tie_point = 0.0;
1572  double y_tie_point = 0.0;
1573  ossim_uint32 tieCount = (ossim_uint32)theTiePoint.size()/6;
1574 
1575  if( (theScale.size() == 3) && (tieCount == 1))
1576  {
1577  //---
1578  // Shift the tie point to the (0, 0) position if it's not already.
1579  //
1580  // Note:
1581  // Some geotiff writers like ERDAS IMAGINE set the "GTRasterTypeGeoKey"
1582  // key to RasterPixelIsArea, then set the tie point to (0.5, 0.5).
1583  // This really means "RasterPixelIsPoint" with a tie point of (0.0, 0.0).
1584  // Anyway we will check for this blunder and attempt to do the right
1585  // thing...
1586  //---
1587  x_tie_point = theTiePoint[3] - theTiePoint[0]*theScale[0];
1588  y_tie_point = theTiePoint[4] + theTiePoint[1]*theScale[1];
1589  }
1590  else if(tieCount > 1)
1591  {
1592  //---
1593  // Should we check the model type??? (drb)
1594  // if (theModelType == ModelTypeGeographic)
1595  //---
1596  if(tieCount >= 4)
1597  {
1598  ossimTieGptSet tieSet;
1599  getTieSet(tieSet);
1600 
1601  if(tieCount > 4)
1602  {
1603  // create a cubic polynomial model
1604  //ossimRefPtr<ossimPolynomProjection> proj = new ossimPolynomProjection;
1605  //proj->setupOptimizer("1 x y x2 xy y2 x3 y3 xy2 x2y z xz yz");
1607  proj->optimizeFit(tieSet);
1608  proj->saveState(kwl, prefix);
1609  if(traceDebug())
1610  {
1612  << "ossimGeoTiff::addImageGeometry: "
1613  << "Creating a Cubic polynomial projection" << std::endl;
1614  }
1615 
1616  return true;
1617 
1618  }
1619  else if(tieCount == 4)
1620  {
1621  // create a bilinear model
1622 
1623  // Should we check the model type (drb)
1624  // if (theModelType == ModelTypeGeographic)
1625 
1627  proj->optimizeFit(tieSet);
1628  proj->saveState(kwl, prefix);
1629 
1630  if(traceDebug())
1631  {
1633  << "ossimGeoTiff::addImageGeometry: "
1634  << "Creating a bilinear projection" << std::endl;
1635  }
1636  return true;
1637  }
1638  else
1639  {
1641  << "ossimGeoTiff::addImageGeometry: "
1642  << "Not enough tie points to create a interpolation model"
1643  <<std::endl;
1644  }
1645  return false;
1646  }
1647  }
1648  else if (usingModelTransform())
1649  {
1650  if(traceDebug())
1651  {
1653  << "ossimGeoTiff::addImageGeometry: Do not support rotated "
1654  << "map models yet. You should provide the image as a sample "
1655  << "and we will fix it" << std::endl;
1656  }
1657  }
1658 
1659  if ((theRasterType == PIXEL_IS_AREA))
1660  {
1661  // Since the internal pixel representation is "point", shift the
1662  // tie point to be relative to the center of the pixel.
1663  if (theScale.size() > 1)
1664  {
1665  x_tie_point += (theScale[0])/2.0;
1666  y_tie_point -= (theScale[1])/2.0;
1667  }
1668  }
1669 
1670  if( thePcsCode && (thePcsCode != USER_DEFINED) )
1671  {
1672  ossimString epsg_spec ("EPSG:");
1673  epsg_spec += ossimString::toString(thePcsCode);
1676  if (proj.valid())
1677  {
1678  proj->saveState(kwl, prefix);
1679  }
1680  // Should be some else "WARNING" here maybe. (drb)
1681  }
1682  else if (getOssimProjectionName() == "unknown")
1683  {
1684  //---
1685  // Get the projection type. If unknown no point going on, so get out.
1686  //---
1687  return false;
1688  }
1689  else
1690  {
1691  // No PCS code but we do have a projection name
1692  // Add these for all projections.
1695  }
1696 
1697  // Now set the image-specific projection info (scale and image tiepoint):
1699  {
1701  {
1703  << "WARNING ossimGeoTiff::addImageGeometry:"
1704  << "\nNot coded yet for unit type: "
1705  << theAngularUnits << endl;
1706  return false;
1707  }
1708 
1709  //---
1710  // Tiepoint
1711  // Have data with tie points -180.001389 so use ossimGpt::wrap() to handle:
1712  //---
1713  ossimGpt tieGpt(y_tie_point, x_tie_point, 0.0);
1714  tieGpt.wrap();
1715  ossimDpt tiepoint(tieGpt);
1716  kwl.add(prefix, ossimKeywordNames::TIE_POINT_XY_KW, tiepoint.toString(), true);
1717  kwl.add(prefix, ossimKeywordNames::TIE_POINT_UNITS_KW, "degrees", true);
1718 
1719  // scale or gsd
1720  if (theScale.size() > 1)
1721  {
1722  ossimDpt scale (theScale[0], theScale[1]);
1723  kwl.add(prefix, ossimKeywordNames::PIXEL_SCALE_XY_KW, scale.toString(), true);
1724  kwl.add(prefix, ossimKeywordNames::PIXEL_SCALE_UNITS_KW, "degrees", true);
1725 
1726  // origin
1727  if ( ossim::isnan(theOriginLat) )
1728  {
1729  //---
1730  // Put the origin lat at the center of the image so the meters per
1731  // pixel is somewhat real.
1732  //---
1733  double centerY = theLength/2.0;
1734  theOriginLat = tieGpt.lat - theScale[1]*centerY;
1735  }
1736 
1737  if ( ossim::isnan(theOriginLon) )
1738  {
1739  theOriginLon = 0.0;
1740  }
1741  }
1742 
1744  {
1747  }
1748 
1749  }
1750  else // Projected
1751  {
1752  // Tiepoint
1753  ossimDpt tiepoint (convert2meters(x_tie_point),convert2meters(y_tie_point));
1754  kwl.add(prefix, ossimKeywordNames::TIE_POINT_XY_KW, tiepoint.toString(), true);
1755  kwl.add(prefix, ossimKeywordNames::TIE_POINT_UNITS_KW, "meters", true);
1756 
1757  // scale or gsd
1758  if (theScale.size() > 1)
1759  {
1761  kwl.add(prefix, ossimKeywordNames::PIXEL_SCALE_XY_KW, scale.toString(), true);
1762  kwl.add(prefix, ossimKeywordNames::PIXEL_SCALE_UNITS_KW, "meters", true);
1763  }
1764 
1765  // origin
1770 
1771  // std paralles for conical projections
1774 
1775  // false easting and northing.
1778 
1779  // Based on projection type, override/add the appropriate info.
1780  if (getOssimProjectionName() == "ossimUtmProjection")
1781  {
1782  // Check the zone before adding...
1783  if (theZone > 0 && theZone < 61)
1784  {
1785  kwl.add(prefix, ossimKeywordNames::ZONE_KW, theZone);
1786  }
1787  else
1788  {
1790  << "FATAL ossimGeoTiff::addImageGeometry: "
1791  << "UTM zone " << theZone << " out of range!\n"
1792  << "Valid range: 1 to 60" << endl;
1793  return false;
1794  }
1795 
1796  // Check the hemisphere before adding.
1797  if (theHemisphere == "N" || theHemisphere == "S")
1798  {
1800  }
1801  else
1802  {
1804  << "FATAL ossimGeoTiff::addImageGeometry: "
1805  << "UTM hemisphere " << theHemisphere << " is invalid!\n"
1806  << "Valid hemisphere: N or S" << std::endl;
1807  return false;
1808  }
1809 
1810  //---
1811  // Must set the central meridian even though the zone should do it.
1812  // (in decimal degrees)
1813  //---
1814  double central_meridian = ( 6.0 * abs(theZone) ) - 183.0;
1816  kwl.add(prefix, ossimKeywordNames::ORIGIN_LATITUDE_KW, 0.0, true);
1817 
1818  } // End of "if (UTM)"
1819 
1820  else if (getOssimProjectionName() == "ossimTransMercatorProjection")
1822 
1823  } // end of projected CS
1824 
1825  //---
1826  // Get the model transformation info if it's present.
1827  //---
1828  if (usingModelTransform())
1829  {
1830  std::vector<double> v = getModelTransformation();
1831  std::ostringstream out;
1832  out << std::setprecision(15); // To avoid truncating.
1833  ossim_uint32 idx = 0;
1834  for(idx =0; idx < 16; ++idx)
1835  {
1836  out << v[idx] << " ";
1837  }
1838  kwl.add(prefix, ossimKeywordNames::IMAGE_MODEL_TRANSFORM_MATRIX_KW, out.str().c_str(), true);
1839  ossimUnitType modelTransformUnitType = OSSIM_UNIT_UNKNOWN;
1840  if(theModelType == ModelTypeGeographic)
1841  {
1842  switch(theAngularUnits)
1843  {
1844  case ANGULAR_DEGREE:
1845  modelTransformUnitType = OSSIM_DEGREES;
1846  break;
1847  case ANGULAR_ARC_MINUTE:
1848  modelTransformUnitType = OSSIM_MINUTES;
1849  break;
1850  case ANGULAR_ARC_SECOND:
1851  modelTransformUnitType = OSSIM_SECONDS;
1852  break;
1853  default:
1854  return false;
1855  }
1856  }
1857  else if(theModelType == ModelTypeProjected)
1858  {
1859  switch(theLinearUnitsCode)
1860  {
1861  case LINEAR_METER:
1862  modelTransformUnitType = OSSIM_METERS;
1863  break;
1864  default:
1865  return false;
1866  }
1867  }
1869  ossimUnitTypeLut::instance()->getEntryString(modelTransformUnitType), true);
1870  }
1871 
1872  if(theScaleFactor > 0.0)
1873  {
1875  }
1876 
1877  if (traceDebug())
1878  {
1880  << "DEBUG ossimGeoTiff::addImageGeometry: Keyword list dump:\n"
1881  << kwl << std::endl;
1882  }
1883 
1884  return true;
1885 }
1886 
1887 double ossimGeoTiff::convert2meters(double d) const
1888 {
1889  switch(theLinearUnitsCode)
1890  {
1891  case LINEAR_METER:
1892  return d;
1893  case LINEAR_FOOT:
1894  return ossim::ft2mtrs(d);
1895  case LINEAR_FOOT_US_SURVEY:
1896  return ossim::usft2mtrs(d);
1897  default:
1898  if (traceDebug())
1899  {
1901  << "DEBUG ossimGeoTiff::convert2meters: "
1902  << "Linear units code was undefined!\n"
1903  << "No conversion was performed." << std::endl;
1904  }
1905  break;
1906  }
1907 
1908  return d;
1909 }
1910 
1912 {
1913  return theProjectionName;
1914 }
1915 
1917 {
1918  //---
1919  // The "parsePcsCode" method can also set the projection name. So check
1920  // it prior to looking in the lookup table.
1921  //---
1922  if (theProjectionName == "unknown")
1923  {
1924  ossimString name = COORD_TRANS_LUT.getEntryString(theCoorTransGeoCode);
1925 
1926  if (name.size())
1927  {
1928  theProjectionName = name;
1929  }
1930  }
1931 
1932  // If still unknown check for the model type.
1933  if (theProjectionName == "unknown")
1934  {
1935  if (theModelType == ModelTypeGeographic)
1936  {
1937  theProjectionName = "ossimEquDistCylProjection";
1938  }
1939  }
1940 
1941  if (traceDebug())
1942  {
1944  << "DEBUG ossimGeoTiff::setOssimProjectionName: "
1945  << "theProjectionName: "
1947  << std::endl;
1948  }
1949 
1950 }
1951 
1952 void ossimGeoTiff::setOssimProjectionName(std::shared_ptr<ossim::TiffHandlerState> state, ossim_int32 entryIdx)
1953 {
1954  ossimString value;
1955  //---
1956  // The "parsePcsCode" method can also set the projection name. So check
1957  // it prior to looking in the lookup table.
1958  //---
1959  if (theProjectionName == "unknown")
1960  {
1961  if(state->getValue(value, entryIdx, "tifftag.coord_trans_code"))
1962  {
1963  ossimString name = COORD_TRANS_LUT.getEntryString(value.toInt32());
1964 
1965  if (name.size())
1966  {
1967  theProjectionName = name;
1968  }
1969  }
1970  }
1971 
1972  // If still unknown check for the model type.
1973  if (theProjectionName == "unknown")
1974  {
1975  if (theModelType == ModelTypeGeographic)
1976  {
1977  theProjectionName = "ossimEquDistCylProjection";
1978  }
1979  }
1980 
1981  if (traceDebug())
1982  {
1984  << "DEBUG ossimGeoTiff::setOssimProjectionName: "
1985  << "theProjectionName: "
1987  << std::endl;
1988  }
1989 
1990 }
1991 
1993 {
1994  return theDatumName;
1995 }
1996 
1998 {
1999  //---
2000  // The "parsePcsCode" method can also set the datum name. So check
2001  // it prior to trying to assign.
2002  //---
2003  if (theDatumName == "unknown")
2004  {
2005  ossimString name = DATUM_LUT.getEntryString(theDatumCode);
2006 
2007  if (!name.empty())
2008  {
2009  theDatumName = name;
2010  }
2011  else
2012  {
2013  // Try the GCS code.
2014  name = DATUM_LUT.getEntryString(theGcsCode);
2015  if (name.size())
2016  {
2017  theDatumName = name;
2018  }
2019  }
2020  }
2021  if (traceDebug())
2022  {
2024  << "DEBUG ossimGeoTiff::setOssimDatumName: "
2025  << "theDatumName: "
2026  << theDatumName
2027  << std::endl;
2028  }
2029 }
2030 
2031 void ossimGeoTiff::setOssimDatumName(std::shared_ptr<ossim::TiffHandlerState> state,
2032  ossim_int32 entryIdx)
2033 {
2034  ossimString value;
2035  //---
2036  // The "parsePcsCode" method can also set the datum name. So check
2037  // it prior to trying to assign.
2038  //---
2039  if (theDatumName == "unknown")
2040  {
2041  if(state->getValue(value, entryIdx, "tifftag.datum_code"))
2042  {
2043  ossimString name = DATUM_LUT.getEntryString(value.toInt32());
2044 
2045  if (!name.empty())
2046  {
2047  theDatumName = name;
2048  }
2049  else
2050  {
2051  // Try the GCS code.
2052  if(state->getValue(value, entryIdx, "tifftag.gcs_code"))
2053  {
2054  name = DATUM_LUT.getEntryString(value.toInt32());
2055  if (name.size())
2056  {
2057  theDatumName = name;
2058  }
2059  }
2060  }
2061  }
2062  }
2063  if (traceDebug())
2064  {
2066  << "DEBUG ossimGeoTiff::setOssimDatumName: "
2067  << "theDatumName: "
2068  << theDatumName
2069  << std::endl;
2070  }
2071 }
2072 
2073 //*************************************************************************************************
2075 //*************************************************************************************************
2077 {
2078  // key 3072 Section 6.3.3.1 codes
2079  ossimString epsg_spec (ossimString("EPSG:") + ossimString::toString(thePcsCode));
2082  ossimMapProjection* map_proj = PTR_CAST(ossimMapProjection, proj.get());
2083  if (!parseProjection(map_proj))
2084  thePcsCode = 0;
2085 
2086  return (thePcsCode != 0);
2087 }
2088 
2089 //*************************************************************************************************
2091 //*************************************************************************************************
2093 {
2094  if (map_proj == NULL)
2095  return false;
2096 
2097  // Initialize parameters from base-class ossimMapProjection:
2098  if (map_proj->isGeographic())
2099  {
2100  theModelType = ModelTypeGeographic;
2103  }
2104  else
2105  theModelType = ModelTypeProjected;
2106 
2107  theProjectionName = map_proj->getProjectionName();
2108  theFalseEasting = map_proj->getFalseEasting();
2109  theFalseNorthing = map_proj->getFalseNorthing();
2110  theStdPar1 = map_proj->getStandardParallel1();
2111  theStdPar2 = map_proj->getStandardParallel2();
2112  thePcsCode = map_proj->getPcsCode();
2113 
2114  // GCS code only defined for geographic projections and is basically the same as PCS in OSSIM:
2115  theGcsCode = 0;
2116 
2117  ossimGpt origin (map_proj->origin());
2118  theOriginLat = origin.lat;
2119  theOriginLon = origin.lon;
2120 
2121  const ossimDatum* datum = map_proj->getDatum();
2122  if (datum)
2123  {
2124  theDatumName = datum->name();
2125  theDatumCode = datum->epsgCode();
2126  }
2127 
2128  // Now intercept a select few that have additional parameters specific to their derived type:
2129  if (map_proj->getProjectionName() == "ossimEquDistCylProjection")
2130  {
2131  theGcsCode = map_proj->getPcsCode();
2132  }
2133  else if (theProjectionName == "ossimUtmProjection")
2134  {
2135  ossimUtmProjection* utm_proj = PTR_CAST(ossimUtmProjection, map_proj);
2136  theHemisphere = utm_proj->getHemisphere();
2137  theScaleFactor = 0.9996; // UTM fixed
2138  theZone = utm_proj->getZone();
2139  }
2140  else if (theProjectionName == "ossimTransMercatorProjection")
2141  {
2143  theScaleFactor = tm_proj->getScaleFactor();
2144  }
2145  else if (theProjectionName == "ossimBngProjection")
2146  {
2147  // ### LEGACY HACK ###
2148  theDatumName = "OGB-M";
2149  theFalseEasting = 400000.0;
2150  theFalseNorthing = -100000.0;
2151  theScaleFactor = .9996012717;
2152  theOriginLat = 49.0;
2153  theOriginLon = -2.0;
2154  theHemisphere = "N";
2155  }
2156 
2157  return true;
2158 }
2159 
2161 {
2162  return theZone;
2163 }
2164 
2165 void ossimGeoTiff::getScale(std::vector<double>& scale) const
2166 {
2167  scale = theScale;
2168 }
2169 
2170 void ossimGeoTiff::getTiePoint(std::vector<double>& tie_point) const
2171 {
2172  tie_point = theTiePoint;
2173 }
2174 
2175 void ossimGeoTiff::getModelTransformation(std::vector<double>& transform) const
2176 {
2177  transform = theModelTransformation;
2178 }
2179 
2180 const std::vector<double>& ossimGeoTiff::getScale() const
2181 {
2182  return theScale;
2183 }
2184 
2186 {
2188  return OSSIM_PIXEL_IS_AREA;
2189 
2190  return OSSIM_PIXEL_IS_POINT;
2191 }
2192 
2193 const std::vector<double>& ossimGeoTiff::getTiePoint() const
2194 {
2195  return theTiePoint;
2196 }
2197 
2198 const std::vector<double>& ossimGeoTiff::getModelTransformation() const
2199 {
2200  return theModelTransformation;
2201 }
2202 
2204 {
2205  return theWidth;
2206 }
2207 
2209 {
2210  return theLength;
2211 }
2212 
2214 {
2215  // Capture stream flags.
2216  std::ios_base::fmtflags f = out.flags();
2217 
2218  out << setiosflags(ios::fixed) << setprecision(15)
2219  << "ossimGeoTiff::print" << std::endl;
2220 
2221  if (theScale.size())
2222  {
2223  std::vector<double>::const_iterator i = theScale.begin();
2224  ossim_uint32 index = 0;
2225  while ( i < theScale.end() )
2226  {
2227  out << "theScale[" << index << "]: " << (*i) << std::endl;
2228  ++index;
2229  ++i;
2230  }
2231  }
2232  else
2233  {
2234  out << "theScale is empty..." << endl;
2235  }
2236 
2237  if (theTiePoint.size())
2238  {
2239  std::vector<double>::const_iterator i = theTiePoint.begin();
2240  ossim_uint32 index = 0;
2241  while ( i < theTiePoint.end() )
2242  {
2243  out << "theTiePoint[" << index << "]: " << (*i) << std::endl;
2244  ++index;
2245  ++i;
2246  }
2247  }
2248  else
2249  {
2250  out << "theTiePoint is empty..." << endl;
2251  }
2252 
2253  if (theModelTransformation.size())
2254  {
2255  std::vector<double>::const_iterator i = theModelTransformation.begin();
2256  ossim_uint32 index = 0;
2257  while ( i < theModelTransformation.end() )
2258  {
2259  out << "theModelTransformation[" << index << "]: "
2260  << (*i) << std::endl;
2261  ++index;
2262  ++i;
2263  }
2264  }
2265  else
2266  {
2267  out << "theModelTransformation is empty..." << endl;
2268  }
2269 
2270  // Set the flags back.
2271  out.flags(f);
2272 
2273  return out;
2274 }
2275 
2276 
2278 {
2279  //---
2280  // If we have 16 model points do we always use them? (drb)
2281  //
2282  // In other word should the check just be if size == 16?
2283  //---
2284  if (getModelTransformation().size() == 16)
2285  {
2286  if (theScale.size() == 0)
2287  {
2288  // Need at least 24 (which is four ties) to use bilinear.
2289  if (theTiePoint.size() < 24)
2290  {
2291  return true;
2292  }
2293  }
2294  }
2295  return false;
2296 }
2297 
2299 {
2300  ossim_uint32 idx = 0;
2301  ossim_uint32 tieCount = (ossim_uint32)theTiePoint.size()/6;
2302  const double* tiePointsPtr = &theTiePoint.front();
2303  double offset = 0;
2304  if (hasOneBasedTiePoints())
2305  {
2306  offset = -1.0;
2307  }
2308 
2309  for(idx = 0; idx < tieCount; ++idx)
2310  {
2311  ossimDpt xyPixel(tiePointsPtr[0]+offset, tiePointsPtr[1]+offset);
2312  // tie[3] = x, tie[4]
2313  ossimGpt gpt(tiePointsPtr[4], tiePointsPtr[3], tiePointsPtr[5]);
2314 
2315  tieSet.addTiePoint(new ossimTieGpt(gpt, xyPixel, .5));
2316  tiePointsPtr+=6;
2317  }
2318 }
2319 
2321 {
2322  bool result = false;
2323 
2324  // Assuming ties of (x,y,z,lat,lon,hgt) so size should be divide by 3.
2325  if (theTiePoint.size()%6)
2326  {
2327  return result;
2328  }
2329 
2330  ossim_float64 minX = 999999.0;
2331  ossim_float64 minY = 999999.0;
2332  ossim_float64 maxX = 0.0;
2333  ossim_float64 maxY = 0.0;
2334 
2335  const ossim_uint32 SIZE = (ossim_uint32)theTiePoint.size();
2336  ossim_uint32 tieIndex = 0;
2337 
2338  while (tieIndex < SIZE)
2339  {
2340  if ( theTiePoint[tieIndex] < minX ) minX = theTiePoint[tieIndex];
2341  if ( theTiePoint[tieIndex] > maxX ) maxX = theTiePoint[tieIndex];
2342  if ( theTiePoint[tieIndex+1] < minY ) minY = theTiePoint[tieIndex+1];
2343  if ( theTiePoint[tieIndex+1] > maxY ) maxY = theTiePoint[tieIndex+1];
2344  tieIndex += 6;
2345  }
2346 
2347  if ( (minX == 1) && (maxX == theWidth) &&
2348  (minY == 1) && (maxY == theLength) )
2349  {
2350  result = true;
2351  }
2352 
2353 #if 0
2354  if (traceDebug())
2355  {
2357  << "ossimGeoTiff::hasOneBasedTiePoints DEBUG:"
2358  << "\nminX: " << minX
2359  << "\nmaxX: " << maxX
2360  << "\nminY: " << minY
2361  << "\nmaxY: " << maxY
2362  << "\ntheWidth: " << theWidth
2363  << "\ntheLength: " << theLength
2364  << "\none based: " << (result?"true":"false")
2365  << std::endl;
2366  }
2367 #endif
2368 
2369  return result;
2370 }
2371 
2372 //*************************************************************************************************
2373 // ArcMAP 9.2 bug workaround.
2374 // ESH 05/2008 -- ArcMap 9.2 compatibility hack
2375 // If the PCS code is for a HARN state plane and the implied pcs
2376 // code's units is feet (us or intl), we find the equivalent code
2377 // for units of meters. We're doing this because ArcMap (9.2 and
2378 // less) doesn't understand the non-meters HARN codes. However,
2379 // the units are left unchanged in this process, so the units can
2380 // be different than the user-specified pcs code. ArcMap 9.2
2381 // seems to understand the mixed definition just fine.
2382 // OLK 04/2010 -- Converted to vector<pair> scheme after refactoring EPSG factory.
2383 //*************************************************************************************************
2385 {
2386  static const ossim_uint16 harn_feet[] =
2387  {
2388  2867, 2868, 2869, 2870, 2871, 2872, 2873, 2874, 2875, 2876, 2877, 2878, 2879, 2880, 2881, 2882,
2389  2883, 2884, 2885, 2886, 2887, 2888, 2891, 2892, 2893, 2894, 2895, 2896, 2897, 2898, 2899, 2900,
2390  2901, 2902, 2903, 2904, 2905, 2906, 2907, 2908, 2909, 2910, 2911, 2912, 2913, 2914, 2915, 2916,
2391  2917, 2918, 2919, 2920, 2921, 2922, 2923, 2924, 2925, 2926, 2927, 2928, 2929, 2930, 2967, 2968
2392  };
2393  static const ossim_uint16 harn_meters[] =
2394  {
2395  2761, 2762, 2763, 2766, 2767, 2768, 2769, 2770, 2771, 2772, 2773, 2774, 2775, 2776, 2777, 2778,
2396  2779, 2780, 2781, 2787, 2788, 2789, 2798, 2799, 2804, 2805, 2806, 2807, 2808, 2809, 2813, 2814,
2397  2818, 2825, 2826, 2827, 2828, 2829, 2830, 2831, 2832, 2833, 2836, 2837, 2838, 2839, 2843, 2844,
2398  2845, 2846, 2847, 2848, 2849, 2850, 2851, 2853, 2854, 2855, 2856, 2859, 2860, 2861, 2792, 2793
2399  };
2400 
2401  ossim_uint16 result = 0;
2402  int index = 0;
2403  while ((result == 0) && (index < 64))
2404  {
2405  if (harn_feet[index] == feet_harn_code)
2406  result = harn_meters[index];
2407  ++index;
2408  }
2409  return result;
2410 };
2411 
double theScaleFactor
Definition: ossimGeoTiff.h:290
double theOriginLat
Definition: ossimGeoTiff.h:287
ossim_uint32 theWidth
Definition: ossimGeoTiff.h:271
double convert2meters(double d) const
Converts double passed in to meters if needed.
std::vector< double > theScale
Definition: ossimGeoTiff.h:265
static bool remove(const ossimFilename &pathname)
Removes pathname from filesystem if supported by platform.
std::basic_ostringstream< char > ostringstream
Class for char output memory streams.
Definition: ossimIosFwd.h:35
static const char * DATUM_KW
virtual double getFalseNorthing() const
ossim_uint32 getEpsgCode() const
double lond() const
Will convert the radian measure to degrees.
Definition: ossimGpt.h:97
static const char * CENTRAL_MERIDIAN_KW
bool addImageGeometry(ossimKeywordlist &kwl, const char *prefix=0) const
Add geometry info from tags to keword list.
ossimString beforeRegExp(const char *regularExpressionPattern) const
Returns from start of string up to but not including found pattern.
ossimUnitType
Represents serializable keyword/value map.
ossim_uint16 theAngularUnits
Definition: ossimGeoTiff.h:279
std::basic_ifstream< char > ifstream
Class for char input file streams.
Definition: ossimIosFwd.h:44
bool valid() const
Definition: ossimRefPtr.h:75
bool usingModelTransform() const
const char * find(const char *key) const
ossim_int64 fileSize() const
const std::vector< double > & getModelTransformation() const
virtual ossim_uint32 getPcsCode() const
Returns the EPSG PCS code or 32767 if the projection is a custom (non-EPSG) projection.
virtual ossimString getEntryString(ossim_int32 entry_number) const
virtual ossim_uint32 epsgCode() const
Definition: ossimDatum.h:59
double nan()
Method to return ieee floating point double precision NAN.
Definition: ossimCommon.h:135
static const char * IMAGE_MODEL_TRANSFORM_MATRIX_KW
std::vector< double > theModelTransformation
Definition: ossimGeoTiff.h:267
double y
Definition: ossimDpt.h:165
ossimDpt getDecimalDegreesPerPixel() const
Returns the decimal degrees per pixel.
virtual const ossimString & code() const
Definition: ossimDatum.h:57
static ossimString toString(bool aValue)
Numeric to string methods.
bool parseProjection(ossimMapProjection *map_proj)
Initializes data members given a projection. Returns TRUE if successful.
virtual bool isGeographic() const
virtual const ossimString & name() const
Definition: ossimDatum.h:58
ossim_uint16 theGcsCode
Definition: ossimGeoTiff.h:277
virtual bool saveState(ossimKeywordlist &kwl, const char *prefix=0) const
double acosd(double x)
Definition: ossimCommon.h:264
static const ossimErrorCode OSSIM_ERROR
virtual ossimString getClassName() const
Definition: ossimObject.cpp:64
std::vector< double > theTiePoint
Definition: ossimGeoTiff.h:266
void setPixelType(ossimPixelType type)
Sets the data member "thePixelType".
ossim_uint16 theDatumCode
Definition: ossimGeoTiff.h:278
int mapZone() const
Returns the map zone as an interger.
static std::mutex theMutex
Definition: ossimGeoTiff.h:294
virtual double getStandardParallel2() const
Derived classes should implement as needed.
virtual double getB() const
const std::vector< double > & getScale() const
TIFF * theTiffPtr
Definition: ossimGeoTiff.h:254
unsigned short ossim_uint16
#define abs(a)
Definition: auxiliary.h:74
double latd() const
Will convert the radian measure to degrees.
Definition: ossimGpt.h:87
void wrap()
Wrap method to maintain longitude between -180 and +180 and latitude between -90 and +90...
Definition: ossimGpt.h:305
static const char * TYPE_KW
ossimString getOssimDatumName() const
Returns an ossimString representing the ossim datum name code.
ossimPixelType getPixelType() const
Returns data member "thePixelType".
virtual double optimizeFit(const ossimTieGptSet &tieSet, double *targetVariance=0)
virtual double getStandardParallel1() const
Derived classes should implement as needed.
double theStdPar2
Definition: ossimGeoTiff.h:285
ossim_uint16 theLinearUnitsCode
Definition: ossimGeoTiff.h:283
ossim_int32 toInt32() const
ossimString afterRegExp(const char *regularExpressionPattern) const
Returns from position after found pattern to end of string.
static const char * PIXEL_SCALE_XY_KW
void setOssimDatumName()
Attempts to set the ossim datum code.
double theFalseEasting
Definition: ossimGeoTiff.h:288
double ossim_float64
void add(const char *prefix, const ossimKeywordlist &kwl, bool overwrite=true)
double theStdPar1
Definition: ossimGeoTiff.h:284
static const char * ZONE_KW
ossimProjection * createProjection(const ossimFilename &filename, ossim_uint32 entryIdx) const
ossim_uint16 theModelType
Definition: ossimGeoTiff.h:275
yy_size_t size
virtual const ossimDatum * getDatum() const
static const char * TIE_POINT_XY_KW
ossim_uint16 theCoorTransGeoCode
Definition: ossimGeoTiff.h:282
static const char * FALSE_NORTHING_KW
std::string::size_type size() const
Definition: ossimString.h:405
static ossimEpsgProjectionDatabase * instance()
Instantiates singleton instance of this class:
static bool writeJp2GeotiffBox(const ossimFilename &tmpFile, const ossimIrect &rect, const ossimProjection *proj, std::vector< ossim_uint8 > &buf, ossimPixelType pixelType)
Writes a geotiff box to a buffer.
ossim_uint16 theRasterType
Definition: ossimGeoTiff.h:276
virtual ossimString getProjectionName() const
Returns the projection name.
virtual ossimGpt origin() const
ossim_uint16 getMetersEquivalentHarnCode(ossim_uint16 feet_harn_code)
unsigned int ossim_uint32
virtual const ossimEllipsoid * ellipsoid() const
Definition: ossimDatum.h:60
double toDouble() const
#define PTR_CAST(T, p)
Definition: ossimRtti.h:321
virtual bool saveState(ossimKeywordlist &kwl, const char *prefix=0) const
Method to save the state of an object to a keyword list.
ossimString getOssimProjectionName() const
Returns an ossimString representing the ossim projection name.
virtual std::ostream & print(std::ostream &out) const
Prints data members.
ossimGpt ulGroundPt() const
Returns the upper left ground point.
static const char * STD_PARALLEL_1_KW
static const char * FALSE_EASTING_KW
ossimString theAsciiParam
Definition: ossimGeoTiff.h:269
struct tiff TIFF
ossimPixelType getRasterType() const
int getWidth() const
static const char * IMAGE_MODEL_TRANSFORM_UNIT_KW
bool parsePcsCode()
Initializes data members given a projection code.
double ft2mtrs(double feet)
Definition: ossimCommon.h:372
const ossimMapProjection * getProjection() const
Returns reference to "theProjection".
storage class for tie point between ground and image based on ossimGpt
Definition: ossimTieGpt.h:24
ossimDpt getMetersPerPixel() const
Returns the pixel size in meters.
static ossimProjectionFactoryRegistry * instance()
ossim_uint32 theLength
Definition: ossimGeoTiff.h:272
static const char * ORIGIN_LATITUDE_KW
const std::vector< double > & getTiePoint() const
ossimString thePcsCitation
Definition: ossimGeoTiff.h:281
virtual double getFalseEasting() const
ossimString theHemisphere
Definition: ossimGeoTiff.h:259
bool hasOneBasedTiePoints() const
Attempts to detect if tie points are one or zero based.
ossimDpt ulEastingNorthingPt() const
Returns the upper left easting and northing as a ossimDpt.
void addTiePoint(ossimRefPtr< ossimTieGpt > aTiePt)
operations
short ossim_int16
ossimPixelType
ossimString toString(ossim_uint32 precision=15) const
Definition: ossimDpt.cpp:160
double theFalseNorthing
Definition: ossimGeoTiff.h:289
ossimUnitType getProjectionUnits() const
OSSIM considers all map projection coordinates (including false eastings and northings) to be in mete...
double usft2mtrs(double feet)
Definition: ossimCommon.h:373
storage class for a set of geographic tie points, between master and slave images ...
void setOssimProjectionName()
Attempts to set the ossim projection name from keys read.
static const char * HEMISPHERE_KW
double x
Definition: ossimDpt.h:164
const char * c_str() const
Returns a pointer to a null-terminated array of characters representing the string&#39;s contents...
Definition: ossimString.h:396
bool empty() const
Definition: ossimString.h:411
virtual ossimProjection * createProjection(const ossimFilename &filename, ossim_uint32 entryIdx) const
STUB. Not implemented.
virtual bool saveState(ossimKeywordlist &kwl, const char *prefix=0) const
ossim_int32 getZone() const
std::vector< double > theDoubleParam
Definition: ossimGeoTiff.h:268
void getTieSet(ossimTieGptSet &tieSet) const
Initializes tieSet from theTiePoints.
int getLength() const
ossim_float64 lat
Definition: ossimGpt.h:265
double mtrs2usft(double meters)
Definition: ossimCommon.h:375
static const char * SCALE_FACTOR_KW
double mtrs2ft(double meters)
Definition: ossimCommon.h:374
static const char * STD_PARALLEL_2_KW
bool theGeoKeysPresentFlag
Definition: ossimGeoTiff.h:257
bool readTags(const ossimFilename &file, ossim_uint32 entryIdx=0)
Reads tags.
ossimGeoTiff()
default constructor
static ossimUnitTypeLut * instance()
Returns the static instance of an ossimUnitTypeLut object.
static bool writeTags(TIFF *tiffOut, const ossimRefPtr< ossimMapProjectionInfo > projectionInfo, bool imagineNad27Flag=false)
double central_meridian(double xmin, double xmax)
virtual double getA() const
ACCESS METHODS:
unsigned char ossim_uint8
static ossimEpsgProjectionFactory * instance()
Implements singleton pattern.
OSSIMDLLEXPORT std::ostream & ossimNotify(ossimNotifyLevel level=ossimNotifyLevel_WARN)
ossimString theDatumName
Definition: ossimGeoTiff.h:263
Projection Factory for coded projections defined in database.
std::basic_ostream< char > ostream
Base class for char output streams.
Definition: ossimIosFwd.h:23
static const char * PIXEL_SCALE_UNITS_KW
ossim_uint32 thePcsCode
Definition: ossimGeoTiff.h:280
static const char * TIE_POINT_UNITS_KW
ossimString theProjectionName
Definition: ossimGeoTiff.h:262
int ossim_int32
static int getPcsUnitType(ossim_int32 pcsCode)
double theOriginLon
Definition: ossimGeoTiff.h:286
ossimPrivateGtifDef * thePrivateDefinitions
Definition: ossimGeoTiff.h:292
bool isnan(const float &v)
isnan Test for floating point Not A Number (NAN) value.
Definition: ossimCommon.h:91