Georeferenced 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 geoimages is GeoTIFF. GeoTIFF images are just regular TIFF (tagged image file format) images with some additional tags embedded within the file. There are other geoimage formats: Sometimes, JPEG images come with a jpeg world ("JGW") file that provides georeferencing information. This article focuses on the GeoTIFF format and provides some conceptual background and a sample implementation for visualizing geoimages. 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. Georeferencing 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 geoimages, 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 georeferencing 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–firstperson 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 documentview 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 DeviceIndependent 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:
 Working with the Map Projection Library
 Map projections
 Using the map projection library
 Forward projections
 Inverse projections
 Projections without the library
 Working with the GeoTIFF Library
 Common tags and variants
 Displaying Georeferenced Images
 2D image display
 3D image display
 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 semitransparent 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.
Working with the GeoTIFF Library
Common Tags
Extracting the georeferencing 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 upperleftmost 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 georeferencing 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 georeferencing 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 georeferencing 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, >_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 < (argSize1)) && (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;
}
Displaying GeoReferenced Images
2D Image Display
GeoTIFF images are laid out such that the upper left (north west) corner has the UTM coordinates equal to the tie points obtained in the code snippet above. Moving left to right the X coordinate (easting) increases by the number of pixels times the Xresolution. Moving top to bottom, the Y coordinate decreases by the Yresolution for every pixel. So, if the Yresolution is 10 meters and the mouse is 10 pixels below the top of the image Y = tieY  10 x 10 = tieY  100 meters. UTM coordinate are usually measured in meters. Once you have the tie points, the UTM zone, and the resolution, it's relatively simple to set up a projection using the projection library. So, the projection library takes in (X,Y) coordinates measured in meters and spits out (latitude, longitude) measured in radians (in the sample, I convert to degrees). So, taking the above computations to get UTM (X,Y) for each point on the image and the projection library, you now have a way to calculate (latitude, longitude) for each pixel on the screen.
Just to recap: To calculate the (X, Y) position of the mouse on a UTMprojected GeoTIFF, use these calculations:
X = TIE_POINT.X + PIXEL_SCALE.X * MOUSE_POINT.X
Y = TIE_POINT.Y  PIXEL_SCALE.Y * MOUSE_POINT.Y
From this (X, Y) point, use the UTMProjection class above to calculate the latitude and longitude of the point. The UTMProjection class needs to be initialized from the tags within the TIFF image (see the code snippet above).
3D Image Display
Displaying a georeferenced image in 3D takes slightly more effort than a 2D representation. The extra effort is really in two areas: coordinate transformations and texture mapping.
Coordinate transformations pick up where the 2D coordinate transformations leave off. For the 2D case, you calculated the UTM coordinates (X, Y) for each pixel, and then used the projection class to calculate (latitude, longitude) for each pixel...now what? Well, latitude and longitude coordinates are actually part of a spherical coordinate system. Normally, a spherical coordinate system is expressed in terms of the Greek letters r, F, q (Rho, phi, and theta). Rho is the radius and it normally omitted in geographicspherical coordinates because it's understood to be the radius of the Earth, phi is simply longitude, and theta is latitude. So, you can take each point on the image convert to UTM coordinates, then to spherical coordinates (in other words, latitude and longitude). Now, youneed to convert from spherical to 3D cartesian coordinates so that you can lay out your 3D Earth model with the GeoTIFF image textured onto the correct portion of the model. Given latitude and longitude (q, F), you calculate the (x, y, z) position of a point the following way:
x = R * cos((F)) * cos((q));
y = R * sin((q));
z = R * sin((F)) * cos((q));
Because of the orientation of the (x, y, z) coordinate system used by OpenGL, it's necessary to take the negative of the longitude in these calculations. If you care to visualize this image, you are in space looking down at 0° longitude. To your left are negative longitudes and to the right positive. However, for an OpenGL image, you are looking down the X axis towards the origin and the Z axis points to your left, intersecting the Earth at 90° longitude (west longitude  south of Texas). So basically, the Z coordinate axis is reversed for the normal longitude convention, which is to take western hemisphere longitudes as negative and eastern longitudes as positive. This reversal from normal spherical coordinate transformations can be confusing if you're not careful. It's tempting to write the calculations wherever you need them in your code, but I suggest you set up a single function for transformation and use it everywhere. Here is the one I used in the sample project:
void Earth3D::GeoToXYZ(double lat, double lon, double& x, double& y,
double& z, double R)
{
x = R * cos((lon*GEO::DE2RA)) * cos((lat*GEO::DE2RA));
y = R * sin((lat*GEO::DE2RA));
z = R * sin((lon*GEO::DE2RA)) * cos((lat*GEO::DE2RA));
}
Another complication of moving to 3D coordinates is texture mapping. Texture mapping GeoTIFF images can be confusing because the GeoTIFF tie points are at the upper right corner of the image but OpenGL texture coordinates (for 2D textures) considers the lower left corner of the image to be the origin (0, 0) increasing to (1, 1) at the upper right corner of the image. Also, it helps to keep in mind that a texturemapped point is like a spot weld between the texture image and the underlying model (triangle mesh). If the spot welds are too far apart, the areas in between may distort slightly when you look closely. Texture mapping each and every pixel will reduce your rendering speed so it's a trade off depending on the application. A more interesting approach is mipmapping, which essentially uses different textureimage resolutions depending on the scale of the model. Unfortunately, there isn't time or space to cover that here.
The viewing approach used in the sample project essentially builds the model around the origin, with the actual Earth radius corresponding to an OpenGL model radius of 1.0. Viewing orientation is accomplished by using the convenient OpenGL method: gluLookAt. This method allows you to specify the spot on the model (or anywhere) that you wish to view (in the sample, I choose the center of the image), the position where you are viewing from, and an up vector. OpenGL takes that information and conveniently transforms the coordinates so that you get exactly what you want. It's great for firstperson fly throughs or just general viewing because you can just draw your model wherever you like, without having to offset any of your objects. Here is a code snippet where I set the "look" information:
void Earth3D::SetOrientation(void)
{
if (m_flyMode && m_flyFirstPerson)
{
gluLookAt(m_flyPt[0], m_flyPt[1], m_flyPt[2],
m_flyTo[0], m_flyTo[1], m_flyTo[2],
m_flyTo[0], m_flyTo[1], m_flyTo[2]);
}
else
{
double upVec[3];
// dividing by the Earth radius essentially scales from a world
// coordinate system to this OpenGL coordinate system where the
// radius of the Earth equals 1.0
double D = (m_eyeDist / m_model_matrix.XScale())/GEO::ERAD;
V3D p(m_unitSouthPt);
V3D u(m_unitCenterPt);
V3D q, v, w, x;
// rotate the south vector around the center vector
ArbitraryRotate(p,(m_model_matrix.YRot())*GEO::DE2RA,u,q);
v.Create(m_unitEastPt);
// rotate the east vector around the center vector
ArbitraryRotate(v,(m_model_matrix.YRot())*GEO::DE2RA,u,w);
p = q;
// rotate the (formerly) south vector around the (formerly) east vector
ArbitraryRotate(p,m_model_matrix.XRot()*GEO::DE2RA,w,q);
// calculate the eye point
m_eyePt[0] = m_centerPt[0] + D * q.x;
m_eyePt[1] = m_centerPt[1] + D * q.y;
m_eyePt[2] = m_centerPt[2] + D * q.z;
p.Create(m_unitCenterPt);
// rotate the center vector around the (formerly) east vector
ArbitraryRotate(p,m_model_matrix.XRot()*GEO::DE2RA,w,q);
upVec[0] = q.x;
upVec[1] = q.y;
upVec[2] = q.z;
gluLookAt( m_eyePt[0], m_eyePt[1], m_eyePt[2],
m_centerPt[0], m_centerPt[1], m_centerPt[2],
upVec[0], upVec[1], upVec[2]);
}
}
So, the eye point is calculated by going some azimuth and distance from the center of the image and then going up in altitude. The azimuth is initialized to 180 degrees, making the initial eye point due south from the geographic center of the image.
[scrn6.jpg]
Additional Information about the Sample Project
Create a folder for the sample files, download the sample files, and then unzip them all into the folder you created. That should insure that you have the correct folder structure to build everything. First, build the TIFF libraries, then projection, then JPEG, then GeoTIFF, and finally the sample project.
In the vc6\GTIFFSample folder, there is a sample GeoTIFF file for a small area in Southern California. Use the sample to try the fly through. Keep in mind that the sample project will not work with map projections other than UTM and only do 3D viewing for GeoTIFF images.
The file Earth3D.cpp has a useful method called ArbitraryRotate that rotates a vector about an arbitrary axis. This method was taken from a book called 3D Math Primer for Graphics and Game Development by Fletcher Dunn and Ian Parberry. It's a great book for understanding the mathematics behind 3D graphics!
For the sample project, I really wanted to show the flythrough capabilities that you see on the cable news shows. It really is not that difficult to show...the stuff you learned about above is really all you need to understand: coordinate transformations and texture mapping.
To do a fly through, load the sample map, or any UTM GeoTIFF, into the application and while you are in 2D mode, click the flag button on the toolbar and click a point near one of the corners of the image. Then, click the flag again and select a point at the opposite corner of the image. Make the altitude for the first point 10,000 meters and the altitude for the second point 500 meters (flying nosedown so you can see the image more clearly). Then, select View>3D, then Tools>Fly Click Points. If you choose the First Person check, you will experience flying the click points; if not, you will see a little airplane (made of three triangles) fly the click points.
Otherwise, the sample application behaves as you might expect. Keep in mind that you can only go 3D on GeoTIFF images, but the application will load other types of images. Also, you can fly more than two click points. The application is intended to explain the concepts; its purpose is not for people to hand in as homework, so no whining about functionality.