MetaCRS - Coordinates & Projections

MetaCRS is a project started by Frank Warmerdam of PROJ.4 fame.

The initial basis for this project is to act as an anchor for an OSGeo Project encompassing several projections, and coordinate system related technologies.

The first technologies or libraries that are being used in this project are:

proj4js (Javascript: Rich Greenwood and Mike Adair )
proj.4
libproj4 (the projection-only library maintained by Gerald Evenden)
OSGSpatialReference (GDAL coordinate system translation classes)
CS-Map (the recently open sourced library from Norm Olsen/Autodesk)

So where do non-programmers fit into this?

Frank has covered that as well and there are areas where a little geodesist like myself can contribute. If you feel you can, by all means join the mailing list.

He states on the wiki the following suggestions:

• Common Spatial Reference System or Coordinate Reference System Names and Descriptions
• Coordinate System (and CRS related object) dictionaries. Stuff like the EPSG dictionary.
• Datum shift lists (towgs84), and datum grid shift files (NTv1, etc).
• Transformations, calculations, and algorithms written in pseudocode that can be edited in different languages.
• Descriptions of spatial reference systems that can be used by developers in different programming languages.
• Notes on transformation from different representations of a CRS (WKT, PROJ.4, GCTP, GML,...).
• Test suites with test points in a variety of coordinate systems and their lat/long and WGS84 equivelents).
• Articles on spatial reference systems and translations useful for programmers interested in spatial reference system implementations. For example: Understanding The Difference Between National Vertical Datum of 1929 and the North American Datum of 1988
  

This will be an interesting project to be involved in.

I also hope to explain in my blog a little later, some answers to some of the suggestions above.

The Robinson Projection: Not shaped like an Egg

A Pleasing View of the World

The Robinson Projection came into being in 1963 and was introduced by Dr. Arthur H. Robinson.

This projection can be classified as a pseudo-cylindrical projection because of its straight parallels, along each of which the meridians are spaced evenly. The central meridian is also a straight line and all other meridians are curved.

The projection is neither equal-area nor conformal, therefore abandoning both for a compromise for creating what Dr. Robinson felt produced a better overall view of the world. This was the first map projection to be developed for commercial interests. Rand McNally felt that many of the map projections in use did not present the earth as a whole very well. With Mercator the poles were distorted. Robinson, was essentially contracted to develop a map projection that did not maintain angle, direction, or limit distortion, but was sanctioned to produce a map projection that "looked good" for books and atlases.

Remember maps are designed for one of the following 4 reasons:

• Conformality - the shapes of places are accurate
• Distance - measured distances are accurate
• Area/Equivalence - the areas represented on the map are proportional to their area on the earth
• Direction - angles of direction are portrayed accurately

Dr. Robinson specified the projection to be constructed by referring to a table of cartesian coordinate values at specific intersections of latitude and longitude. The intermediate locations are to be found by interpolation. Dr. Robinson developed the projection through a series of trials, continually iterating till he settled upon the meridian shapes and parallel spacing most pleasing to the eye. In comparison with other Map Projections, they are mainly developed or are formulated as mathematical equations.

Parallels are straight parallel lines, equally spaced between latitudes 38 degrees north and south. Space decreases beyond these limits. The Equator is 0.8487 times as long as the circumference of a sphere of equal area. The central meridian is a straight line 0.5072 as long as the Equator. Other meridians are equally spaced elliptical arcs and concave toward the central meridian. The scale is true along latitudes 38 degrees north and south, constant along any given latitude, and the same for the latitude of opposite sign (Robinson 1974; Snyder and Voxland 1989).

This map projection is also based on a sphere, not an ellipsoid. This is an important point to remember, as the earth is being modelled differently.

What is the true shape of the Earth?

The Earth, in actual fact, is shaped more like an egg, hence, even an ellipsoid is not the best model. In the past, when datum's were more local (such as NAD27, SAD69, Cape Datum, etc.), various ellipsoids were tied to local points on Earth. For NAD27, this was Meades Ranch, Kansas.

But back to the Earth's shape; it is very close to an oblate spheroid — a rounded shape with a bulge around the equator — although the precise shape (the geoid) varies from this by up to 100 metres.

The rotation of the Earth creates the equatorial bulge so that the equatorial diameter is 43 km larger than the pole to pole diameter.

What about height? A whole other story.

Interesting, eh? There are so many different ways to see the Earth. When we look at height systems and describe the geoid, we are describing an equipotential surface. But height is a whole other story, as gravity is involved and it is not as simple as getting the "Zed" or "Zee" measurement from your GPS.

I'll explain height later and how we can determine MSL (and where the "mean" actually comes from).

