A basic introduction to Python, its core concepts as well as problem solving strategies based on certain geospatial packages and internet research.

  • instructor for the day: fpl

Session 1 - Python basics

  • ~90% of packages are ported to Python3 to this date; check if your a package you need is already ported: Py3Readiness
1
chmod a+x - executable for all groups
1
2
3
4
5
6
7
#!/usr/bin/python
#~ # This is the well-known Fibonacci series
a, b = 0, 1
while b < 2000:
print a
a, b = b, a + b
  • keyword arguments can be useful if you want pass some, but not all, arguments to a function which takes multiple parameters
1
2
3
4
5
6
7
8
9
10
11
12
13
'''
Keyword arguments in calling functions
'''
def fibonacci(n=2000):
a, b = 0, 1
f = []
while b<n:
f.append(a)
a, b = b, a+b
return f
s = fibonacci(n=10000)
print s

Session 2 - Python OGR

  • adding a custom path for modules:
1
sys.path.append('/home/user/my_module')
  • useful Python OGR functions
1
2
3
GetFeatureCount() - returns feature count
GetSpatialRef().ExportToProj4() - returns a proj4 string
GetPointCount() - returns point count
  • a simple script for checking out some OGR Python functions
    • attention (for me) to be paid to geometry.GetGeometryName() == 'POINT' for checking the Geometry type
    • also an interesting (& crude) way to get CLI arguments: args = [] /n args.append(sys.argv)
    • open a shape file shp = ogr.Open(shpFile)
    • get geometry of feature - geometry = feature.GetGeometryRef()
    • get field name - feature.GetField(fieldName)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#Examine a shapefile with ogr
from osgeo import ogr
import os
import sys
args = []
args.append(sys.argv)
# set working dir
os.chdir('../files')
# check if file name was given
try:
shpFile = args[0][1]
except:
print 'No input file specified.'
sys.exit(1)
# check if field was given
try:
fieldName = args[0][2]
except:
print 'No field specified.'
sys.exit(1)
# open the shapefile
shp = ogr.Open(shpFile)
# Get the layer
try:
layer = shp.GetLayer()
except:
print 'File not found.'
sys.exit(1)
# Loop through the features
# and print information about them
for feature in layer:
geometry = feature.GetGeometryRef()
# check if the field name exists
try:
feature.GetField(fieldName)
except:
print 'Wrong field name given.'
sys.exit(1)
if geometry.GetGeometryName() == 'POINT':
# print the info
print geometry.GetX(), geometry.GetY(), feature.GetField(fieldName)
else:
print 'Only works for point geometries.'
sys.exit(1)

Session 3 - Python OGR

  • get all functions for a object in OGR, open iypthon
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
from osgeo import ogr
shp = ogr.Open('point.shp')
shp. + Press TAB
> shp.CommitTransaction shp.GetDriver shp.GetMetadata_List shp.SetDescription
shp.CopyLayer shp.GetLayer shp.GetName shp.SetMetadata
shp.CreateLayer shp.GetLayerByIndex shp.GetRefCount shp.SetMetadataItem
shp.DeleteLayer shp.GetLayerByName shp.GetStyleTable shp.SetStyleTable
shp.Dereference shp.GetLayerCount shp.GetSummaryRefCount shp.StartTransaction
shp.Destroy shp.GetMetadata shp.Reference shp.SyncToDisk
shp.ExecuteSQL shp.GetMetadataDomainList shp.Release shp.TestCapability
shp.FlushCache shp.GetMetadataItem shp.ReleaseResultSet shp.name
shp.GetDescription shp.GetMetadata_Dict shp.RollbackTransaction shp.this
layer = shp.GetLayer()
layer. + Press TAB
>
layer.AlterFieldDefn layer.GetFeature layer.GetSpatialRef layer.SetNextByIndex
layer.Clip layer.GetFeatureCount layer.GetStyleTable layer.SetSpatialFilter
layer.CommitTransaction layer.GetFeaturesRead layer.Identity layer.SetSpatialFilterRect
layer.CreateFeature layer.GetGeomType layer.Intersection layer.SetStyleTable
layer.CreateField layer.GetGeometryColumn layer.Reference layer.StartTransaction
layer.CreateFields layer.GetLayerDefn layer.ReorderField layer.SymDifference
layer.CreateGeomField layer.GetMetadata layer.ReorderFields layer.SyncToDisk
layer.DeleteFeature layer.GetMetadataDomainList layer.ResetReading layer.TestCapability
layer.DeleteField layer.GetMetadataItem layer.RollbackTransaction layer.Union
layer.Dereference layer.GetMetadata_Dict layer.SetAttributeFilter layer.Update
layer.Erase layer.GetMetadata_List layer.SetDescription layer.next
layer.FindFieldIndex layer.GetName layer.SetFeature layer.schema
layer.GetDescription layer.GetNextFeature layer.SetIgnoredFields layer.this
layer.GetExtent layer.GetRefCount layer.SetMetadata
layer.GetFIDColumn layer.GetSpatialFilter layer.SetMetadataItem