Graticule Hacking with OGR
MapServer has a built-in graticule object that can be configured in the same way as its scalebar and reference map. PCL's MapServer-based engine, however, is not going to use the built-in graticule. Users will instead be required to render a graticule layer from a set of parallel and meridian lines that I do not just happen to have lying around. In the past I have scripted such processes using the shapelib module, but for some practice with OGR I decided to write a simple script using ogr.py to generate my annotated graticule lines. OGR's pymod/samples/tigerpoly.py was my reference for creating and writing a shapefile of features. Right after I finished, it was pointed out to me that I'd overlooked the pymod/mkgraticule.py script recently committed to CVS. Still, mine has attributes for labeling. After a few more tweaks, I'll use the following script to generate parallel and meridian data to be distributed with PCL.
The graticule must consist of individual short line segments for use with MapServer, which throws out features that cannot be completely projected. This occurs in the case of Orthographic, and possibly other hemispheric type projections. The multitude of lines complicates labeling a bit, but that may be something I can overcome as the script is refined.
We begin by creating shapefile to contain our graticule line segments:
import math import ogr # Create an ESRI shapefile of parallels and meridians for a MapServer # world map. filename = 'graticule.shp' driver = ogr.GetDriverByName('ESRI Shapefile') driver.DeleteDataSource(filename) ds = driver.CreateDataSource(filename) layer = ds.CreateLayer('graticule', geom_type=ogr.wkbLineString)
Next, define several attribute fields. A 'TYPE' field will let me classify and render differently the parallels and meridians if I so choose. 'VALUE' and 'ABS_VALUE' fields will be useful for labeling the parallels and meridians:
# Create an integer field for classification # 0: parallel # 1: meridian fd = ogr.FieldDefn('TYPE', ogr.OFTInteger) fd.SetWidth(1) fd.SetPrecision(0) layer.CreateField(fd) # Create two fields for labeling fd = ogr.FieldDefn('VALUE', ogr.OFTInteger) fd.SetWidth(4) fd.SetPrecision(0) layer.CreateField(fd) fd = ogr.FieldDefn('ABS_VALUE', ogr.OFTInteger) fd.SetWidth(4) fd.SetPrecision(0) layer.CreateField(fd)
Next, create parallels and insert them into the layer:
# First, the parallels at 10 degree intervals, with one degree resolution for j in range(1,18): y = 10*float(j-9) for i in range(0, 360): line = ogr.Geometry(type=ogr.wkbLineString) # hack: MapServer has trouble within .1 decimal degrees of the # dateline if i == 0: x1 = -179.9 x2 = -179.0 elif i == 359: x1 = 179.0 x2 = 179.9 else: x1 = float(i-180) x2 = x1 + 1.0 line.AddPoint(x1, y) line.AddPoint(x2, y) f = ogr.Feature(feature_def=layer.GetLayerDefn()) f.SetField(0, 0) f.SetField(1, y) f.SetField(2, math.fabs(y)) f.SetGeometryDirectly(line) layer.CreateFeature(f) f.Destroy()
And now the meridians:
# Next, the meridians at 10 degree intervals and one degree resolution for i in range(0, 37): x = 10*float(i-18) # hack: MapServer has trouble within .1 decimal degrees of the # dateline if i == 0: x = -179.9 if i == 36: x = 179.9 for j in range(10, 170): line = ogr.Geometry(type=ogr.wkbLineString) y1 = float(j - 90) y2 = y1 + 1.0 line.AddPoint(x, y1) line.AddPoint(x, y2) f = ogr.Feature(feature_def=layer.GetLayerDefn()) f.SetField(0, 1) f.SetField(1, x) f.SetField(2, math.fabs(x)) f.SetGeometryDirectly(line) layer.CreateFeature(f) f.Destroy() # destroying data source closes the output file ds.Destroy()
Because the Python Cartographic Library doesn't yet support non-EPSG coordinate systems (coming soon), I'm going to use MapServer's scripting module to draw the graticules:
import mapscript mo = mapscript.mapObj() mo.setFontSet('/home/sean/projects/ms_44/mapserver/tests/fonts.txt') mo.selectOutputFormat('image/png') layer = mapscript.layerObj() layer.status = mapscript.MS_DEFAULT layer.data = 'graticule.shp' layer.setProjection('+init=epsg:4326') layer.type = mapscript.MS_LAYER_LINE layer.labelitem = 'VALUE' layer.classitem = 'TYPE' li = mo.insertLayer(layer) layer = mo.getLayer(li) # render parallels and meridians in same style style = mapscript.styleObj() style.color.setHex('#808080') co = mapscript.classObj() co.label.type = mapscript.MS_TRUETYPE co.label.font = 'Vera' co.label.size = 6 co.label.position = mapscript.MS_AUTO co.label.mindistance = 200 co.insertStyle(style) co.setExpression('([TYPE] EQ 0)') co.label.color.setRGB(0, 255, 0) layer.insertClass(co) style.color.setHex('#0000FF') co.insertStyle(style) co.setExpression('([TYPE] EQ 1)') co.label.color.setRGB(255, 0, 0) layer.insertClass(co)
Next, draw the graticule using a Mollweide projection:
# Mollweide Projection mo.setSize(600, 300) mo.setProjection('+proj=moll +lon_0=0.0') mo.setExtent(-18100000, -9050000, 18100000, 9050000) im = mo.draw() im.save('graticule-moll.png')
Looking good, graticule-moll.png:

And finally, an Orthographic projection:
mo.setSize(400, 400) mo.setProjection('+proj=ortho +lon_0=0.0 +lat_0=40.0') mo.setExtent(-6400000,-6400000,6400000,6400000) im = mo.draw() im.save('graticule-ortho.png')
Looks like we have an ortho-related bug in MapServer or PROJ.4, graticule-ortho.png:
