Sponsored Links

Polar Maps and Projections: Part 2, Implementation PDF Print E-mail
Written by Richard Marsden   
Monday, 18 January 2010 09:27

The first part of this article looked at different ways of producing polar maps and surveyed a number of different azimuthal projections that are often used for polar maps.

In this second part, I produce a working implementation using UMN MapServer and OpenLayers. The working implementations can be found in the Polar Map section of Equal-Area-Maps.com .

 

This implementation is based on my previous article that used MapServer to implement non-Mercator projections for Equal-Area-Maps.com. This article can be found here.

This example uses the Stereographic Projection centered on the north pole, but it is easily modified to use the south pole, or to use other azimuthal projections that were covered in the first part of this article. Oblique projections are possible, but these will require more work when it comes to map data preparation and clipping.

Converting the equal-area code to produce a polar map is actually quite simple. The biggest problem involves data clipping. Previously we were interested in global maps. Polar projections are often only valid for one hemisphere. Even those that can project both projections (eg. the Stereographic Projection) produce seriously distorted antipodal hemispheres. Therefore, we need to clip our map and map data to the hemisphere of interest.

As before, we use MapServer to produce a base map consisting of coastlines, national boundaries, and a graticule. The graticule is also used to create an "ocean background" in contrast with the surrounding background of the map.

Previously, the graticule was created using a Python script which called various shapelib utilities. This script created a series of 30x30 degree 'squares' which tiled the entire globe. This script is readily modified to only tile the northern hemisphere. Here is the modified script:

ERROR [include_code_listing plugin]: File Not Found (/usr/www/users/winwaed/geowebguru/img/2010/north_graticule.py)

 

Only two changes have been made to the script. First, the loop which controls the latitude of the squares is modified so that it only loops from the equator to the north pole. Secondly, the line segments on each square are produced with more data points. This is required because the azimuthal projections tend to result in meridians and parallels that have sharper curves than the cylindrical and pseudo-cylindrical projections. More data points result in a smoother, less-jagged final rendering.

The lakes shape file is very easy to modify: All of our lakes (the Great Lakes and the Caspian Sea) are in the northern hemisphere, so they do not need to be clipped.

Unfortunately our world borders shape files cross the equator in a number of places. These could be clipped using a script that read the shape file and processed it a coordinate at a time. This is probably the most efficient way if you had to clip large numbers of files to one hemisphere or another. I only had two shape files (world_borders.shp and world_borders_simple.shp), so I chose to use existing tools.

The GDAL library contains the very useful OGR2OGR utility. This is intended to convert between different vector formats, but it has a growing number of transformation operations. The current development version, GDAL 1.7, contains a new set of clipping operations. The new -clipsrc option will clip a shapefile to a rectangle defined in geographic coordinates. Therefore, the following command clips our world_borders.shp file to the northern hemisphere:

ogr2ogr -f "ESRI Shapefile" world_borders_north.shp world_borders.shp -clipsrc -180.0 0.0 180.0 90.0

This is a quick way of clipping shapes, although the equatorial line segments introduced by the clipping are simple line segments. OGR2OGR does not insert extra points. As with the graticule, we have to insert a number of data points along these equatorial line segments in order to produce a smooth line in the final projection. This is where the custom programming solution would come in to its own. Only a few shapes need to be modified (Brazil and the Democratic Republic of Congo being the largest), so I chose to manually add them using Quantum GIS. Quantum GIS is an open source GIS system. This is definitely over-kill for adding a few data points, but it was a good excuse to familiarize myself with the system.

The shapes have all been clipped and modified for the northern hemisphere. We are now ready to create the MapServer MAP file. This is almost identical to the Equal-Area-Maps.com MAP files, except the we use the new northern hemisphere shapefiles, and we use different projection parameters. Here is the MAP file for the Stereographic Projection:

ERROR [include_code_listing plugin]: File Not Found (/usr/www/users/winwaed/geowebguru/img/2010/stere.map)

 

Similarly, the web page and client JavaScript are virtually unchanged. New Proj4JS projection definitions are required for the azimuthal projections:

ERROR [include_code_listing plugin]: File Not Found (/usr/www/users/winwaed/geowebguru/img/2010/stereographic.html)

 

Proj4JS supports most of the azimuthal projections discussed in Part 1. The Near-Sided Perspective Projection (+proj=nsper) is not supported. The Gnomonic (+proj=gnom) Projection was not supported in Proj4JS when the azimuthal maps were being implemented, but I wrote and submitted my own implementation. This is now listed as an untested projection supported Proj4JS.

In these examples, Proj4JS is used to re-project KML and GeoRSS overlays. Due to web browser security restrictions, these files must be from the same domain as the host webpage.

Another problem with these overlays is that they are not clipped to the northern hemisphere. Some projections (eg. the Gnomonic) will not plot antipodal hemisphere data points, but others (eg. the Stereographic) will. This looks unsightly, but there are currently no easy fixes. One ugly approach would be to produce a special version of Proj4JS which projects the out-of-range points to infinity (or close enough). A better approach would be to override an OpenLayers class to do the clipping. I shall revisit this issue in a future article.

Conclusions

And there we have it. The existing equal area maps were easily modified to use new polar azimuthal projections. The only area that we have to be careful with, is that we have to clip the data to the area of interest, and exclude data from the opposing hemisphere.

Although only one of these projections is an equal area projection, I have posted the resulting maps in a new Polar Projection section on the Equal Area Maps website.