Polar Maps and Projections: Part 2, Implementation 2

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.

 

This implementation is based on my previous article that used MapServer to implement non-Mercator projections. 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 30×30 degree ‘squares’ which tiled the entire globe. This script is readily modified to only tile the northern hemisphere. Here is the modified script:

#!/usr/local/bin/python

# Shapefile Graticule Creation Script for the Northern Hemisphere
# (C) Copyright 2009 Winwaed Software Technology LLC
# (Berkeley License applies if not stated explicitly)
#
# Based on graticule.py (whole world graticule), this version 
# only creates a graticule for the northern hemisphere.
# This is intended for polar projections (lat=90).
# The script can be easily modified for the southern hemisphere.
#
# usage: north_graticule.py file 30
# - Creates a shapefile with a 30deg graticule
# Shapefile is polygon with rectangular graticule shapes covering the 
# entire northern hemisphere (lat=0->90)
# Filename should NOT include the .shp extension

import sys
import os
import string


fname = sys.argv[1]
angular_size = string.atoi( sys.argv[2] )

print "Creating Graticule (N.Hemisphere only), size: %d" % (angular_size)

# Create the shapefile and a matching dbf file

os.system("./shpcreate "+fname+".shp polygon")
os.system("./dbfcreate "+fname+".dbf -s NAME 6")

half_ny = 90 / angular_size
half_nx = 180 / angular_size

# Grid is defined to a higher resolution for a polar projection compared to 
# the cylindrical and pseudo-cylindrical projections of before
sub_ang = angular_size / 19.0
sub_list = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19];

ix = -half_nx
while ix < half_nx:
   iy = 0  # start from the equator and iterate to the North Pole

   while iy < half_ny:
      # Create the new shape coords in a string
      cmd_string = "./shpadd " + fname + ".shp"
      basex = ix * angular_size
      basey = iy * angular_size
   
      cmd_string = cmd_string + " %f %f" % ( basex, basey );

      for i in sub_list:
        cmd_string = cmd_string + " %f %f" % ( (basex), (basey+sub_ang*i) );
      for i in sub_list:
        cmd_string = cmd_string + " %f %f" % ( (basex+sub_ang*i),  (basey+angular_size) );
      for i in sub_list:
        cmd_string = cmd_string + " %f %f" % ( (basex+angular_size),  (basey+angular_size - i*sub_ang) );
      for i in sub_list:
        cmd_string = cmd_string + " %f %f" % ( basex+angular_size-i*sub_ang,  basey);

      # Update the shapefile with this shape
      os.system( cmd_string );
      os.system( "./dbfadd "+fname+".dbf \"GRID\"" );

      iy = iy + 1

   ix = ix + 1 

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:

# POLAR Azimuthal Stereographic (stere) Projection Base Map
MAP
NAME "outline"
UNITS meters
EXTENT -13000000 -13000000 13000000 13000000
SIZE 640 480
IMAGECOLOR 255 255 255
#TRANSPARENT on
IMAGETYPE png
SHAPEPATH "/path_to_your_shapefiles/"
FONTSET "/your_font_directory/fontset.txt"
CONFIG "PROJ_LIB" "/your_local_directory/share/proj/"
DEBUG on
STATUS ON

##################################### 
# Web object
#
WEB
     IMAGEPATH "/your_web_dir/tmp/"
     IMAGEURL "/tmp/"
     METADATA
       "wms_title"     "WMS Azimuthal Map Server"
       "wms_onlineresource" "http://www.yourdomain.com/mapserv.cgi?map=/your_map_directory/stere.map"
       "wms_srs"       "EPSG:3411"
     END
     LOG "/your_local_directory/mapserv.log"
END

PROJECTION
  "proj=stere"
  "lat_0=90"
  "lat_ts=70"
  "lon_0=0"
  "k=1"
  "units=m"
END


##################################### 
# Ocean solid fill

LAYER
     NAME "outline"
     METADATA
       "wms_title"   "outline"
     END
     DATA "grid30_north"
     PROJECTION
        "init=epsg:4326"
     END
     STATUS default
     TYPE polygon
     CLASS
          STYLE
               COLOR 192 192 255
          END
     END
END


##################################### 
# Coloured Fill: Land, Zoomed out
#
LAYER
     NAME "outline"
     METADATA
       "wms_title"   "outline"
     END
     DATA "world_borders_simple_north"
     PROJECTION
        "init=epsg:4326"
     END
     MINSCALE 20000001
     STATUS default
     TYPE polygon
     CLASS
          STYLE
               COLOR 241 238 232
          END
     END
END

