Calibrating a map

It is important understand that there is more than one way to draw a map of any area. If you look in an Atlas you may see some maps of the world which are rectangular, with the poles stretched out so they are the same size as the equator. You may also see maps of the world with the surface layed out like segments of orange peel. These are just two examples of the many different map projections used by cartographers.

You need to remember that an algorithm which correctly converts map coordinates to actual locations in projection X will not work for maps created using a different projection. There's no "one size fits all" method of doing this.

As well as projection, you should also be aware of the concept of ellipsoid. The world is not perfectly spherical. It is slightly squashed and its shape is not precisely regular. To a cartographer, an ellipsoid is a mathematical model of the earth's shape and again there are several different ones to choose from. The reason for the differences is not just that cartographers cannot agree. Because of the irregularities of the earth's shape, an ellipsoid which works well for one country or region may not work so well for another.

The differences between ellipsoids are slight and you may be able to ignore them if you don't need to be too precise, but ellipsoids are mentioned in the following notes so you need to understand what they are.

Most maps of small areas use projections which have the following properties:

Distances between objects on the map are directly proportional to the distances between those objects in the real world. For example if your map is described as "1 inch to the mile" then it is reasonable to assume that an inch in any direction represents a mile in the real world. This is not true of all maps. Going back to our rectangular world map in the atlas, a horizontal inch on such a map would represent hundreds of miles in the real world whilst a horizontal inch near the pole might only represent a mile or two.
The map is (approximately) "orthogonal". By that I mean that if you take any arbitrary point on the map then a line pointing due North from that point would be at right angles to a line pointing due East. Be aware that this is an approximation, which is why I am assuming your map is of a small area. The North pole is a point so, in actual fact, a line from the bottom left corner of your map, pointing due North, would probably lean slightly to the right whilst a line from the bottom right will lean slightly left. If your map is of a small area, the angles will be so slightly off 90 degrees that we can effectively ignore the difference.

In my own mapping software, I assume that all maps have the above properties. You just need to be aware that you are making assumptions like these when devising your calibration algorithm.

If you want to convert to latitude and longitude, you have another problem to overcome. The distance between lines of longitude is not constant. The further you go from the equator, the closer they get - eventually meeting at the pole. The implications of this are twofold. Firstly, the scale multiplier for converting a map distance in the EW direction into degrees of longitude will be different from the multiplier for distances in the NS direction. Secondly, the EW multiplier will be different at the bottom of your map from the top (assuming top is North). The difference in EW multipliers from top to bottom of map can probably be ignored if the map is of a small area and not too near either pole. The difference between the NS and EW multipliers cannot be ignored.

In MapMan, I invite the user to supply two reference points for which both the map coordinates and real world coordinates are known. The first thing I do is to calculate an aspect ratio - that is the ratio between the multipliers to convert EW map distances into degrees East and the multiplier to convert NS map distances into degrees North. Note that I can't yet calculate the actual multipliers but I can calculate the ratio between them based upon the average latitude of the two reference points. Here's a code snippet for calculating this ratio:

double NSRadius = EarthRadiusAtThisLatitude();
double EWRadius = NSRadius * cos(DegreesToRadians(Northings));
// If the EW radius is close to zero then we have a problem -
// we are too near the pole.
if (EWRadius < 100)        Report an error

return NSRadius / EWRadius;

The calculation EarthRadiusAtThisLatitude() is where the concept of ellipsoid rears its ugly head. To calculate this value, you need to know the parameters of the ellipsoid you are using. Here's another snippet:

double Ellipsoid::EarthRadiusAtLatitude(double DegreesNorth) const
{
    double cosLat = cos(DegreesToRadians(DegreesNorth));
    return sqrt( (EquatorialRadiusKm * EquatorialRadiusKm * (1 - e2))
                    / (1 - e2 * cosLat * cosLat));
}

The value e2 is the eccentricity squared. This depends upon the ellipsoid used to create the map you are working from. For the GRS80 ellipsoid the value is approximately 0.00669. The equatorial radius in km also depends upon the ellipsoid. For GRS80 it is 6378.137

Let's now consider the two lines, one on the map and one in the real world, running from reference point 1 to reference point 2. These lines can be defined as vectors x2-x1, y2-y1. Where x1,y1 are the coordinates of ref pt 1 and x2,y2 are the coordinates of ref pt 2

Next I work out the bearing (angle) of the line from ref pt 1 to 2 both on the map and in the real world. In both cases, the angle can be calculated as the arctangent of horizontal and vertical distances but in the case of the real world vector you need to adjust the eastings distance by the aspect ratio calculated above. The difference between the map bearing and the real world bearing is the amount of rotation (if any) you need to apply to your map to orient it so that North is directly up the page.

If you know your map is correctly oriented to start with, then you can skip this calculation but if your map is, say, a scan of a paper map then you need to calculate the rotation because you cannot guarantee you put the paper down exactly square on the scanner.

Next I rotate reference point 2 (using reference point 1 as the centre of rotation) so that the map bearing matches the real world bearing. Here's the rotation code:

RotatedMapDistance.dx = MapDistance.dx * CosR + MapDistance.dy * SinR;
RotatedMapDistance.dy = MapDistance.dy * CosR - MapDistance.dx * SinR;

Now I can work out the scale factor needed to convert map distances to real world coordinates. Remember that this is two numbers, not one, because our target units are degrees and degrees are different sizes depending on whether they are Northings or Eastings. The scale calculation is simply:

Scale.Eastings = RealWorldDistance.dx / RotatedMapDistance.dx;
Scale.Northings = RealWorldDistance.dy / RotatedMapDistance.dy;

Finally I am in a position to calculate my transformation matrix. This is a matrix by which you can multiply Screen coordinates to convert them into real world coordinates taking into account the rotation and scaling. At the same time I work out the inverse matrix which does conversion in the other direction

m_TransformationMatrix[0][0] = Scale.Eastings * CosR;
m_TransformationMatrix[0][1] = Scale.Eastings * SinR;
m_TransformationMatrix[1][0] = 0 - Scale.Northings * SinR;
m_TransformationMatrix[1][1] = Scale.Northings * CosR;

// Calculate the inverse matrix
double Determinant = m_TransformationMatrix[0][0] * m_TransformationMatrix[1][1]
                     - m_TransformationMatrix[1][0] * m_TransformationMatrix[0][1];

m_InverseMatrix[0][0] = m_TransformationMatrix[1][1] / Determinant;
m_InverseMatrix[0][1] = 0 - m_TransformationMatrix[0][1] / Determinant;
m_InverseMatrix[1][0] = 0 - m_TransformationMatrix[1][0] / Determinant;
m_InverseMatrix[1][1] = m_TransformationMatrix[0][0] / Determinant;