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

Geo-referenced images are regular map images, aerial photos, or satellite imagery with geographic referencing information embedded within them or contained within an external file. They are used for a variety of purposes, such as planning, navigation, and so forth. A common format for geo-images is GeoTIFF. GeoTIFF images are just regular TIFF (tagged image file format) images with some additional tags embedded within the file. There are other geo-image formats: Sometimes, JPEG images come with a jpeg world (“JGW”) file that provides geo-referencing information. This article focuses on the GeoTIFF format and provides some conceptual background and a sample implementation for visualizing geo-images. In a nutshell, the extra tags in a GeoTIFF file enable us to determine the precise geographic coordinates of each and every pixel in the image. Knowing the geographic layout of the pixels makes it easy to show data overlaid on the map or to show the map and data overlaid on a 3D perspective view of the Earth. Geo-referencing also makes it possible to do virtually realistic 3D fly throughs, as demonstrated in the sample application.

Unfortunately, this subject can be a little unwieldy because, to visualize geo-images, you need to understand image manipulation, map projections, geographic coordinate systems, and 2D & 3D graphics. Normally, that would be too much to cover in one article, but fortunately you can build on previous articles covering these topics. So, this article fills in some of the details of map projections and GeoTIFF files not covered in other CodeGuru articles.

The goal is to implement the following capabilities:

  • Load a GeoTIFF image file into a DIB section and decode the geo-referencing information
  • Show the image in 2D and display the geographic and UTM coordinates of the mouse location and enable the use to select points
  • Show a 3D perspective image of the GeoTIFF and enable the user to select points on the 3D view
  • Display a first person and non–first-person fly through of the image, similar to those shown on cable news shows

There are only four requirements, but this really is quite a bit to put in one example. If you look at the sample program, you will notice that there is a lot going on in there. I’ve tried to keep it as simple and straightforward as possible but considering the goals its not as clean as I would like. However, it is essentially a document-view architecture MFC application. The view class handles the 2D image viewing while two separate classes handle the 3D rendering: GLView3D and Earth3D. GLView3D is a basic OpenGL class that draws output to DIBs. Making use of some virtual functions the Earth3D class provides the GeoTIFF related 3D rendering. DIBSection (implemented in dibsect.h & dibsect.cpp) encapsulates a DIB section. The files: dibtiff.(h+cpp), dibbmp.(h+cpp) and dibijg.(h+cpp) handling reading and writing of TIFF, BMP, and JPEG images, respectively. These files and their corresponding libraries where all demonstrated and explained in a previous article: Working with TIFF Images. GLView3D would be familar if you read the article: Drawing and Printing OpenGL Graphics Using Device-Independent Bitmaps. Finally, some of the geographic coordinate systems code may be familar to if you read the article: Geographic Distance and Azimuth Calculations.

So, this article covers the following topics:

  1. Working with the Map Projection Library
    1. Map projections
    2. Using the map projection library
    3. Forward projections
    4. Inverse projections
    5. Projections without the library
  2. Working with the GeoTIFF Library
    1. Common tags and variants
  3. Displaying Geo-referenced Images
    1. 2D image display
    2. 3D image display
  4. Additional information about the sample project

Working with the Map Projection Library

Map Projections

A map projection is a systematic representation of a curved surface, typically the surface of the Earth, on a plane. That means taking coordinates referenced to the Earth (in other words, latitude, longitude, and sometimes height above mean sea level), and systematically generating horizontal and vertical (x and y) coordinates. This process is necessary because the surface of a globe cannot be represented on a flat map without some distortion, and paper maps are much more convenient to carry around than globes. Can you imagine the size of globe you would need to accurately portray city street information? The term projection comes into play because the approach for converting coordinates from one system to another (globe to paper map) is basically like projecting light through a semi-transparent globe and onto an opaque paper map. This brings up two questions: 1) what direction is the light shining from, and 2) how to position the paper with respect to the globe. Keep in mind there is not much you can do to change the geometry of a globe, about the only choice you have is to rotate it and after rotating the globe you still have the same basic shape. Paper maps, on the other hand, can be rolled to make a cylinder or cone, or they can be rotated to different orientations.

Figure 1: Geometries for conical projections.

Figure 2: Geometries for cylindrical projections.

Figure 3: Geometries for azimuthal projections.

