Data, examples and exercises

   https://github.com/kjordahl/SciPy-Tutorial-2015

Slides

   http://kjordahl.github.io/SciPy-Tutorial-2015

The emphasis of this tutorial will be on:

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

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

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

• Reading vector data with Fiona
• Geometry with Shapely
• Brief interlude on map projections
• Rasterio
• GeoPandas

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

GDAL python bindings are not very "pythonic"

geo = gdal.Open(raster_file)
drv = geo.GetDriver()
print(drv.GetMetadataItem('DMD_LONGNAME'))

Most of the GDAL bindings are thin wrappers of their C++ counterparts.

See examples/gdal_example.py

OGR → Fiona

GDAL → Rasterio

Usage: fio [OPTIONS] COMMAND [ARGS]...

Commands:
  cat      Concatenate and print the features of datasets
  collect  Collect a sequence of features.
  dump     Dump a dataset to GeoJSON.
  info     Print information about a dataset.
  insp     Open a dataset and start an interpreter.
  load     Load GeoJSON to a dataset in another format.

exercises/Fiona exercise.ipynb

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

exercises/Shapely exercise.ipynb

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)

examples/pyproj example.ipynb

with rasterio.open('manhattan.tif') as f:
    img = f.read(1)
imshow(img, cmap='gray')

$ rio insp --ipython manhattan.tif
Rasterio 0.24.1 Interactive Inspector (Python 2.7.9)
Type "src.meta", "src.read_band(1)", or "help(src)" for more information.
In [1]: src.crs
Out[1]: {'init': u'epsg:26918'},

In [2]: src.height, src.width
Out[2]: (5205, 5757)

In [3]: grd = src.read(1)

In [4]: grd.shape
Out[4]: (5205, 5757)

examples/warp radar.ipynb

see also rasterio docs

from rasterio.features import rasterize
mask = rasterize([poly], transform=src.transform, out_shape=src.shape)

examples/rasterio_mask.py

utm18n = ccrs.UTM(18)
ax = plt.axes(projection=utm18n)
ax.imshow(src.read(1), origin='upper',
          extent=(left, right, bottom, top), cmap='gray')
ax.coastlines(resolution='10m', linewidth=4, color='red')
gl = ax.gridlines(linewidth=2, color='lightblue', alpha=0.5, linestyle='--')

examples/rasterio_cartopy.py (see also Cartopy docs)

exercises/landsat exercise.ipynb

GeoPandas is pure python (2.6, 2.7, 3.3+)

GeoPandas depends on:

• Pandas (0.13 and up)
• Shapely (GEOS)
• Fiona (GDAL/OGR)
• Pyproj (PROJ.4)
• Matplotlib (and Descartes)
• psycopg2, sqlalchemy, geopy, rtree (optional)

• Geometry operations (Shapely)
• Data alignment (pandas)
• Coordinate transformations (pyproj)
• Read/write GIS file formats (Fiona)
• Create a GeoDataFrame from PostGIS table
• Output any object as geoJSON
• Plotting

examples/GeoDataFrame example.ipynb

see also:

• GeoSeries notebook
• Texas Natural Regions notebook
• Geocoding example
• CitiBike example

exercises/Idaho lakes.ipynb

joined = sjoin(points, boros, how='inner')

examples/sjoin.py

• FOSS4G talk by Sean Gillies
• PySAL
• PostGIS
• Leaflet
• MPLLeaflet
• d3.geo