opengeodata.de

Matera - Day 3

2017-06-21

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

Session 1 - Python basics

#!/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 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

sys.path.append('/home/user/my_module')
GetFeatureCount() - returns feature count
GetSpatialRef().ExportToProj4() - returns a proj4 string
GetPointCount() - returns point count
#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

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