##################################### 
# Coloured Fill: Land, Zoomed in
#
LAYER
     NAME "outline"
     METADATA
       "wms_title"   "outline"
     END
     DATA "world_borders_north"
     PROJECTION
        "init=epsg:4326"
     END
     MAXSCALE 20000000
     STATUS default
     TYPE polygon
     CLASS
          STYLE
               COLOR 241 238 232
          END
     END
END

##################################### 
# Coloured Fill: Lakes and Inland Seas
#
LAYER
     NAME "outline"
     METADATA
       "wms_title"   "outline"
     END
     DATA "lakes_simple"
     PROJECTION
        "init=epsg:4326"
     END
     STATUS default
     MINSCALE 20000001
     TYPE polygon
     CLASS
          STYLE
               COLOR 192 192 255
          END
     END
END
LAYER
     NAME "outline"
     METADATA
       "wms_title"   "outline"
     END
     DATA "lakes"
     PROJECTION
        "init=epsg:4326"
     END
     STATUS default
     MAXSCALE 20000000
     TYPE polygon
     CLASS
          STYLE
               COLOR 192 192 255
          END
     END
END


##################################### 
# Shape boundaries to be drawn (black)
#
LAYER
     NAME "outline"
     METADATA
       "wms_title"   "outline"
     END
     DATA "world_borders_simple_north"
     PROJECTION
        "init=epsg:4326"
     END
     STATUS default
     MINSCALE 20000001
     TYPE line
     CLASS
          STYLE
               SIZE 1
               COLOR 0 0 0
          END
     END
END
LAYER
     NAME "outline"
     METADATA
       "wms_title"   "outline"
     END
     DATA "world_borders_north"
     PROJECTION
        "init=epsg:4326"
     END
     STATUS default
     MAXSCALE 20000000
     TYPE line
     CLASS
          STYLE
               SIZE 1
               COLOR 0 0 0
          END
     END
END


##################################### 
# Graticule

LAYER
     NAME "outline"
     METADATA
       "wms_title"   "outline"
     END
     DATA "grid30_north"
     PROJECTION
        "init=epsg:4326"
     END
     STATUS default
     TYPE line
     CLASS
          STYLE
               COLOR 95 95 95
          END
     END
END


END # mapfile

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

<html>
<head>
<title>Stereographic Polar Map Projection</title>

<link rel="stylesheet" href="http://www.equal-area-maps.com/theme/default/style.css" type="text/css" />
<link rel="stylesheet" href="http://www.equal-area-maps.com/theme/default/framedCloud.css" type="text/css" />
<link rel="stylesheet" href="http://www.equal-area-maps.com/theme/default/scalebar-fat.css" type="text/css" />

<style type="text/css">
    #map {
        width: 100%;
        height: 512px;
        border: 1px solid black;
    }
    .olPopup p { margin:0px; font-size: .9em;}
    .olPopup h2 { font-size:1.2em; }
</style>

<script src="http://www.equal-area-maps.com/proj4js/lib/proj4js.js"> </script>
<script src="http://www.equal-area-maps.com/OpenLayers.js"></script>
<script src="http://www.equal-area-maps.com/control/ScaleBar.js"></script>

<script type="text/javascript">

// Projection definitions for Proj4JS

// Standard geographic coords
Proj4js.defs["EPSG:4326"] = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";

// Lambert Azimuthal Equal Area
//Proj4js.defs["EPSG:3408"] = "+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
Proj4js.defs["EPSG:3408"] = "+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +units=m +no_defs";
// Stereographic
Proj4js.defs["EPSG:3411"] = "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
// UPS (special variant of the Stereographic)
Proj4js.defs["EPSG:32661"] = "+proj=stere +lat_0=90 +lat_ts=90 +lon_0=0 +k=0.994 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";

/////////////////////////
// Our own invented codes
// Orthographic
Proj4js.defs["EPSG:102015"] = "+proj=ortho +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
// Azimuthal Equidistant
Proj4js.defs["EPSG:102016"] = "+proj=aeqd +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
// Near-Sided Perspective (DUMMY)
Proj4js.defs["EPSG:102017"] = "+proj=ortho +lat_0=90 +lon_0=0 +h=7000000 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
// Gnomonic
Proj4js.defs["EPSG:102018"] = "+proj=gnom +lat_0=90 +lon_0=0 +x_0=6300000 +y_0=6300000 +ellps=WGS84 +datum=WGS84 +units=m +no_defs";



var icon_list = new Array();
icon_list[0] = 'http://www.equal-area-maps.com/img/marker.png';
icon_list[1] = 'http://www.equal-area-maps.com/img/marker-blue.png';
icon_list[2] = 'http://www.equal-area-maps.com/img/marker-gold.png';
icon_list[3] = 'http://www.equal-area-maps.com/img/marker-green.png';
var nRSSLayers=-1;

