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

GDAL

GDAL (Geospatial Data Abstraction Library) is a library for reading and writing raster geospatial data formats, and is released under the permissive X/MIT style free software license by the Open Source Geospatial Foundation. As a library, it presents a single abstract data model to the calling application for all supported formats. It may also be built with a variety of useful command-line utilities for data translation and processing. The related OGR library (OGR Simple Features Library[2]), which is part of the GDAL source tree, provides a similar capability for simple features vector data.

Source: https://en.wikipedia.org/wiki/GDAL

GDAL 1.9.1 (Python)

This Python package and extensions are a number of tools for programming and manipulating the GDAL Geospatial Data Abstraction Library. Actually, it is two libraries -- GDAL for manipulating geospatial raster data and OGR for manipulating geospatial vector data -- but we'll refer to the entire package as the GDAL library for the purposes of this document.

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

Exercise geotiff.py

In [11]:
"""
Exercise: Read a geotiff file as a numpy array and display with matplotlib.

1. Download the data file from http://public.enthought.com/~kjordahl/scipy/manhattan.tif

2. Open the file with GDAL.  What projection is it in?

3. Read the image data into a numpy array and plot as a 3-color image
   with matplotlib.

4. Set the plot axis limits to the proper map coordinates.

BONUS

5. plot the locations of the Citibike stations in the file citibike.json
   (or real-time from http://citibikenyc.com/stations/json)
   
Author: Kelsey Jordahl, Enthought
Scipy 2013 geospatial tutorial

"""

import os
from osgeo import gdal
import matplotlib.pyplot as plt
# GDAL does not use python exceptions by default
gdal.UseExceptions()
In [13]:
! wget http://public.enthought.com/~kjordahl/scipy/manhattan.tif
--2013-07-13 00:59:44--  http://public.enthought.com/~kjordahl/scipy/manhattan.tif
Resolving public.enthought.com (public.enthought.com)... 50.17.225.20
Connecting to public.enthought.com (public.enthought.com)|50.17.225.20|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 119987605 (114M) [image/tiff]
Saving to: ‘manhattan.tif’

100%[======================================>] 119.987.605 2,03MB/s   in 51s    

2013-07-13 01:00:35 (2,26 MB/s) - ‘manhattan.tif’ saved [119987605/119987605]

2- Open the file with GDAL. What projection is it in?

In [14]:
gtif = gdal.Open('manhattan.tif')
gtif.GetProjectionRef()
Out[14]:
'PROJCS["UTM",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.2572221010002,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["meters",1],AUTHORITY["EPSG","26918"]]'

3- Read the image data into a numpy array and plot as a 3-color image with matplotlib.

4- Set the plot axis limits to the proper map coordinates.

In [15]:
arr = gtif.ReadAsArray()
trans = gtif.GetGeoTransform()
extent = (trans[0], trans[0] + gtif.RasterXSize*trans[1],
          trans[3] + gtif.RasterYSize*trans[5], trans[3])

plt.imshow(arr[:3,:,:].transpose((1, 2, 0)), extent=extent)
plt.show()