from IPython.display import YouTubeVideo
YouTubeVideo('1fzQKMp_tdE')
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
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).
from pyproj import Proj
p = Proj(init='epsg:3857')
p.srs
p(-97.75, 30.25)
p(-10881480.225042492, 3535725.659799159, inverse=True)
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.
Exercise basemap_ex.py¶
"""
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.
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)
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)
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)
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)
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.
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.
Exercise geotiff.py¶
"""
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()
1- Download the data file from http://public.enthought.com/~kjordahl/scipy/manhattan.tif
! wget http://public.enthought.com/~kjordahl/scipy/manhattan.tif
2- Open the file with GDAL. What projection is it in?
gtif = gdal.Open('manhattan.tif')
gtif.GetProjectionRef()
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.
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.
Exercise geology.py¶
"""
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).
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))
2- Calculate the total area of each primary rocktype (ROCKTYPE1) by summing the property AREA.
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']
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.
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()
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.
Exercise nyc.py¶
"""
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.
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)
2- Calculate the area of each borough. What are the units?
for boro, geom in nyc_geom.iteritems():
print '%s area %10.0f square feet (%5.2f sq miles)' % (boro,
geom.area,
(geom.area / 27878400))
3 - Calculate the fraction of the area of each borough that lies more than 1 km (3281 feet) from its boundary
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)
4- Extract the simple polygon representing the island (not the borough) of Manhattan. HINT: It will be the largest polygon in the borough.
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)
