GIS Project in Python from OSGeo4W - Part 2

Here are the working codes of the examples provided in Chapter 2 of of Python Geospatial Development Second Edition. I have tried to stick to the text book version of the programs and to the links provided in my previous post

Example 1 and Example 2

import os
from osgeo import ogr

daShapefile = r"tl_2012_us_state.shp"
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open(daShapefile, 0) # 0 means read-only. 1 means writeable.

# Check to see if shapefile is found.
if dataSource is None:
    print 'Could not open %s' % (daShapefile)
    print 'The driver is %s' % driver
else:
    print 'Opened %s and the driver used is %s' % ((daShapefile), driver)
    layer = dataSource.GetLayer()
    numLayers = dataSource.GetLayerCount()
    featureCount = layer.GetFeatureCount()
    print "Number of Layers in %s: %d and had %d features" % (os.path.basename(daShapefile),numLayers, featureCount)

    for featureNum in range(featureCount) :
        feature = layer.GetFeature(featureNum)
        featureName = feature.GetField("NAME")
        geometry = feature.GetGeometryRef()
        geometryName = geometry.GetGeometryName()
       
        print "Feature %d . %s with a %s in geometry data and has the following attributes: " % (featureNum, featureName, geometryName)

        attributes = feature.items()

        for key, value in attributes.items():
            print " %s = %s \t" %(key, value)
           
Example 3:
import osgeo
import os
from osgeo import ogr

daShapefile = r"tl_2012_us_state.shp"
driver = ogr.GetDriverByName('ESRI Shapefile')
shapefile = driver.Open(daShapefile, 0)

def analyzeGeometry(geometry, indent=0) :
    s = []
    s.append(" " * indent)
    s.append(geometry.GetGeometryName())
    if geometry.GetPointCount() > 0 :
        s.append(" with %d data points " %geometry.GetPointCount())
    if geometry.GetGeometryCount() > 0 :
        s.append(" containing:")

    print "".join(s)

    for i in range(geometry.GetGeometryCount()):
           
        analyzeGeometry(geometry.GetGeometryRef(i), indent+1)

layer = shapefile.GetLayer(0)
feature = layer.GetFeature(55)
geometry = feature.GetGeometryRef()

analyzeGeometry(geometry)

Example 4:
import osgeo
import os
from osgeo import ogr

daShapefile = r"tl_2012_us_state.shp" # this line assumes that the program and data are in the same directory

# In case the python program and shapefile are in two different directories  then :
# Please provide the full directory path of the shapefile to be analysed

driver = ogr.GetDriverByName('ESRI Shapefile')
shapefile = driver.Open(daShapefile, 0)

def findPoints(geometry, results) :
    for i in range(geometry.GetPointCount()):
        x,y,z = geometry.GetPoint(i)
        if results['north'] == None or results['north'] [1] < y:
            results['north'] = (x,y)
        if results['south'] == None or results['south'] [1] > y:
            results['south'] = (x,y)
    for i in range(geometry.GetGeometryCount()):
        findPoints(geometry.GetGeometryRef(i), results)

layer = shapefile.GetLayer(0)
feature = layer.GetFeature(55)
geometry = feature.GetGeometryRef()

results = { 'north' : None,
            'south' : None}

findPoints(geometry, results)

print "Northern most point is (%0.4f, %0.4f)" % results['north']
print "Southern most point is (%0.4f, %0.4f)" % results['south']

Example 5 and 6:
import osgeo
import os
from osgeo import ogr
import math
from decimal import *

daShapefile = r"tl_2012_us_state.shp" # this line assumes that the program and data are in the same directory

# In case the python program and shapefile are in two different directories  then :
# Please provide the full directory path of the shapefile to be analysed

driver = ogr.GetDriverByName('ESRI Shapefile')
shapefile = driver.Open(daShapefile, 0)

def findPoints(geometry, results) :
    for i in range(geometry.GetPointCount()):
        x,y,z = geometry.GetPoint(i)
        if results['north'] == None or results['north'] [1] < y:
            results['north'] = (x,y)
        if results['south'] == None or results['south'] [1] > y:
            results['south'] = (x,y)
    for i in range(geometry.GetGeometryCount()):
        findPoints(geometry.GetGeometryRef(i), results)

layer = shapefile.GetLayer(0)
feature = layer.GetFeature(55)
geometry = feature.GetGeometryRef()

results = { 'north' : None,
            'south' : None}

findPoints(geometry, results)

print "Northern most point is (%0.4f, %0.4f)" % results['north']
print "Southern most point is (%0.4f, %0.4f)" % results['south']

for i in range(len(results)):
    print "%d %0.4f" %(i, results['north'] [i])

# Text book exercise: Calculate the mid point and distance between the north and south extreme points
# print "%s" % type(results) # check the data type of results

latlong = results.values()
latlongrad = [[0.0] *2 for i in range(2)] # initialising and creating a 2 dimensional array of floats

for i in range(2):
    for j in range(2):
        latlongrad[i][j] = math.radians(latlong[i][j])

dLat = latlongrad[0][1]-latlongrad[1][1]
dLong = latlongrad[0][0]-latlongrad[1][0]

# haversine formula:

a = math.sin(dLat/2)**2 + math.cos(latlongrad[0][1]) * math.cos(latlongrad[1][1]) * math.sin(dLong/2)**2

c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))

distance = 6371 * c

print " Great Circle distance between northern and southern most points is %0.0f kilometers" % distance