Monthly Archives: July 2011

Peony in rain

Projections for Programmers – One Projection (Part 1)

Unfamiliar with the terms in this post?   If so, refer to my other post that defines geospatial terms.

If you’re a programmer new to GIS,  one of the first things you’ll bump into is map projections.   Understanding map projections is a significant part of the GIS domain of expertise, and in this post I’ll do my best to explain the parts of a projection so that you’re not scratching your head the next time you look at one.

A good way to look at a the various pieces that make up a projection is with a textual representation of the projection parameters.   There are a couple of ways to to use text to describe a projection – PROJ4 parameters, an ESRI prj file, a TIFF world file, and others. The OGC has an open Well Known Text (WKT) format for describing projections that’s very similar to an ESRI PRJ file.

Here’s the OGC WKT for UTM Zone 14N, downloaded from http://spatial-reference.org/.

PROJCS["NAD83 / UTM zone 14N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]
            ],
            AUTHORITY["EPSG","6269"]
        ],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]
        ],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]
        ],
        AUTHORITY["EPSG","4269"]
    ],
    UNIT["metre",1,
       AUTHORITY["EPSG","9001"]
    ],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-99],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","26914"],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]
]

What is all this?  At a glance, there’s a lot to understand, but it’s not so bad if you break it down.

First of all, let’s get that keyword AUTHORITY out of the way. It’s good to understand what it’s about, but not really pertinent to the understanding of a projection.  All AUTHORITY tells you is who came up with the definition for that part of the projection (in this example, everything is defined by the EPSG – the European Petroleum Survey Group),   and where you can look in the authority’s database to get more information.  The most useful authority setting is the value 26914 on line 26, which is the EPSG code for this projection.   If you want to look up these projection details yourself,  you can just plug EPSG:26914 into the search field at http://www.spatialreference.org.   If you’re using something like GeoTools or GDAL for your geospatial coding,  you can easily refer to this projection with the EPSG code.

Enough about AUTHORITY.   Let’s grey out the AUTHORITY definitions in the WKT so they don’t distract from the real content.

PROJCS["NAD83 / UTM zone 14N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]
            ],
            AUTHORITY["EPSG","6269"]
        ],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]
        ],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]
        ],
        AUTHORITY["EPSG","4269"]
    ],
    UNIT["metre",1,
       AUTHORITY["EPSG","9001"]
    ],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-99],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","26914"],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]
]

With the AUTHORITY definitions less visible, let’s start on the outside and work our way in.

PROJCS["NAD83 / UTM zone 14N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]
            ],
            AUTHORITY["EPSG","6269"]
        ],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]
        ],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]
        ],
        AUTHORITY["EPSG","4269"]
    ],
    UNIT["metre",1,
       AUTHORITY["EPSG","9001"]
    ],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-99],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","26914"],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]
]

The outer element is the PROJCS – the “Projected Coordinate System”.  It’s just a wrapper for everything that defines the projection.  There’s the name “NAD83 / UTM zone 14N” that tells you a little something about the projection, but don’t take that name as some canonical reference.  Each authority usually has a slightly different name for the exact same projection.

Simple enough – now let’s move on to the next element in the WKT, the geographic coordinate system.

PROJCS["NAD83 / UTM zone 14N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]
            ],
            AUTHORITY["EPSG","6269"]
        ],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]
        ],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]
        ],
        AUTHORITY["EPSG","4269"]
    ],
    UNIT["metre",1,
       AUTHORITY["EPSG","9001"]
    ],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-99],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","26914"],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]
]

Every projection has a curved surface from which it was projected that’s called the geographic coordinate system.    In the WKT above, this is defined by the GEOGCS attribute.   Like the PROJCS, the GEOGCS has a name – NAD83 – that tells you a little something about it, and like the PROJCS this shouldn’t be considered the canonical name for the coordinate system, since other authorities may have a slightly different name for it.   If we ignore the name, the GEOGCS is really defined by the SPHEROID, PRIMEM (prime meridian) and UNIT.

The SPHEROID element defines the shape of the earth used for this particular projection.  As in most cases, it’s an ellipsoid, and in this particular case it’s the ellipsoid defined by “Geodetic Reference System 1980”.  The ellipse is defined by its semi-major axis of  6378137 and its inverse flattening of 298.257222101.

The PRIMEM tells you where the Y-axis is (at longitude 0, in Greenwich, England).   The X-axis isn’t given, and is assumed to be at the equator (latitude 0).

And UNIT tells you that the GEOGCS has its units defined in degrees.  Since the standard angular measure of unit in EPSG is the radian, the UNIT element also provides a conversion factor for converting a degree to a radian (multiply by 0.01745329251994328, or 2π/360)

Wow!   That’s a lot of info for one post.   Let’s take a break and cover the other projection elements in a second post.