Focusing on the cylindrical projections of Figure 2; the first image (in Figure 2) is essentially the model used for the famous Mercator projection and the second image (in Figure 2) shows the geometry of the Transverse Mercator projection. The main difference being the Mercator cylinder touches the equator everywhere making the distortion near the equator minimal, and the Transverse Mercator touches the central meridian of the projection everywhere minimizing the distortion close to a particular meridian. The Universal Transverse Mercator is a specific implementation of the Transverse Mercator projection. UTM partitions meridians into 60 zones each 6 degrees of longitude in width, running from the South to the North poles. Zone 1 has central meridian 177° W going East and covers everything North and South of the equator from 180° W to 174° W longitude. Washington, D.C. (approximately 77° W) is located in zone 18 which runs from 78° W to 72° W with central meridian 75° W. Coordinates that have been projected into UTM are generally referred to as northing and easting. Northing is essentially the distance in meters along the transverse cylinder (Figure 2, second image) from the equator. Easting is the horizontal distance along the transverse cylinder measured from 500,000 meters west of the central meridian. So if a longitude is exactly on the central meridian and gets converted to UTM its easting value should be 500,000. Why 500,000? Why not zero for the central meridian? Not sure, but possibly its to keep the northing and easting numbers positive. Three degrees from the central meridian corresponds to approximately 256,000 meters at 40° latitude and roughly 333,000 meters at the equator. So, generally speaking you will only see northing and easting numbers for UTM between about 160,000 and 840,000. Beyond those limits the coordinates would most likely be referenced to the next 6° zone. So to specify a point on the Earth in UTM you need to provide a northing, easting, and zone (technically speaking you need to specify an ellipsoid also). Sometimes you may see the zone omitted, this is because the zone is “understood,” usually by everyone except the guy looking at the numbers. Generally if a map of Washington, D.C. uses UTM projection then zone 18 is understood. Using another zone is theoretically possible but it would lead to much too much distortion.

So why so much coverage of UTM? It is used extensively and especially because GeoTIFF images from the USGS (United States Geological Survey) are provided in UTM projection, so its helpful to understand the projection.

Using the Map Projection Library

Luckily for us programmers, a free map projections library, known as proj, is freely available from www.remotesensing.org. The sample project has a stripped down version of the projections library. The wonderful thing about this library is that it is so simple to use, all you have to do is provide the library with the central meridian and zone (and ellipsoid), then you can call the forward and inverse functions as much as you like.

The forward projection function converts geographic coordinates to northing and easting. The inverse function converts northing and easting coordinates to geographic coordinates. Here is a little UTM wrapper around the projections library that demonstrates how to initialize and make forward/inverse calls.

// utmproj.h
//

#ifndef UTM_H_
#define UTM_H_

namespace GEO {
   const double PI                = 3.14159265359;
   const double TWOPI             = 6.28318530718;
   const double DE2RA             = 0.01745329252;
   const double RA2DE             = 57.2957795129;
   const double ERAD              = 6378.135;
   const double ERADM             = 6378135.0;
   const double AVG_ERAD          = 6371.0;
   const double FLATTENING        = 1.0/298.257223563;    // Earth
                                                          // flattening
                                                          // (WGS '84)
   const double EPS               = 0.000000000005;
   const double KM2MI             = 0.621371;
   const double GEOSTATIONARY_ALT = 35786.0;    // km
}

#include <projects.h>
#include <proj_api.h>

class UTMProjection {
   protected:
   PJ * m_pj;
public:
   UTMProjection(void)
   {
      m_pj = 0;
   }
   ~UTMProjection(void)
   {
      if (m_pj)
      {
         pj_free(m_pj);
      }
      m_pj = 0;
   }
   void forward(double lat, double lon, double& x, double& y)
   {
      projUV sUV;

      sUV.u = lon * GEO::DE2RA;
      sUV.v = lat * GEO::DE2RA;

      sUV = pj_fwd( sUV, m_pj );

      x = (int)(sUV.u + 0.50);
      y = (int)(sUV.v + 0.50);
   }

   void inverse(double x, double y, double& lat, double& lon)
   {
      projUV sUV;

      sUV.u = x;
      sUV.v = y;

      sUV = pj_inv( sUV, m_pj );

      lon = sUV.u * GEO::RA2DE;
      lat = sUV.v * GEO::RA2DE;
   }
   void setup(INT32 zone)
   {
      char pjData[4][24];
      char * pjPtrs[4];

      strcpy(pjData[0],"proj=utm");
      sprintf(pjData[1],"zone=%d",zone);
      strcpy(pjData[2],"ellps=clrk66");
      strcpy(pjData[3],"units=m");

      pjPtrs[0] = pjData[0];
      pjPtrs[1] = pjData[1];
      pjPtrs[2] = pjData[2];
      pjPtrs[3] = pjData[3];

      m_pj = pj_init(4,pjPtrs);
   }
};

#endif

Conic projections work a little differently than cylindrical projections like UTM due to the geometry of the surface that geographic information is being projected onto. Similar to the cylindrical projections, conic projection require information about where the flat surface (an imaginary paper map rolled into a cone or cylinder) intersects the Earth, but unlike the cylindrical case, the conical geometry may intersect the surface in one place, two places, or not at all. Typical map projections use two intersections between the cone and the Earth and those intersections are designated by two latitudes.

More by Author

Get the Free Newsletter!

Subscribe to Developer Insider for top news, trends & analysis

Must Read