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

Fiona

Fiona provides a minimal, uncomplicated Python interface to the open source GIS community's most trusted geodata access library and integrates readily with other Python GIS packages such as pyproj, Rtree, and Shapely.

How minimal? Fiona can read feature records as mappings from shapefiles or other GIS vector formats and write mappings as records to files using the same formats. That's all. There aren't any feature or geometry classes. Records and their geometries are just data.

Source: https://github.com/sgillies/Fiona

Exercise geology.py

In [16]:
"""
Exercise: Reading a shapefile with Fiona

1. Create a list of all the unique rock types in the data
   (in properties ROCKTYPE1 and ROCKTYPE2).

2. Calculate the total area of each primary rocktype (ROCKTYPE1) by summing
   the property AREA.

3. Plot the polygons in the data with basemap, coloring by primary rock type.

BONUS:

4. Calculate the total area again, this time by using Shapely to calculate the
   area of each polygon.  HINT: You may want to use New Jersey State Plane
   coordinates, EPSG:32111


Author: Kelsey Jordahl, Enthought
Scipy 2013 geospatial tutorial

"""

import os
import zipfile

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from collections import defaultdict
import fiona

1- Create a list of all the unique rock types in the data (in properties ROCKTYPE1 and ROCKTYPE2).

In [19]:
rocks = []
with fiona.open('njgeol_poly_dd.shp') as f:
    for rec in f:
        rocks.append( rec['properties']['ROCKTYPE1'] )
        rocks.append( rec['properties']['ROCKTYPE2'] )
len(set(rocks))
Out[19]:
48

2- Calculate the total area of each primary rocktype (ROCKTYPE1) by summing the property AREA.

In [20]:
from collections import defaultdict
d = defaultdict(float)
with fiona.open('njgeol_poly_dd.shp') as f:
    for rec in f:
        d[rec['properties']['ROCKTYPE1']] += rec['properties']['AREA']
In [21]:
fig, ax = subplots(figsize=(12,12))
ax.bar(arange(len(d.values())), d.values())
_ = ax.set_xticks(arange(len(d.keys())))
_ = ax.set_xticklabels( d.keys(), rotation='vertical' )

3- Plot the polygons in the data with basemap, coloring by primary rock type.

In [22]:
fig, ax = subplots(figsize=(12,12))
west, east, south, north = -75.6, -73.5, 38.5, 41.5
m = Basemap(projection='merc', llcrnrlat=south, urcrnrlat=north,
           llcrnrlon=west, urcrnrlon=east, lat_ts=0, resolution='l')

colormap = defaultdict(lambda: np.random.random(3))
with fiona.open('njgeol_poly_dd.shp') as f:
    for idx,rec in enumerate(f):
        coords = rec['geometry']['coordinates'][0]
        rocktype = rec['properties']['ROCKTYPE1']
        x, y = m(*zip(*coords))
        poly = Polygon(zip(x,y), facecolor=colormap[rocktype])
        ax.add_patch(poly)
        
m.drawmapboundary()
m.drawcoastlines()
Out[22]:
<matplotlib.collections.LineCollection at 0xa2dd710>

Shapely

Shapely is a BSD-licensed Python package for manipulation and analysis of planar geometric objects. It is based on the widely deployed GEOS (the engine of PostGIS) and JTS (from which GEOS is ported) libraries. This C dependency is traded for the ability to execute with blazing speed. Shapely is not concerned with data formats or coordinate systems, but can be readily integrated with packages that are.

Source: https://github.com/sgillies/shapely

Exercise nyc.py

In [24]:
"""

Exercise: Boroughs of New York City

1. Read the borough boundaries shapefile from data/nybb_13a/nybb.shp
   into a dictionary of Shapely geometries keyed by the property BoroName.
   
2. Calculate the area of each borough.  What are the units?

3. Calculate the fraction of the area of each borough that lies more than 1 km
   (3281 feet) from its boundary

BONUS

4. Extract the simple polygon representing the island (not the borough) of
   Manhattan.  HINT: It will be the largest polygon in the borough.

Author: Kelsey Jordahl, Enthought
Scipy 2013 geospatial tutorial

"""
import os
import numpy as np
import fiona
from shapely.geometry import shape

def plot_polygon(ax, poly, color='red'):
    a = np.asarray(poly.exterior)
    ax.add_patch(Polygon(a, facecolor=color, alpha=0.3))
    ax.plot(a[:, 0], a[:, 1], color='black')

def plot_multipolygon(ax, geom, color='red'):
    """ Can safely call with either Polygon or Multipolygon geometry
    """
    if geom.type == 'Polygon':
        plot_polygon(ax, geom, color)
    elif geom.type == 'MultiPolygon':
        for poly in geom.geoms:
            plot_polygon(ax, poly, color)

1- Read the borough boundaries shapefile from data/nybb_13a/nybb.shp into a dictionary of Shapely geometries keyed by the property BoroName.

In [27]:
nyc_geom = defaultdict()
colors = ['red', 'green', 'orange', 'brown', 'purple']

fig, ax = subplots(figsize=(12,12))

with fiona.open('nybb.shp') as f:
    crs = f.crs
    units = crs['units']
    print 'units', units
    for rec in f:
        color = colors.pop()
        print rec['geometry']['type']
        boro = rec['properties']['BoroName']
        nyc_geom[boro] = shape(rec['geometry'])
        plot_multipolygon(ax, nyc_geom[boro], color=color)

labels = ax.get_xticklabels() 
for label in labels: 
    label.set_rotation(90) 
units us-ft
MultiPolygon
MultiPolygon
MultiPolygon
MultiPolygon
MultiPolygon

2- Calculate the area of each borough. What are the units?

In [28]:
for boro, geom in nyc_geom.iteritems():
    print '%s area %10.0f square feet (%5.2f sq miles)' % (boro,
                                                           geom.area,
                                                           (geom.area / 27878400))
Bronx area 1186805997 square feet (42.57 sq miles)
Brooklyn area 1959433451 square feet (70.29 sq miles)
Staten Island area 1623855480 square feet (58.25 sq miles)
Manhattan area  636441883 square feet (22.83 sq miles)
Queens area 3049948268 square feet (109.40 sq miles)

3 - Calculate the fraction of the area of each borough that lies more than 1 km (3281 feet) from its boundary

In [29]:
figure, ax = subplots(figsize=(12,12))
for boro, geom in nyc_geom.iteritems():
    interior = geom.buffer(-3281)
    plot_multipolygon(ax, interior, color='gray')
    print '%4.1f%% of %s more than 1 km from boundary' % (100 * interior.area / geom.area,
                                                          boro)
46.0% of Bronx more than 1 km from boundary
50.2% of Brooklyn more than 1 km from boundary
62.7% of Staten Island more than 1 km from boundary
25.0% of Manhattan more than 1 km from boundary
56.7% of Queens more than 1 km from boundary

4- Extract the simple polygon representing the island (not the borough) of Manhattan. HINT: It will be the largest polygon in the borough.

In [30]:
fig, ax = subplots(figsize=(12,12))
idx = argmax([s.area for s in nyc_geom['Manhattan'].geoms])
manhattan_i = nyc_geom['Manhattan'].geoms[idx]
plot_polygon(ax, manhattan_i)