Dec 11 2005
Tissot Indicatrix - Examining the distortion of map projections
The Tissot Indicatrix is a valuable tool for showing the distortions caused by map projections. It is essentially a series of imaginary polygons that represent perfect circles of equal area on a 3D globe. When projected onto a 2D map, their shape, size and/or angles will be distorted accordingly allowing you to quickly assess the projection’s accuracy for a given part of the globe.
I’ve seen great Tissot diagrams in text books but I wanted to create the indicatrix as a polygon dataset so that I could project and overlay it with other data in a GIS. To do this I wrote a python script using the OGR libraries, which I will revist in a minute. But first the visually interesting part:
Here is a world countries shapefile overlaid with the Tissot circles in geographic (unprojected lat-long) coordinates:
Next I reprojected the datasets to the Mercator projection using ogr2ogr:
ogr2ogr -t_srs “+proj=merc” countries_merc.shp countries_simpl.shp countries_simpl
ogr2ogr -t_srs “+proj=merc” tissot_merc.shp tissot.shp tissot
Note that the angles are perfectly preserved (the trademark feature of the Mercator projection) but the size is badly distorted.
Now lets try Lambert Azimuthal Equal Area (in this case the US National Atlas standard projection - EPSG code 2163).
ogr2ogr -t_srs “epsg:2163″ countries_lambert.shp countries_simpl.shp countries_simpl
ogr2ogr -t_srs “epsg:2163″ tissot_lambert.shp tissot.shp tissot
This is a great projection for preserving area but get outside the center and shapes become badly distorted:
The best way to experiment with this is to bring the tissot.shp file into ArcMap (or another program that supports on-the-fly projection) and play with it in real time. The distortions of every projection just leap off the screen…
OK, now for the geeky part. Here’s the python/OGR script used to create the tissot shapefile. The basic process is to lay out a grid of points across the globe in latlong, loop through the points and reproject each one to an orthographic projection centered directly on the point, buffer it, then reproject to latlong. The end result is a latlong shapefile representing circles of equal area on a globe.
#!/usr/bin/env python
# Tissot Circles
# Represent perfect circles of equal area on a globe
# but will appear distorted in ANY 2d projection.
# Used to show the size, shape and directional distortion
# by Matthew T. Perry
# 12/10/2005
import ogr
import os
import osr
output = 'tissot.shp'
debug = False
# Create the Shapefile
driver = ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists(output):
driver.DeleteDataSource(output)
ds = driver.CreateDataSource(output)
layer = ds.CreateLayer(output, geom_type=ogr.wkbPolygon)
# Set up spatial reference systems
latlong = osr.SpatialReference()
ortho = osr.SpatialReference()
latlong.ImportFromProj4('+proj=latlong')
# For each grid point, reproject to ortho centered on itself,
# buffer by 640,000 meters, reproject back to latlong,
# and output the latlong ellipse to shapefile
for x in range(-165,180,30):
for y in range (-60,90,30):
f= ogr.Feature(feature_def=layer.GetLayerDefn())
wkt = 'POINT(%f %f)' % (x, y)
p = ogr.CreateGeometryFromWkt(wkt)
p.AssignSpatialReference(latlong)
proj = '+proj=ortho +lon_0=%f +lat_0=%f' % (x,y)
ortho.ImportFromProj4(proj)
p.TransformTo(ortho)
b = p.Buffer(640000)
b.AssignSpatialReference(ortho)
b.TransformTo(latlong)
f.SetGeometryDirectly(b)
layer.CreateFeature(f)
f.Destroy()
ds.Destroy()
This is a neat way of immediately understanding what a Tissot Indicatrix is actually is used for!
I am not familiar with Python but don’t you have to do layer.destroy() and output.close() or something like that?
I’m fairly sure that the ds.Destroy() at the end takes care of closing everything out.
Hey Matt,
Super impressive script! What a clever idea. I am trying to reproduce this on a student version of ArcGIS 9. I think I am missing a library (ImportError: No module named ogr) I guess this is not available on the student version? I will try at work and see if it works.
you would need GDAL/OGR for the library to be present.
awsome program: I was able to modify it for a quick visual analysis of the distortion in one of our online projects:
http://169.237.35.250/~dylan/temp/aea-tissot_1×1_deg.png
thanks.
Fantastic! I’ve got to do a lecture on projections in a very basic maps class - will load and run your script into ArcGIS later but don’t have time right now - in the meantime, may I use your graphics? Properly credited, of course?
That is a very useful shapefile (and script). I can see using this to show a client or students why we use different projections.
Many thanks for your very interesting blog. The particular post helped me a lot, because I needed to create a “mapplet” for my Mapserver application to overlay a distortion estimation. I would like to post my method because it is possible that someone will need an integration with mapserver.
def tissot_layer(path, outproj, map):
”’ This function creates a tissot Indicatrix layer and add this
layer to a mapscript layer (layerObj)
@param path: the path for the Tissot shapefile to be created
@param outproj: The proj4 definition of the Mapserver map level
projection. Accessible through mapscript’s mapObj.getProjection()
@param map: The mapscript mapObj to hold the new layer
@rtype: String
”’
output = path
debug = False
# Create the Shapefile
driver = ogr.GetDriverByName(’ESRI Shapefile’)
if os.path.exists(output):
driver.DeleteDataSource(output)
ds = driver.CreateDataSource(output)
layer = ds.CreateLayer(output, geom_type=ogr.wkbPolygon)
# Set up spatial reference systems
latlong = osr.SpatialReference()
ortho = osr.SpatialReference()
map_projection = osr.SpatialReference()
latlong.ImportFromProj4(’+proj=latlong’)
# The output projection we need
map_projection.ImportFromProj4(outproj)
# For each grid point, reproject to ortho centered on itself,
# buffer by 640,000 meters, reproject back to latlong,
# and output the latlong ellipse to shapefile
for x in range(-180,180,30):
for y in range (-90,90,30):
f= ogr.Feature(feature_def=layer.GetLayerDefn())
wkt = ‘POINT(%f %f)’ % (x, y)
p = ogr.CreateGeometryFromWkt(wkt)
p.AssignSpatialReference(latlong)
proj = ‘+proj=ortho +lon_0=%f +lat_0=%f’ % (x,y)
ortho.ImportFromProj4(proj)
p.TransformTo(ortho)
b = p.Buffer(640000)
b.AssignSpatialReference(ortho)
# Here the original script reprojects the feature back to
# latlong. But we must reproject it to the map’s output
# projection using its Proj4 definition.
# TESTED WITH QGIS: WORKS (EXAMPLE LAEA PROJECTION 2163)
b.TransformTo(map_projection)
f.SetGeometryDirectly(b)
layer.CreateFeature(f)
f.Destroy()
ds.Destroy()
# By now we have a “polygon” layer consisting of features in
# the map’s proj4 definition. Now we must open the dataset
# and populate the layer with the Tissot features
_layer = mapscript.layerObj(map)
_layer.name = “tissot”
_layer.type = mapscript.MS_LAYER_POLYGON
_class = mapscript.classObj(_layer)
_style = mapscript.styleObj(_class)
_style.color.setRGB(255,0,0)
_style.outlinecolor.setRGB(0,0,0)
_layer.status = mapscript.MS_ON
_layer.opacity = 50
_layer.setProjection(outproj)
_layer.data = path
From now on the map can be drawn and include the layer
ooops as i see tabs don’t work. sorry!
Just curious as to why you call
CreateLayer(output, geom_type = ogr.wkbPolygon)
instead of
CreateLayer(output, ogr.wkbPolygon)
which seems to cause an assertion
Don, I believe the CreateLayer method has a number of other optional arguments so to avoid having to specify the details for all of them in order, I just set geom_type explicitly.
Sorry, but this is not how you generate the Tissot indicatrix. The parameters of the ellipse are defined by the properties of the projection AT a point. The crucial properties are the maximum and minimum POINT scales and the direction of the direction of maximum scale. This is enough to construct an indicatrix which is always an exact ellipse. For convenience this exact ellipse is magnified and overlain on the projection map. The distorted shapes shown in the outer regions of the azimuthal example are quite impossible. Projecting finite circles (which are not described by point parameters) is simply not correct. Of course your shapes are reasonably ok until you get to large distortions.
In most current 2D GIS data models (ie based on the Simple Features spec), a polygon area must be represented as a closed series of vertices. In order to allow for magnifying and overlaying the ellipse in real coordinate space, you must approximate the tissot ellipses using polygons. So while it may not be a mathematically precise tissot indicatrix, it allows GIS users to simulate it which was the entire point of this exercise.
Hi
I’m a newbie so forgive the obvious question :
import ogr
import os
import osr
How do I get this to work, from what I can gather I need to down load this (Which I’ve done from http://fwtools.maptools.org/ )
I’m using Arcview (info licence)
Windows7 (64bit version)
but I’m getting the error
Running script Script…
: No module named ogr
Failed to execute (Script).
End Time: Mon May 03 10:15:09 2010
Justin,
First off, ArcView has nothing to do with it .. this is a completely ESRI-free solution built on open source software.
Secondly, the script was written 5 years ago so the libraries may have changed slightly. The one instance I am sure of is that ogr and osr are now part of the osgeo library so they are imported like so:
from osgeo import ogr
from osgeo import osr
If those dont work, you need to confirm that you have installed the GDAL python bindings for the version of python that you are running. If using FWTools, you must use the version of python that comes with FWTools.