This section describes some of the inner workings of the coordinate transforms used in the map objects. Here is a diagram showing the modules and their sources:
This Python module provides a number of functions and classes for representing and manipulating terrestrial coordinates. Functions include:
Given an angle in degrees, this function returns the corresponding angle in radians.
Given an angle in radians, this function returns that angle as degrees.
Given a distance on the earth's surface in feet, this function returns the angle subtended by that distance from the center of the earth, in radians.
The inverse function of feetToAngle.
Given an angle in degrees, this function returns the equivalent angle as a tuple (d,m,s) where d is whole degrees, m is whole minutes, and s is seconds as a float.
Inverse function of degDMS; the argument can be either a 2-tuple (d,m) or a 3-tuple (d,m,s), and each value can be either int or float type.
There are three classes in module terrapos.py:
This class represents an arbitrary position on the earth's surface as a latitude and longitude in radians. It may optionally include an elevation above sea level in feet and a map datum value (such as WGS-84 and NAD-27) defining the overall coordinate system in use. (Most GPS coordinates are given as WGS-84, while most topo maps use NAD-27. The difference can be as much as 200 feet.)
Objects of this class support a number of operations in spherical geometry, such as the ability to find the surface (rhumbline) distance between two points, or to find a new position given a starting position, bearing, and distance.
This class is derived from the TerraPosition class. It differs from the base class in that it explicitly represents lat-long coordinates in degrees (as opposed to UTM coordinates).
You might ask, why doesn't the TerraPosition class suffice for representing lat-lon coordinates, since we can easily convert the radian values to degrees? Well, we certainly could do it that way, but there is an annoying little pathology of such conversions that the LatLon class is designed to circumvent.
Suppose we take an angle, such as 107° 50' 0" and convert it to radians. For some values, converting back to mixed units, and then rounding the individual units, can cause this value to be displayed as 107° 49' 60.0".
To prevent this kind of ugliness, the LatLon class remembers the original units and regurgitates those units when displaying the value in string form, rather than converting them from radians. The class constructor also takes care of converting the original values to radians for computational convenience in the base class.
An instance of class TerraBox represents a rectangular area in lat-long coordinates. This is useful for applications such as producing a map of a specific area. The class also supports methods that can tell you whether a point is inside such an area; whether two areas overlap; and the rectangles defining the intersection and the union of two rectangles.
Module mapbase.py contains mechanisms for building up a composite image from individual tiles and supporting transformations between image coordinates and terrestrial coordinates.
Principle classes include:
Represents the entire image pyramid: all sizes and kinds of tiles.
Represents all the available map information for a specific rectangular area, that is, for a range of latitudes and longitudes. Such an area will appear as a rectangle in the usual Mercator projection with latitude running vertically and longitude running horizontally.
A MagBox object manages the geometry of assembling a set of tiles into a single rectangle, and can produce either photo or topo images.
An object of this class represents the corners file for a given magnification. Its purpose is to manage just the geometry of tile slots, the square holes into which tiles fit. It is a container for TileCorner and TileSlot objects.
A TileSlot object describes one tile slot: its corners, and its tile number, which is called tileCR (for tile row and column) internally.
A tile slot is not considered to exist unless all four of its corners are known.
Each TileCorner object describes one corner of a tile, corresponding to one line of the corners file for a given magnification.
An object of this class represents one UTM coordinate. In the current version, there is no obvious way to convert between UTM and lat-long coordinates. When such algorithms become available, they will be moved into the terrapos.py module.
Working with tiles requires that we convert back and forth between two coordinate systems:
Display coordinates are the coordinates of pixels on the screen. As usual for windows, y increases from top to bottom, an inversion of the usual Cartesian y coordinate.
Each map tile has an upper left (northwest) corner at display coordinate (0,0) and occupies a 200 x 200 pixel area on the screen. So, for example, the pixel in the southwest corner is at display coordinate (0,199).
Terrestrial coordinates are latitude and longitude (although UTM coordinates may be supported in the future).
For assorted complex historical reasons, longitude and latitude lines are not exactly parallel to the x and y coordinates of display space. Accordingly, the algorithm for converting between display coordinates and terrestrial coordinates is not trivial, but it is reasonably straightforward.
Here is a diagram of an idealized map tile:
In this figure, the four corners are labeled NW, NE, SW, and SE. Consider a point P within the tile. Point W is the point on the west edge of the tile with the same latitude as point P, and point E is the point on the east edge with the same latitude as P. Assuming that the tile is an accurate Mercator projection of the terrain, there should be a straight line EPW connecting all points in the tile with the same latitude. Similarly, line NPS connects the points with the same longitude.
Therefore, deriving the display coordinate of P from its lat-long coordinates is straightforward:
Find the coordinates of point W by interpolating the latitude of P between the latitudes of points NW and SW, and then project that fraction onto the display length of 200.
Find the coordinates of point E by interpolating along the east side similarly.
Find the coordinates of points N and S by interpolating P's longitude along the north and south edges respectively.
Derive the coordinates of P by finding the intersection of lines WE and NS. We use the two-point form of the equation of a line, (y-ya)/(x-xa)=(yb-ya)/(xb-xa), where the coordinates of the two points are (xa,ya) and (xb,yb). This gives us equations for lines WE and NS, and it is simple linear algebra to solve the system of these two equations in two variables.
The inverse transform uses the same procedure, simply exchanging the two coordinate systems.