2D & 3D Visualization Techniques for Geo-Referenced Images

Working with the GeoTIFF Library

Common Tags

Extracting the geo-referencing tags from a TIFF file is nearly as simple as getting an image's width and height. To get the resolution and tie points, simple calls to TIFFGetField are used with the tags: GTIFF_PIXELSCALE and GTIFF_TIEPOINTS. Generally speaking, these two tags tell you the projected coordinates for the upper-leftmost pixel on the image and how many meters of ground distance correspond to each pixel (in each direction).

Other useful pieces of information are the map projection details. The most common map projection (in GeoTIFF files) seems to be UTM, although it wouldn't be surprising to see Lambert Conformal Conic. There seem to be two common ways for obtaining map projection information: 1) by calling GTIFKeyInfo and GTIFKeyGet with the key = GTCitationGeoKey, or 2) calling TIFFGetField with the tag = GTIFF_ASCIIPARAMS. The following code is intended to work for both format variants. Unfortunately, the projection information is usually provided in the form of a string like: UTM Zone 17 N WGS84. GeoTIFF is a wonderful standard enabling geo-referencing tasks, but it's unfortunate that the standard allows different string formats for projection information. It's unfortunate that the GeoTIFF standard does not use some form of XML for the geo-referencing information. Another silly variation that you may run across is that often USGS GeoTIFF images say WGS84 in their projection string but the horizontal datum actually used is NAD27, which implies the Clarke 1866 ellipsoid. If anyone understands what the deal is, please educate me.

// GEORef is useful for storing geo-referencing information about an image
class GEORef {
   void initGeoRef(void)
   {
   }
public:
   enum Projections {
      UnknownPrj,
      UTM
   };
   enum Ellipsoids {    // some common ellipsoids
      UnknownEllipsoid,
      Clarke1866,
      GRS80,
      IAU76,
      AIRY,
      WGS72,
      WGS84,
      Sphere
   };
   double m_tieX;
   double m_tieY;
   double m_resX;
   double m_resY;
   INT32 m_zone;
   Projections m_proj;
   Ellipsoids m_ellipse;
   GEORef(void)
   {
      m_tieX = m_tieY = m_resX = m_resY = 0.0;
      m_zone = 0.0;
      m_proj = UnknownPrj;
      m_ellipse = UnknownEllipsoid;
   }
   GEORef(const GEORef& gr)
   {
      m_tieX = m_tieY = m_resX = m_resY = 0.0;
      m_zone = 0.0;
      m_proj = UnknownPrj;
      m_ellipse = UnknownEllipsoid;
      Copy(gr);
   }
   ~GEORef(void){}
   void Copy(const GEORef& gr)
   {
      m_tieX = gr.m_tieX;
      m_tieY = gr.m_tieY;
      m_resX = gr.m_resX;
      m_resY = gr.m_resY;
      m_zone = gr.m_zone;
      m_proj = gr.m_proj;
      m_ellipse = gr.m_ellipse;
   }
   GEORef& operator = (const GEORef& gr)
   {
      Copy(gr);
      return *this;
   }
};

bool ReadGeoInfo(const char * filename, GEORef& gr)
{
   bool result = false;

   TIFF * tiff = 0;
   GTIF * gtif = 0;

   tiff = XTIFFOpen(filename,"r");
   if (tiff)
   {
      gtif = GTIFNew(tiff);
      if (gtif)
      {
         const int key_size = 2048;
         char key_data[key_size];
         memset(key_data,0,(size_t)key_size);
         tagtype_t cit_type;
         int cit_size;
         INT32 ellps = 1;

         string geoCitationKey;

         double * d_list = 0;
         int d_list_count = 0;

         int cit_len = GTIFKeyInfo(gtif,GTCitationGeoKey, &cit_size,
                                   &cit_type);
         if (cit_len)
         {
            if (GTIFKeyGet(gtif,GTCitationGeoKey,key_data,0,cit_len))
            {
               geoCitationKey = (char *)key_data;
            }
         }
         if (geoCitationKey.length() == 0)
         {
            char * gt_ascii = 0;

            if (TIFFGetField(tiff, GTIFF_ASCIIPARAMS, &gt_ascii))
            {
               geoCitationKey = (char *)gt_ascii;
            }
         }

         if (geoCitationKey.length() > 0)
         {
            UtilityParser p;
            p.ParseString(geoCitationKey);

            if (p.GetSize() > 0)
            {
               vector<string> * pArgs = p.GetData();

               UINT32 argSize = pArgs->size();

               bool utmproj = false;

               for (UINT32 i = 0; i < argSize; i++){
                  string iArg = (*pArgs)[i];
                  StringUtil::Uppercase(iArg);

                  if (iArg == "UTM")
                  {
                     utmproj = true;
                  }

                  if ((iArg == "ZONE") && (i < (argSize-1)) &&    (utmproj))
                  {
                     iArg = (*pArgs)[i+1];

                     char cZone[3];
                     cZone[0] = iArg.c_str()[0];
                     cZone[1] = iArg.c_str()[1];
                     cZone[2] = 0;

                     int zone = atoi(cZone);

                     gr.m_proj = GEORef::UTM;
                     gr.m_zone = zone;
                     result = true;
                  }

                  if (StringUtil::StartsWith(iArg,"WGS84"))
                  {
                     gr.m_ellipse = GEORef::WGS84;
                  }
                  if (StringUtil::StartsWith(iArg,"NAD27"))
                  {
                     gr.m_ellipse = GEORef::Clarke1866;
                  }
               }
            }
         }

         if (TIFFGetField(tiff, GTIFF_TIEPOINTS, &d_list_count,
                          &d_list))
         {
            gr.m_tieX = d_list[3];
            gr.m_tieY = d_list[4];
         }
         if (TIFFGetField(tiff, GTIFF_PIXELSCALE, &d_list_count,
                          &d_list))
         {
            gr.m_resX = d_list[0];
            gr.m_resY = d_list[1];
         }

         XTIFFClose(tiff);
         GTIFFree(gtif);
      }
   }

   return result;
}

Andy McGovern - Software Developer/Engineer. Special interests: astronomy, image processing, orbital mechanics, mathematics. Currently work at the Johns Hopkins University Applied Physics Laboratory on the science operations team for the CRISM instrument on the Mars Reconnaissance Orbiter.