In [1]:
from IPython.display import YouTubeVideo

YouTubeVideo('1fzQKMp_tdE')
Out[1]:

Installation instructions

#create virtual environment virtualenv --distribute geoenv source geoenv/bin/activate

pyproj

#install pyproj (python interface to PROJ.4 library) pip install pyproj

Installing basemap

# install geos sudo apt-get install libgeos-dev cd /usr/lib sudo ln -s libgeos-3.3.3.so libgeos.so sudo ln -s libgeos-3.3.3.so libgeos.so.1 wget https://github.com/matplotlib/basemap/archive/master.zip pip install master.zip

Source: https://github.com/SciTools/cartopy/issues/46, http://blog.thefrontiergroup.com.au/2012/11/wherefore-art-thou-libgeos/

Installing GDAL

sudo apt-get install libgdal-dev

Installing GDAL 1.9.1

pip install --no-install GDAL then specify where the headers are: python setup.py build_ext --include-dirs=/usr/include/gdal/ then install it: pip install --no-download GDAL

Source: http://gis.stackexchange.com/questions/28966/python-gdal-package-missing-header-file-when-installing-via-pip

Installing Fiona

pip install Fiona

Installing Shapely

pip install Shapely

pyproj

Pyrex wrapper to provide python interfaces to PROJ.4 (http://proj.maptools.org) functions.

Performs cartographic transformations (converts from longitude, latitude to native map projection x,y coordinates and vice versa, or from one map projection coordinate system directly to another).

Source: https://pypi.python.org/pypi/pyproj/

In [3]:
from pyproj import Proj
p = Proj(init='epsg:3857')
p.srs
Out[3]:
'+units=m +init=epsg:3857 '
In [4]:
p(-97.75, 30.25)
Out[4]:
(-10881480.225042492, 3535725.659799159)
In [5]:
p(-10881480.225042492, 3535725.659799159, inverse=True)
Out[5]:
(-97.75, 30.24999999999999)

Basemap

The matplotlib basemap toolkit is a library for plotting 2D data on maps in Python. It is similar in functionality to the matlab mapping toolbox, the IDL mapping facilities, GrADS, or the Generic Mapping Tools. PyNGL and CDAT are other libraries that provide similar capabilities in Python.

Basemap does not do any plotting on it’s own, but provides the facilities to transform coordinates to one of 25 different map projections (using the PROJ.4 C library). Matplotlib is then used to plot contours, images, vectors, lines or points in the transformed coordinates. Shoreline, river and political boundary datasets (from Generic Mapping Tools) are provided, along with methods for plotting them. The GEOS library is used internally to clip the coastline and polticial boundary features to the desired map projection region.

Basemap provides facilities for reading shapefiles.

Basemap is geared toward the needs of earth scientists, particular oceanographers and meteorologists. I originally wrote Basemap to help in my research (climate and weather forecasting), since at the time CDAT was the only other tool in python for plotting data on map projections. Over the years, the capabilities of Basemap have evolved as scientists in other disciplines (such as biology, geology and geophysics) requested and contributed new features.

Source: http://matplotlib.org/basemap/users/intro.html

Exercise basemap_ex.py

In [6]:
"""
Exercise: Plotting with basemap

1. Draw a world map centered on Austin, Texas (if possible) 
   in the following projections:
    a) Mercator
    b) Robinson
    c) Orthographic
    d) Azimuthal equidistant

2. Plot the following great circle routes on a global map:
    a) Hawaii to Hong Kong
    b) Hong Kong to Moscow
    c) Moscow to Havana, Cuba
    d) Havana to Quito, Ecuador
   Coordinates of these locations are given below.  Try to choose
   projection parameters that allow you to see all of the great circles at once.
   Plot black dots at the start and end points as well.

Author: Kelsey Jordahl, Enthought
Scipy 2013 geospatial tutorial

"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# (lon, lat)
austin = (-97.75, 30.25)
hawaii = (-157.8, 21.3)
hongkong = (114.16, 22.28)
moscow = (37.62, 55.75)
havana = (-82.38, 23.13)
quito = (-78.58, -0.25)

land_color = 'lightgray'
water_color = 'lightblue'
# or choose your own colors...

1- Draw a world map centered on Austin, Texas (if possible) in the following projections: a) Mercator b) Robinson c) Orthographic d) Azimuthal equidistant

2- Plot the following great circle routes on a global map: a) Hawaii to Hong Kong b) Hong Kong to Moscow c) Moscow to Havana, Cuba d) Havana to Quito, Ecuador Coordinates of these locations are given below. Try to choose projection parameters that allow you to see all of the great circles at once. Plot black dots at the start and end points as well.

In [7]:
fig, ax = subplots(figsize=(12,12))
map = Basemap(projection='merc', llcrnrlat=-80, urcrnrlat=80,
            llcrnrlon=-180, urcrnrlon=180, resolution='l')
# draw great circle route

land_color = 'lightgray'
water_color = 'lightblue'

map.fillcontinents(color=land_color, lake_color=water_color)
map.drawgreatcircle(hawaii[0],hawaii[1],hongkong[0],hongkong[1],color='b')
map.drawgreatcircle(hongkong[0],hongkong[1],moscow[0],moscow[1],color='b')
map.drawgreatcircle(moscow[0],moscow[1],havana[0],havana[1],color='b')
map.drawgreatcircle(havana[0],havana[1],quito[0],quito[1],color='b')
map.drawcoastlines()
map.drawparallels(np.arange(-90.,120.,30.))
map.drawmeridians(np.arange(0.,420.,60.))
map.drawmapboundary(fill_color=water_color)
ax.set_title('Mercator')
map.ax = ax


x, y = map(*zip(*[hawaii, hongkong, moscow, havana, quito]))
map.plot(x, y, marker='o', markersize=6, markerfacecolor='black', linewidth=0)
Out[7]:
[<matplotlib.lines.Line2D at 0x4fc2790>]
In [8]:
fig, ax = subplots(figsize=(12,12))

map = Basemap(projection='robin', llcrnrlat=-80, urcrnrlat=80,
            llcrnrlon=-180, urcrnrlon=180, resolution='l', lon_0=austin[0])
# draw great circle route

land_color = 'lightgray'
water_color = 'lightblue'

map.fillcontinents(color=land_color, lake_color=water_color)
map.drawgreatcircle(hawaii[0],hawaii[1],hongkong[0],hongkong[1],color='b')
map.drawgreatcircle(hongkong[0],hongkong[1],moscow[0],moscow[1],color='b')
map.drawgreatcircle(moscow[0],moscow[1],havana[0],havana[1],color='b')
map.drawgreatcircle(havana[0],havana[1],quito[0],quito[1],color='b')
map.drawcoastlines()
map.drawparallels(np.arange(-90.,120.,30.))
map.drawmeridians(np.arange(0.,420.,60.))
map.drawmapboundary(fill_color=water_color)
ax.set_title('Robinson')
map.ax = ax

x, y = map(*zip(*[hawaii, hongkong, moscow, havana, quito]))
map.plot(x, y, marker='o', markersize=6, markerfacecolor='black', linewidth=0)
Out[8]:
[<matplotlib.lines.Line2D at 0x6440690>]
In [9]:
fig, ax = subplots(figsize=(12,12))

map = Basemap(projection='ortho', lon_0=austin[0], lat_0=austin[1])
# draw great circle route

land_color = 'lightgray'
water_color = 'lightblue'

map.fillcontinents(color=land_color, lake_color=water_color)
map.drawgreatcircle(hawaii[0],hawaii[1],hongkong[0],hongkong[1],color='b')
map.drawgreatcircle(hongkong[0],hongkong[1],moscow[0],moscow[1],color='b')
map.drawgreatcircle(moscow[0],moscow[1],havana[0],havana[1],color='b')
map.drawgreatcircle(havana[0],havana[1],quito[0],quito[1],color='b')
map.drawcoastlines()
map.drawparallels(np.arange(-90.,120.,30.))
map.drawmeridians(np.arange(0.,420.,60.))
map.drawmapboundary(fill_color=water_color)
ax.set_title('Orthographic')
map.ax = ax

x, y = map(*zip(*[hawaii, hongkong, moscow, havana, quito]))
map.plot(x, y, marker='o', markersize=6, markerfacecolor='black', linewidth=0)
Out[9]:
[<matplotlib.lines.Line2D at 0x63d0b10>]
In [10]:
fig, ax = subplots(figsize=(12,12))

map = Basemap(projection='aeqd', lon_0=austin[0], lat_0=austin[1])
# draw great circle route

land_color = 'lightgray'
water_color = 'lightblue'

map.fillcontinents(color=land_color, lake_color=water_color)
map.drawgreatcircle(hawaii[0],hawaii[1],hongkong[0],hongkong[1],color='b')
map.drawgreatcircle(hongkong[0],hongkong[1],moscow[0],moscow[1],color='b')
map.drawgreatcircle(moscow[0],moscow[1],havana[0],havana[1],color='b')
map.drawgreatcircle(havana[0],havana[1],quito[0],quito[1],color='b')
map.drawcoastlines()
map.drawparallels(np.arange(-90.,120.,30.))
map.drawmeridians(np.arange(0.,420.,60.))
map.drawmapboundary(fill_color=water_color)
ax.set_title('Azimuthal equidistant')
map.ax = ax


x, y = map(*zip(*[hawaii, hongkong, moscow, havana, quito]))
map.plot(x, y, marker='o', markersize=6, markerfacecolor='black', linewidth=0)
Out[10]:
[<matplotlib.lines.Line2D at 0x61f0390>]