Data, examples and exercises

   https://github.com/kjordahl/SciPy2013

Slides

   http://kjordahl.github.io/SciPy2013

The emphasis of this tutorial will be on:

• Open source tools for using geospatial data in python
• Python interfaces for major libraries
• "Pythonic" tools for working with data
• Interfacing with NumPy and the SciPy stack
• Hands-on exercises

Some big topics we will not be covering (much):

• Geospatial analysis
• Web mapping tools
• PostGIS (or other geospatial databases)
• Full GIS applications (e.g. ArcGIS and QGIS)

• Map projections with basemap
• GDAL/OGR python bindings
• Reading vector data with Fiona
• Geometry with Shapely
• Advanced topics: PostGIS, interfacing with GIS

Examples:

• EPSG:4326 latitude, longitude in WGS-84 coordinate system
• EPSG:900913 and EPSG:3857: Google spherical Mercator
• ESRI:102718: NAD 1983 StatePlane New York Long Island FIPS 3104 Feet

Create an SRS with pyproj:

>>> from pyproj import Proj
>>> p = Proj(init='epsg:3857')
>>> p.srs
'+units=m +init=epsg:3857 '
>>> p(-97.75, 30.25)
(-10881480.225042492, 3535725.659799159)
>>> p(-10881480.225042492, 3535725.659799159, inverse=True)
(-97.75, 30.25)

Try to do better than this one, too:

Cartopy is a Python package designed to make drawing maps for data analysis and visualisation as easy as possible.

• object oriented projection definitions
• point, line, polygon and image transformations between projections
• integration to expose advanced mapping in matplotlib
• integrated shapefile reading and Shapely

Vector data includes points, lines, polygons

     

Raster data includes images, digital elevation models, 2-D fields

     

GDAL (Geospatial Data Abstraction Library) is the open source Swiss Army knife of raster formats. It also includes the OGR simple features library for vector formats.

     

GDAL's python bindings expose most of the functionality of GDAL.

>>> from osgeo import gdal
>>> from osgeo import ogr

>>> for i in range(gdal.GetDriverCount()):
...    print gdal.GetDriver(i).LongName
Virtual Raster
GeoTIFF
National Imagery Transmission Format
Raster Product Format TOC format
ECRG TOC format
Erdas Imagine Images (.img)
CEOS SAR Image
CEOS Image
JAXA PALSAR Product Reader (Level 1.1/1.5)
Ground-based SAR Applications Testbed File Format (.gff)
ELAS
Arc/Info Binary Grid
Arc/Info ASCII Grid
GRASS ASCII Grid
...

Most of the GDAL bindings are straightforward wrappers of their C++ counterparts. One python-specific method is ReadAsArray().

img = 'filename.img'
geo = gdal.Open(img)
arr = geo.ReadAsArray()
imshow(arr, cmap=cm.BrBG_r)

May require further processing of the array data, but now it is just a NumPy array.

See examples/kauai.py

Python bindings do not raise exceptions for common errors. For example:

>>> from osgeo import gdal
>>> f = gdal.Open('foo.img')
>>> print f
None

You can enable exceptions explicitly with UseExceptions():

>>> gdal.UseExceptions()
>>> from osgeo import gdal
>>> f = gdal.Open('foo.img')
...
RuntimeError: `foo.img' does not exist in the file system,
and is not recognised as a supported dataset name.

It's hard to beat ogr2ogr for translating from one vector file format to another.

$ ogr2ogr -f geojson coutries.geojson ne_50m_admin_0_countries.shp

Translate the Natural Earth vector data to GeoJSON.

>>> for i in range(ogr.GetDriverCount()):
...    print ogr.GetDriver(i).name
ESRI Shapefile
MapInfo File
UK .NTF
SDTS
TIGER
S57
DGN
VRT
REC
Memory
...

see examples/ogr_read.py for a full example of reading a vector file with OGR in python.

exercises/geology.py

Shapely is a python library for geometric operations using the GEOS library.

Shapely can perform:

• geometry validation
• geometry creation (e.g. collections)
• geometry operations

>>> line = LineString([(0, 0), (1, 1), (0, 2), (2, 2), (3, 1), (1, 0)])
>>> dilated = line.buffer(0.5)
>>> eroded = dilated.buffer(-0.3)

object.almost_equals(other[, decimal=6])

object.contains(other)

object.crosses(other)

object.disjoint(other)

object.equals(other)

object.intersects(other)

object.touches(other)

object.within(other)

Zip is your friend! Use it with tuple unpacking to change between sequences of (x, y) pairs and seperate x and y sequences.

>>> pts = [(0, 0), (1, 0), (1, 1), (2, 1), (2, 2)]
>>> x, y = zip(*pts)
>>> print x, y
(0, 1, 1, 2, 2) (0, 0, 1, 1, 2)

Also, instead of calling f(x, y), you can just use

>>> f(*zip(*pts))

And when size matters,

>>> from itertools import izip

     

PostGIS is a spatial database extender for PostgreSQL that adds support for geographic objects allowing location queries to be run in SQL.

# select ST_NPoints(the_geom), rocktype1 FROM gis_schema.njgeol;
 st_npoints |         rocktype1
------------+----------------------------
        698 | limestone
        763 | siltstone
...

connect to PostGIS database using the python DB API

import psycopg2 as db
conn = db.connect("host=localhost dbname=gis_test user=kels")
cursor = conn.cursor()
sql = 'SELECT ST_NPoints(the_geom), rocktype1 FROM gis_schema.njgeol'
cursor.execute(sql)
results = cursor.fetchall()
conn.commit()
conn.close()

for n, rocktype in results:
    print n, 'points in a polygon of', rocktype

GEOMETRY

• Cartesian
• Uses projected units for calculations

GEOGRAPHY

• Spherical (Ellipsoidal)
• Distance and area calulations in meters

Three ways that QGIS can interact with python:

• Python console in QGIS application
• Python plugins
• Use PyQGIS in custom application

ArcPy is an official scripting layer for ESRI ArcGIS

• Mapnik
• GeoDjango
• GeoAlchemy
• Sean Gillies' blog
• Exploring climate data notebook