Equal-Area Map Projections with Basemap and matplotlib/pylab

PDF (read as slides)

Org (reproduce)

Plotting global equal area maps with python, matplotlib/pylab and Basemap.

Table of Contents

1 Problem

1.1 lat/lon pixels misrepresent areas

  1. Projected   B_columns
    1. Simple Flat   BMCOL B_block

      sibiria-china-flat.png


    2. Globe   BMCOL

      sibiria-china-globe.png



  2. Note   B_ignoreheading

    Sibiria 13,1 · 106 km² vs. china 9.7 · 106 km²

    Maps thanks to Marble Desktop Globe and Open Street Map, available under CC by-sa and Open Data Commons Open Database License (ODbL).


2 Map-Notes

2.1 Map Projections

  1. Maps   B_columns
    1. Hobo-Dyer   BMCOL B_block
      • Rectangle
      • equidistant longitude
      • longitude/latitude over the Mediterranean Sea (more exactly: \(37.5^\circ\))
      • Similar Maps: Gall-Peters (thinner), Lambert (wider)
      • Basemap: equal area cylindrical (cea) with latts=37.5

    2. Hammer   BMCOL B_block
      • Elliptic
      • Low distortion at the poles
      • 2:1 → 2 per page
      • the earth appears round without making it hard to recognize patterns
      • Similar Maps: Mollweide (more distorted at the poles, parallel latitudes)
      • Basemap: hammer

    3. Flat Polar Quartic   B_block BMCOL
      • Elliptic with polar cuts
      • parallel latitudes
      • Standard parallels at \(33^\circ 45' N/S\)
      • poles are \(\frac{1}{3}\) the equator
      • Similar: Eckart IV (poles are half the equator)
      • Basemap: mbtfpq


3 Plotted

3.1 Hobo-Dyer

m = map.Basemap(projection='cea', lat_ts=37.5)
outfile = "hobo-dyer.png"
pl.title("Hobo Dyer: Cylindric Equal Area at $37.5^\\circ N$")

hobo-dyer.png

3.2 Hammer

m = map.Basemap(projection='hammer', lon_0=0)
outfile = "hammer.png"
pl.title("Hammer") # latex-test: $\frac{1}{2}$

hammer.png

3.3 Flat Polar Quartic

m = map.Basemap(projection='mbtfpq', lon_0=0)
outfile = "flatpolarquartic.png"
pl.title("Flat Polar Quartic: parallels at $33^\\circ 45' N/S$")

flatpolarquartic.png

4 Other maps

4.1 Other Equal Area map types

  • Goode homolosine: Split, focus on land or ocean, straight latitude parallels, approximately preserve most shapes. Not available in matplotlib.
  • Eckert IV: Like Flat Polar Quartic, parallels at 40° 30' N/S, poles are half the equator.
  • Lambert cylindrical equal area: Like Hobo Dyer, very wide, shapes at the equator are correct.
  • Gall Peters: Like Hobo Dyer, appears more distorted than Hobo-Dyer, shapes over europe correct (45°).
  • Mollweide: Like Hammer with straight latitude parallels.
  • Werner: It’s a heart :) - focus on a hemisphere without ignoring the rest. General case: Bonne.
  • Tobler: General case leading to Lambert, Mollweide, Mollington and a few more — also see Tobler1973 after you manage to gnawl through the paywall…
  • Collignon: Triangle, for cosmic microwave background.

5 Conclusing

5.1 Maps I plan to use

  1. Maps   B_columns
    1. Hobo-Dyer   BMCOL B_block

      hobo-dyer.png

      To show regional fluxes and longitudinally constrained regions: Easy to spot on rectangular grid.


    2. Hammer   B_block BMCOL

      hammer.png

      To show a global overview: Helps the understanding of global data because it appears most similar to a real earth while including the whole earth surface.


    3. Flat Polar Quartic   B_block BMCOL

      flatpolarquartic.png

      For mainly latitudinally constrained regions: Straight latitudinal lines and high latitudinal resolution near the poles.



6 Thank you!

6.1 Thank you for listening!

  1. Questions?   B_block BMCOL

7 Appendix: Supporting functions

7.1 Basemap Imports

# basemap, pylab and numpy for plotting
import mpl_toolkits.basemap as map
import pylab as pl
import numpy as np
# netcdf for reading the emission files
import netCDF4 as nc

7.2 Draw a map

<<addmapfeatures>>
<<addindicatrix>>
try:
  <<addemissions>>
  <<addcolorbar>>
except RuntimeError: # recover from missing fluxfile
  m.fillcontinents(color='coral',lake_color='aqua')
pl.savefig(outfile)
return "./" + outfile + ""

7.3 Map features

# add map lines
m.drawcoastlines()
# only fill continents if we do not plot emissions
# m.fillcontinents(color='coral',lake_color='aqua')
m.drawparallels(np.arange(-90.,120.,30.), 
                labels=[False,True,True,False])
m.drawmeridians(np.arange(0.,420.,60.), 
                labels=[True,False,False,True])
m.drawmapboundary(fill_color='aqua')

7.4 Tissots Indicatrix

# draw tissot's indicatrix to show distortion.
for y in np.linspace(m.ymax/20,19*m.ymax/20,9):
    for x in np.linspace(m.xmax/20,19*m.xmax/20,12):
        lon, lat = m(x,y,inverse=True)
        poly = m.tissot(lon,lat,4.0,100,
                        facecolor='green',
                        zorder=10,alpha=0.5)

7.5 Plot emissions

# d = nc.Dataset("/run/media/arne/3TiB/CTDAS-2013-03-07-2years-base-data/"
#                "analysis/data_flux1x1_weekly/flux_1x1.nc")
d = nc.Dataset("UNPUBLISHED")
biocovmean = np.mean(
    d.variables["bio_flux_prior_cov"][:,:,:], axis=0)
# projection: matplotlib.org/basemap/users/examples.html
lons, lats = pl.meshgrid(range(-180, 180), 
                         range(-90, 90))
x, y = m(lons, lats)
# choose my standard color range: vmin = -0.5*vmax
vmax = max(abs(np.max(biocovmean)), 
           2 * abs(np.min(biocovmean)))
vmin = -0.5*vmax
m.pcolor(x, y, biocovmean, shading='flat', 
         vmin=vmin, vmax=vmax) # pcolormesh is faster

7.6 Nice colorbar

pl.rcParams.update({"text.usetex": True, 
                    "text.latex.unicode": True})
colorbar = pl.colorbar(orientation="horizontal", 
                       format="%.2g") # scientific
colorbar.set_label("$CO_{2}$ fluxes [$\\frac{mol}{m^2 s}$]")

Author: Arne Babenhauserheide

Emacs 24.3.1 (Org mode 8.0.2)

Validate XHTML 1.0

AnhangGröße
flatpolarquartic.png127.57 KB
hobo-dyer.png139.26 KB
hammer.png138.15 KB
sibiria-china-flat.png1.14 MB
sibiria-china-globe.png1.14 MB
equal-area-map-projections.pdf3.06 MB
equal-area-map-projections.org10.19 KB