Introduction to Python for geospatial uses

Roger Veciana i Rovira

Twitter: @rveciana

Blog: http://geoexamples.blogspot.com

Python: What is it

An interpreted language
Dynamic types / Strong types
Object oriented
Multi-purpose
Thought to be fast (to learn, to use, to understand)
Thought to create a clean an uniform code

Python: Some advantages

Allows to address simple problems in a simple way, but also complex problems
- Script like coding
- Object oriented programming with exceptions, etc
Big standard library: "batteries included"
Lots of libraries and frameworksthat extend it
- Web creation
- System administration
- Scientific and geographic libraries
- Parallel programming
Python glues libraries coded in different languages.
Software that allows using it interactively: ipython, SAGE

Python basic instructions

a = 3
b = 4
a + b

Conditions

if b > a:
    print "b is bigger than a"

Iterate

for i in range(10):
    print "Iteration number " + str(i)

NumPy for matrix calculations

Open the ipython console

import numpy
a = numpy.ones((20,10))
a * 23
a + 1
a.sum()
a[3:6,4:8]
a[3:6,4:8] = 2
a[[2,3],[4,4]]

Try

Other functions: max(), mean() (try to pass axis=1 as an argument), reshape(), numpy.gradient(a), numpy.arange()
Otros array indexation: a[b], a[b,c]
Logical operations: a>2

GDAL/OGR

Geospatial Data Abstraction Library

Abstract data model for reading raster and vectorial GIS formats
More than 150 supported formats:
- GeoTIFF
- ENVI
- Shapefile
- PostGIS
- GRIB
- NetCDF
- VAISALA IRIS
Used by many applications and written in C, has Python bindings
Reading/Writting/Coordinate changes/GIS Operations

GDAL usage

The examples file is gdal_examples.py

GDAL exercices

Adding two arrays
Zones with more than sis liters in one hour
Zones without rain at 00 and rain at 01 hours
Total number of raining hours
Zones where has rained all the day

Fiona usage example

The file with the examples is fiona_examples.py

Fiona exercices

Find the distance between two points with the distance() function
If a linestring is created as: {'type': 'LineString', 'coordinates': ((1, 40),(2,41))}, use the function crosses and intersection between two linestrings
Calculate hoe many lightnings have fallen in each comarca (use comarcas.shp and llamps.shp)

Projections usage

The file with the examples is osr_examples.py

Projections WKT format

'PROJCS["ETRS89 / UTM zone 31N", Projection name (arbitrary)

GEOGCS["ETRS89", Datum code
DATUM["European_Terrestrial_Reference_System_1989", Datum name
SPHEROID["GRS 1980",6378137,298.2572221010002, AUTHORITY["EPSG","7019"]], Ellipsoid definition

TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6258"]], Aproximate transformatino toWGS84. 0 indicates that PROJ assumes that the projections are similar
PRIMEM["Greenwich",0],Reference meridian
UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4258"]],Transformation to radians

PROJECTION["Transverse_Mercator"],Label defining the projection that proj4 will use

PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",3],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0], UTM projection parameters
UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","25831"]]' Used units

Geotransform

(258047, 1000, 0.0, 4753640, 0.0, -1000)

Xgeo = gt(0) + Xpixel*gt(1) + Yline*gt(2)
Ygeo = gt(3) + Xpixel*gt(4) + Yline*gt(5)

geotransform = ds.GetGeoTransform()
gdal.ApplyGeoTransform(geotransform, x, y)
gdal.InvGeoTransform(geotransform)

Mayavi2 usage

The examples file is mayavi_examples.py