Google Maps: Is the Earth a Sphere or Ellipsoid?

Google Maps sees the Earth as a Spheroid, not an Ellipsoid. This came up through a discussion on the PROJ mailing list and I thought it was interesting to point out how Open Source can even handle projected lat/long systems (such as Google Maps) using a very familiar tool called cs2cs.

Christoper Schmidt wrote about it on his blog and also on his blog he points out the EPSG code to use. The magical number for the Google Mercator Projection (of a lat/long grid based on a sphere) is: 900913

Now onto the fun, showing how we can use Open Source to have our data show within a KML project and Google Maps. Quoting from Frank's FAQ, he provides an excellent example, we see the following the use of cs2cs:

"cs2cs +proj=latlong +datum=WGS84 +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +no_defs"

Notice the sphere is being used? a=b

Because we are dealing with a sphere, the Y values will be greatly different from those on an ellipsoid (30 to 100 metres or more).

Quoting Frank again:

"In this case, and many other cases using spherical projections, the desired approach is to actually treat the lat/long locations on the sphere as if they were on WGS84 without any adjustments when using them for converting to other coordinate systems. The solution is to "trick" PROJ.4 into applying no change to the lat/long values when going to (and through) WGS84. This can be accomplished by asking PROJ to use a null grid shift file for switching from your spherical lat/long coordinates to WGS84.

cs2cs +proj=latlong +datum=WGS84 +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

Note the strategic addition of +nadgrids=@null to the spherical projection definition"

As you can see the value of Open Source and the mailing list and Open Source software is that people are sharing knowledge - whether it be via blogs, lists, or some other means of communication. There is a community out there that supports each other. These are actual users facing everyday problems and looking for solutions. The answers do exist, just the question has to be asked, and the community comes together to help.

Shapefiles and PRJ - Tying them together

ESRI has developed a de-facto standard for Shapefiles, but sometimes we as users' of Shapefiles forget something, a Shapefile is actually a minimum of three files.

The three mandatory/required files are:

• .shp - the file that stores the feature geometry.
• .shx - the file that stores the index of the feature geometry.
• .dbf - the dBASE file that stores the attribute information of features.

A recent thread on the GDAL Mailing List led to the first link about Shapefiles. The key to the thread was a discussion about where Projections are stored in the Shapefile definition. As it can be seen, there is no location in the required files for projections.

Therefore ESRI adopted a new file type (PRJ) with the extension .prj.

So how are projections defined in this file?

As we know coordinate systems in terms of mapping can either geographic (longitude, latitude) or projected (X, Y). In the PRJ definition, the coordinate system is composed of several objects, with every object having a keyword in uppercase. Objects can be composed of other objects. ESRI calls the string in this file a Projection Engine (PE). The Sole purpose of the Projection Engine is to store the metadata for a coordinate system in a string, or in a .prj file. This string, which ESRI also calls a PE (not for Physical Education!) string, must be continuous and not broken.

Now the scary part: You can define your own units, datums, and spheroids!

An example, taken from the ESRI website, shows how we can define in .prj file a projected or geographic definition.

As we know projected coordinate systems (of which maps are) are based upon a geographic coordinate system (latitude, longitude), so the in their sample file, a projected coordinate system first is defined.

For example, UTM zone 10N on the NAD83 datum is defined as

PROJCS["NAD_1983_UTM_Zone_10N",
,
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",500000.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",-123.0],
PARAMETER["Scale_Factor",0.9996],
PARAMETER["Latitude_of_Origin",0.0],
UNIT["Meter",1.0]]

The geographic coordinate system name is followed by the datum, the prime meridian, and the angular unit of measure.

The geographic coordinate system string for UTM zone 10N on NAD 1983 is:

GEOGCS["GCS_North_American_1983",
DATUM["D_North_American_1983",
SPHEROID["GRS_1980",6378137,298.257222101]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.0174532925199433]]

The full string representation of NAD 1983 UTM zone 10N is:

PROJCS["NAD_1983_UTM_Zone_10N",
GEOGCS["GCS_North_American_1983",
DATUM["D_North_American_1983",
SPHEROID["GRS_1980",6378137,298.257222101]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",500000.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",-123.0],
PARAMETER["Scale_Factor",0.9996],
PARAMETER["Latitude_of_Origin",0.0],
UNIT["Meter",1.0]]

As I stated earlier, you can define your own to use the predefined names for map projection and parameter object's.

I hope the above information aids in your understanding of Shapefiles and how projections and coordinate systems are defined.

My best advice for dealing with PRJ files; copy one that you already have and use it as a base, then modify the parameters as you need.

Good luck and have fun!