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
Inhalt abgleichen
Willkommen im Weltenwald!
((λ()'Dr.ArneBab))



Beliebte Inhalte

sn.1w6.org news

Kommunikation im Internet ist in Gefahr!

Liebe Besucher,
leider müssen wir euch kurz unterbrechen:

Das EU-Parlament stimmt nächste Woche über eine Urheberrechtsreform ab,
die alle kleinen Kommunikationsplattformen in ihrer Existenz bedroht.

Offener Brief der Forenbetreiber

Wir bitten euch daher darum, dass ihr

Vielen Dank für eure Unterstützung!

Euer Draketo

Wir sind Teil des #321EUOfflineDay