Scientific & Numeric Applications

Map Projections in C++

David T. Lowerre

David Lowerre is a Senior Engineer at Interstate Electronics, where he works primarily on navigation systems. He has been programming in C for ten years, and is also a C++ enthusiast. He can be reached via e-mail at 71562.2766@compuserve.com or Stives@ao1.com.

Introduction

Since before the days of ancient Rome, people have relied on maps. However, it wasn't until the age of exploration that they understood widely that the world really isn't flat, and therefore cannot be shown on a flat surface without some kind of error. Because maps were now being used for navigation, the nature of this error was of no small concern.

The technique of transferring points from a shape that is not flat to a plane is called projection. When points are projected onto a plane, some kind of distortion always results. Since those early days of exploration, several projections have been developed to deal with this distortion.

Different projections distort different things. For example, some preserve the shape of land masses, while others preserve the area covered by these masses. No projection can do both. A map that preserves shapes and angles is conformal and is useful for navigation. An equal-area map is useful for teaching school children that Greenland is not larger than South America!

So a program that works with maps (either drawing them or using them as input) must be able to handle more than one kind of projection. Since all projections do share some common features, it makes good sense to implement the code to perform them as a C++ class library. This article will establish the framework for just such a class library, and then use that framework to implement the most popular projection of them all: the Mercator projection.

Vectors

A point in a two-dimensional coordinate system requires two components to represent it, one for each dimension. I prefer to work with coordinates in vectors — data structures holding coordinate pairs — rather than having separate variables for each dimension.

The header file vector.hpp, shown in Listing 1, is a stripped-down implementation of a more complete vector handling class. For our purposes here it is used to make it easy to access the desired vector components. The derived classes, Spherical and Planar, are for distinguishing between global and local coordinates.

Using these classes, the following syntax is used for accessing the components of the two types of vector:

Spherical ll;
ll[LAT] = 32.733333;
ll[LON] = -117.25;

Planar xy;
xy[X] = 1.2;
xy[Y] = 5.04;

The nice thing about C++ is that by going through all this bother we have given the compiler what it needs to keep us from making mistakes like the following:

xy[LON] = 1.2;     // index is invalid

Finally, the type definitions for Degrees and Inches are to help me remember which units I am working with. Unfortunately, it is still up to me to keep this straight since the compiler cannot tell the difference.

This is probably a good place to mention sign conventions and units. Spherical coordinates are defined to be in degrees. Latitude is positive above the equator and negative below. Longitude is positive east of Greenwich, England and negative to the west. X and Y values are computed in inches. This makes it easy to convert to and from points on a chart. Also, display and print resolution is usually expressed in Dots Per Inch (DPI), which makes it easy to scale X/Y points for display.

The Projection Class

Projections come in two flavors: the projection for the sphere and the projection for the ellipsoid. A sphere is very easy to work with mathematically, but the earth is not a sphere. It is an ellipsoid. So the math becomes more difficult and we must concern ourselves with the exact shape of the ellipsoid.

An ellipsoid is defined by its polar and equatorial radii. In 1866, Alexander Ross Clarke measured these for the earth to be 3432.275 nautical miles for the polar radius, and 3443.950 nautical miles for the equatorial radius. These measurements are still widely used. More recent measurements by the World Geodetic System in 1972 and 1984 are most commonly used for nautical charts.

These measurements define the geoidal model and are common to all projections. The Projection class constructor takes the model as its only parameter and uses it to compute the eccentricity of the ellipsoid, and its square, and keeps these as protected data members of the class so that they will be available for future computations:

class Projection {
protected:
   int Model;
// eccentricity and eccentricity squared
   double e, e2;
   Planar xyResult;
   Spherical llResult;
public:
   Projection( int model );
   virtual Planar &XY( Spherical &ll ) = 0;
   virtual Spherical &LatLon( Planar &xy ) = 0;
};

The Projection class also defines the two methods for performing the actual projection: XY and LatLon. These methods are defined as pure virtual functions so that the Projection class will be abstract. This means of course that the only way to use this class is to derive a new class from it and provide code for these methods.