var map, scalebar;
var base_layer;
var selectControl;

// Callback for the map initialization
// Creates an OpenLayers map and hooks the basemap layer up with MapServer

function init()
{
   map = new OpenLayers.Map('map', 
                            { 
                   maxExtent: new OpenLayers.Bounds( -26000000,-26000000,
                                                      26000000, 26000000), 
                   maxResolution: 100000, units: 'meters',
                   projection: 'EPSG:3411' 
                            } );
   base_layer = new OpenLayers.Layer.MapServer( "Base Map", 
                       "http://www.yourdomain.com/mapserv.cgi?map=stere.map",
                       {layers: 'outline', minZoomLevel:2 }, 
                       {gutter: 15 }
                    );
   base_layer.buffer=1;
   map.addLayer(base_layer);

   map.setCenter( new OpenLayers.LonLat(0.0, 0.0), 1 );
   map.addControl( new OpenLayers.Control.LayerSwitcher() );

   scalebar = new OpenLayers.Control.ScaleBar( { minWidth:100, 
                                                 maxWidth:300 
                                               } );
   map.addControl(scalebar);

}


 // Call back from the GeoRSS button form
function addGeoRSS()
{
  var urlObj = OpenLayers.Util.getElement('url_georss');
  var value = urlObj.value;
  var parts = value.split("/");

  var size = new OpenLayers.Size(21,25);
  var offset = new OpenLayers.Pixel(-(size.w/2), -size.h);
  nRSSLayers = (nRSSLayers+1) % 4;
  var this_icon = new OpenLayers.Icon( icon_list[nRSSLayers], 
		                       size, offset);
  var newl = new OpenLayers.Layer.GeoRSS( parts[parts.length-1], value, 
                      { icon: this_icon, 
                        projection: new OpenLayers.Projection("EPSG:4326")
                      }  );

  map.addLayer(newl);
  urlObj.value = "";
}

	  
// Call back from the KML button
function addKML()
{
  var urlObj = OpenLayers.Util.getElement('url_kml');
  var value = urlObj.value;
  var parts = value.split("/");

  var newl = new OpenLayers.Layer.GML( parts[parts.length-1], value, 
                 {  format: OpenLayers.Format.KML,
                    formatOptions: {
                       extractStyles: true,
                       extractAttributes: true
                    },
                    projection: new OpenLayers.Projection("EPSG:4326")
                 }  );

  map.addLayer( newl );
  urlObj.value = "";		

  selectControl = new OpenLayers.Control.SelectFeature( newl, 
                   {onSelect: onFeatureSelect, onUnselect: onFeatureUnselect});
  map.addControl(selectControl);
  selectControl.activate();   
}

	  
// Popup Callbacks for the latest KML layer
function onPopupClose(evt)
{
  selectControl.unselect(selectedFeature);
}
function onFeatureSelect(feature)
{
  selectedFeature = feature;
  popup = new OpenLayers.Popup.FramedCloud("chicken", 
                     feature.geometry.getBounds().getCenterLonLat(),
                     new OpenLayers.Size(100,100),
                     "<h2>"+feature.attributes.name + "</h2>" + 
                            feature.attributes.description,
                     null, true, onPopupClose);
  feature.popup = popup;
  map.addPopup(popup);
}
function onFeatureUnselect(feature)
{
   map.removePopup(feature.popup);
   feature.popup.destroy();
   feature.popup = null;
}


</script>

</head>

<body onload="init()" >

<!-- The map object holds the OpenLayers map -->
<div id="map" class="map"></div>

<p style="clear:left" ></p>

<!-- This form can be used by the end user to add the KML or GeoRSS overlays -->

<form onsubmit="return false;">
<table>
<tr><td style="wdith:90px"><b>Add Data</b></td></tr>
<tr><td style="width:90px">GeoRSS URL: </td><td><input type="text" id="url_georss" style="width:100%" value="" />
</td>
<td style="width:100px">
  <input type="submit" onclick="addGeoRSS(); return false;" style="width:100px" value="Load GeoRSS" onsubmit="addGeoRSS(); return false;" />
</td></tr>
<tr><td style="width:90px">KML URL: </td><td><input type="text" id="url_kml" style="width:100%" value="" />
</td>
<td  style="width:100px">
 <input type="submit" onclick="addKML(); return false;" style="width:100px" value="Load KML" onsubmit="addKML(); return false;" />
</td></tr>
</table>
</form>

</body>
</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.