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

Displaying Geo-Referenced 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 X-resolution. Moving top to bottom, the Y coordinate decreases by the Y-resolution for every pixel. So, if the Y-resolution 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 UTM-projected 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 geo-referenced 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 geographic-spherical 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 texture-mapped 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 mip-mapping, which essentially uses different texture-image 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 first-person 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.

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 fly-through 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 nose-down 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.

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.