The two data members xyResult and 11Result are for returning values from these two methods.

The Mercator Projection

In 1569 Gerardus Mercator developed a new method for projecting a sphere (the world) onto a plane (a map) and used it to create a map of the entire known world. Mercator's projection is still one of the most popular methods used to this day, and has been standard since 1910 for nautical charts prepared by the U.S. government.

Mercator's projection is cylindrical. Picture an enormous cylinder wrapped around the earth and touching all the way around the equator (see Figure 1) . Each point on the earth is projected onto this cylinder by extending its radius vector. Then the cylinder is cut vertically and unrolled into a plane. This is the geometric description of a cylindrical projection.

The effect of this projection is that the meridians, or lines of longitude, become parallel. A significant characteristic of the Mercator projection is that lines of latitude and longitude meet at right angles and form a grid. At higher latitudes, the meridians are spread to make them parallel. So the lines of latitude have to be spread as well. As we approach the pole, the spreading of both the meridians and the lines of latitude approaches infinity.

This is why the Mercator projection is not used for high latitudes. The spreading effect makes land masses which are actually small appear to be huge. And it is impossible to show the poles, since you would need a map of infinite size to do it!

The Mercator Class

The Mercator class is derived from the Projection class. It introduces some private members that it needs to do its job:

Spherical MinLL, MaxLL;

These are the boundaries of the area to be dealt with. The defined scale of a Mercator projection is only precise or "true" along one line of latitude. This true-scale latitude is assumed to be midway between the maximum and minimum latitude. The maximum and minimum latitudes can be the same, in which case they define the true-scale latitude.

Points above the minimum latitude will have a positive Y value. Points below will be negative. Points on the right of the minimum longitude will have a positive X value, points on the left will be negative.

double Scale;

This value is computed to convert from radians to inches, using the defined ellipsoidal model and the scale of the chart to be read or displayed.

double Fmin;

This is the unscaled Y value of the minimum latitude, saved as a private member so that it need only be computed when the minimum latitude is defined.

double F( double lat );

This member function computes Y in radians, adjusted for the spreading effect of the cylindrical projection. This function takes a latitude in radians as its input.

void NewScale( double lat, double scale );

This member function computes the value for the data member Scale, which is a function of latitude and the chart scale at that latitude. The scale is specified in inches. So if the scale of the chart is 1:12,000 (one inch on the chart equals 12,000 inches on the ellipsoid) then scale would be 12000.0.

The Mercator class also declares several public members.

Mercator( int model,
   Spherical &min,
   Spherical &max,
   double scale );

This version of the constructor uses two points which define the chart area to compute the bounds and determine the true-scale latitude. This is useful when you are only given the bounds of the chart to work with. The scale is used along with the model identifier to compute the scale factor for the chart.

Mercator( int model,
   Spherical &origin,
   double scale );

This version of the constructor uses just one point, which defines the true-scale latitude and the minimum longitude. This is useful for displaying a chart which is not of a fixed size.

void SetBounds( Spherical &min,
   Spherical &max,
   double scale );

This method recomputes the scale and saves the new boundaries. This is used by the first constructor. It is also used when zooming in or out when displaying a chart.

void SetOrigin( Spherical &origin,
   double scale );

This method is used to change the origin point and the scale. It is used by the second constructor. It can also be used to pan when displaying a chart.

The XY and LatLon methods are defined by the Projection class, but made concrete by the Mercator class.

Calculations

If you don't really care how the math works, skip this part. I am only going to skim over it anyway. For a rigorous discussion of the math, refer to [1] and [3].

Longitude is easy. To go from longitude to X we simply take the difference of the longitude and the minimum longitude, convert to radians and multiply it by the scale factor:

xyResult[X] = RAD( ll[LON] -
   MinLL[LON] ) * Scale;

This is because the lines of longitude are being kept parallel and evenly spaced. To go the other way, we divide X by the scale to get radians, convert to degrees and then add to the minimum longitude:

