Equal-Area Map Projections with Basemap and matplotlib/pylab

Di, 06/25/2013 - 17:23 — Draketo

PDF (read as slides)

Org (reproduce)

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

1. Problem
- 1.1. lat/lon pixels misrepresent areas
2. Map-Notes
- 2.1. Map Projections
3. Plotted
4. Other maps
- 4.1. Other Equal Area map types
5. Conclusing
- 5.1. Maps I plan to use
6. Thank you!
- 6.1. Thank you for listening!
7. Appendix: Supporting functions

1 Problem

1.1 lat/lon pixels misrepresent areas

Projected B_columns
1. Simple Flat BMCOL B_block
2. Globe BMCOL
Note B_ignoreheading

Sibiria 13,1 · 10⁶ km² vs. china 9.7 · 10⁶ 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

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 lat_ts=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$")

3.2 Hammer

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

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$")

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

Maps B_columns
1. Hobo-Dyer BMCOL B_block
  
  To show regional fluxes and longitudinally constrained regions: Easy to spot on rectangular grid.
2. Hammer B_block BMCOL
  
  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
  
  For mainly latitudinally constrained regions: Straight latitudinal lines and high latitudinal resolution near the poles.

6 Thank you!

6.1 Thank you for listening!

Questions? B_block BMCOL

7 Appendix: Supporting functions

7.1 Basemap Imports

Basic imports for plotting with Basemap.

# 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

Basemap: http://matplotlib.org/basemap/
Pylab: http://www.scipy.org/PyLab
Numpy: http://www.numpy.org/
NetCDF: http://code.google.com/p/netcdf4-python/

7.2 Draw a map

Commands to draw a map. If no emission data is available, use the default fill color for continents. Requires previous imports and a basemap setup as done in the plot snippets for the individual maps.

<<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

Commands to add visual orientation to a map.

# 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

Add Tissots Indicatrix to visualize distortion.

# 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

Plot Carbontracker emissions on a basemap (m).

# 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

Add a nice colorbar to the plot.

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

Anhang	Größe
flatpolarquartic.png	127.57 KB
hobo-dyer.png	139.26 KB
hammer.png	138.15 KB
sibiria-china-flat.png	1.14 MB
sibiria-china-globe.png	1.14 MB
equal-area-map-projections.pdf	3.06 MB
equal-area-map-projections.org	10.19 KB

Druckversion
Login to post comments

Use Node:

⚙ Babcom is trying to load the comments ⚙

This textbox will disappear when the comments have been loaded.

If the box below shows an error-page, you need to install Freenet with the Sone-Plugin or set the node-path to your freenet node and click the Reload Comments button (or return).

If you see something like Invalid key: java.net.MalformedURLException: There is no @ in that URI! (Sone/search.html), you need to setup Sone and the Web of Trust

If you had Javascript enabled, you would see comments for this page instead of the Sone page of the sites author.

Note: To make a comment which isn’t a reply visible to others here, include a link to this site somewhere in the text of your comment. It will then show up here. To ensure that I get notified of your comment, also include my Sone-ID.

Link to this site and my Sone ID: sone://6~ZDYdvAgMoUfG6M5Kwi7SQqyS-gTcyFeaNN1Pf3FvY

This spam-resistant comment-field is made with babcom.