Demo: Projection Tradeoffs and Distortion - Tissot

UW Geospatial Data Analysis
CEE467/CEWA567
David Shean

All map projections distort area, distance, and angle. The cartogrophers job is to select the best projection for the data and objectives of the map.

We’re not going to get into cartopy this week, but the following code uses cartopy to reproduce the classic Tissot indicatrix example, which helps to visualize distortion in different global projections.

For the plots below, each circle has a uniform radius of 500 km, and they are distributed uniformly across the globe every 20° of latitude and longitude. Each projection distorts these circles and their relative spacing.

import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import pyproj

#Turn off warnings from cartopy
import warnings
warnings.filterwarnings('ignore')

#Function to create cartopy plot with Tissot circles for input crs
def tissotplot(crs):
    #Create figure and axes
    fig, ax = plt.subplots(subplot_kw={'projection': crs}, figsize=(8,8))
    #Define positions of the circles, spaced every 20 degrees
    lons = range(-180, 180, 20)
    #lats = range(-90, 90, 20)
    lats = range(-80, 81, 20)
    #Draw coastlines
    ax.coastlines()
    #Add gridlines
    gl = ax.gridlines(draw_labels=True, auto_inline=True, alpha=0.5, lw=0.5, linestyle=':')
    gl.xlocator = mticker.FixedLocator(lons)
    gl.ylocator = mticker.FixedLocator(lats)
    #Add tissot circles with 500 km radius
    ax.tissot(facecolor='orange', alpha=0.4, rad_km=500, lons=lons, lats=lats)
    #Title including projection name and proj string
    crs_name = str(crs).split(' ')[0].split('.')[-1]
    ax.set_title('%s EPSG:%s\nproj string: "%s"' % (crs_name, crs.to_epsg(), crs.to_proj4()), fontsize=12)
    #ax.set_global()

#Define projecitons to plot
#Full list is here: https://scitools.org.uk/cartopy/docs/latest/crs/projections.html
crs_list = [ccrs.PlateCarree(), ccrs.LambertAzimuthalEqualArea(), ccrs.LambertCylindrical(), ccrs.Robinson(), \
            ccrs.SouthPolarStereo(), ccrs.Sinusoidal(), ccrs.Orthographic(), ccrs.Mercator(), \
            ccrs.TransverseMercator(), ccrs.TransverseMercator(central_longitude=-123), ccrs.UTM(10)]

for crs in crs_list:
    tissotplot(crs);

Note: The UTM Zone 10N is actually limited to 0° to +84° over Western North America (center portion of the above plot)

pyproj.CRS.from_epsg(32610).area_of_use.bounds

(-126.0, 0.0, -120.0, 84.0)

pyproj.CRS.from_epsg(32611).area_of_use.bounds

(-120.0, 0.0, -114.0, 84.0)

OK, those are nice, now what?

Take some time to review, read a bit about the different projections, and discuss the following with your neighbor

• Which of the above projection has the greatest area distortion? How can you tell?

  • If you skipped the link earlier, maybe go back and read about what the Tissot circles mean :)

• Which has greatest azimuthal distortion?

• If you were crossing the Pacific Ocean in the 1800s, what projection might you want for your paper maps needed for navigation?

This might be a useful resource: https://www.esri.com/arcgis-blog/products/product/mapping/tissots-indicatrix-helps-illustrate-map-projection-distortion/

Other useful resources

• Nice interactive visualization of corresponding areas on sphere and projected map: https://mathigon.org/course/circles/spheres-cones-cylinders#sphere-maps

• Interactive visualization and quantiative information about distortion for many projections: https://bl.ocks.org/syntagmatic/ba569633d51ebec6ec6e

• Map projection lava lamp: https://www.jasondavies.com/maps/transition/

• Also, https://xkcd.com/977/

We will regroup and discuss interpretations as a group