llResult[LON] = MinLL[LON] +
   DEG( xy[X] / Scale );

Latitude is so much harder. This is because the latitude is being spread mathematically to account for spreading of the longitude (which is a function of latitude.) Got it? The workhorse of this operation is the function F( double lat ) which scales latitude exponentially.

double Mercator::F( double lat )
{
   double SinLat = sin( fabs(lat) );
   double eSinLat = e * SinLat;
   double Ym =
      0.5* log((1+SinLat)/(1-SinLat) *
      pow((1-eSinLat)/(1+eSinLat), e) );

   if( (lat < 0 && Ym > 0) ||
      (lat > 0 && Ym < 0 ) )
      Ym = -Ym;

   return Ym;
}

So to get a Y value from latitude, this function is called to get the distance in radians from the equator. The distance that was computed for the true-scale latitude (Fmin) is subtracted from this value and the result is scaled to inches:

xyResult[Y] = ( F( RAD( ll[LAT] ) ) - Fmin ) * Scale;

Even this is simple compared to converting Y back to latitude. For this, no direct computation is possible. We therefore must use successive approximations to arrive at our answer. First we convert Y into radians from the equator:

double Ym = xy[Y] / Scale + Fmin;
double t = exp( -Ym );   // for efficiency

For an initial guess at latitude we use:

    // old and new latitude values
                // in RADIANS
    double lat[2];
lat[1] = (PI / 2) - 2.0 * atan( t );

Now we use the old latitude to compute a more accurate answer:

lat[0] = lat[1];     // save last answer
double eSinLat = e * sin( lat[0] );
double f = pow( (1-eSinLat)/(1+eSinLat), e/2 );
lat[1] = (PI / 2) - 2.0 * atan( t * f );

We repeat until the error is acceptably small:

while( fabs( lat[1] - lat[0] > 0.1e-5 );

Finally, the result is converted from radians to degrees:

llResult[LAT] = DEG( lat[1] );

Applications

Let's say you want to create a geographic database of San Diego bay. First you would want to get the latest revision of the National Oceanic and Atmospheric Administration (NOAA) chart number 18773. Next, compile and link the provided code with MAIN=1 so that the test code is compiled. Now run the program and enter the chart scale: 12000, the true scale latitude: 32 42 (32 degrees, 42 minutes North), and the minimum longitude: -117 15 (117 degrees, 15 minutes West).

Now you can measure the distance in inches from the true scale latitude and the minimum longitude to any feature on the chart. Then use the "Convert X/Y to Lat/Lon" choice to find its latitude and longitude. Write down this answer for later entry into your database. Repeat these steps until you go insane.

To preserve your sanity, write a program to take input from a digitizer tablet on which you position the chart. Your program can then read digitizer inputs, convert them to inches, convert these to latitude and longitude, then store the points in your database automatically.

To display your data using the Mercator projection, put the origin of your chart, and the true-scale latitude, at the center of your display. Converting from pixels to inches using the DPI of your display, determine the X/Y location of the display boundaries. Convert these to latitude and longitude. Now read the points out of your database. If a point falls within the boundaries you just computed, convert it to inches, then convert it to pixels, and put it on the display. Keep in mind that points to the left of the center will have negative values for X in inches, and points below the center will have negative Y values in inches.

Geographic applications are becoming more and more numerous: Navigation, market research, surveying, et cetera. If your application needs to input or display geographic data, you may well find the techniques presented here useful.

References

[1] J. P. Snyder. Map Projections — A Working Manual, U.S. Geological Survey Professional Paper 1395 (U.S. Government Printing Office, Washington, D.C., 1987).

[2] N. Bowditch. The American Practical Navigator, Vols. 1 and 2 (Defense Mapping Agency/Topographic Center, Washington, D.C., 1975).

[3] P. D. Thomas. Conformal Projections in Geodesy and Cartography, U.S. Coast and Geodetic Survey Spec. Pub. 251, 1952.