Introduction to Python for geospatial uses

SIGTE

Roger Veciana i Rovira

meteocat

Twitter: @rveciana

Blog: http://geoexamples.blogspot.com

Python: What is it

Python: Some advantages

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

GDAL/OGR

Geospatial Data Abstraction Library

GDAL usage

The examples file is gdal_examples.py

GDAL exercices

Fiona usage example

The file with the examples is fiona_examples.py

Fiona exercices

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

Some links

Thank you

/

#