Solution Excercise 9 Team Hadochi

Jorn van der Ent

Michiel Voermans

23 January 2017

Load Modules and check for presence of ESRI Shapefile drive



In [ ]:

    
from osgeo import ogr
from osgeo import osr
import os

driverName = "ESRI Shapefile"
drv = ogr.GetDriverByName( driverName )
if drv is None:
    print "%s driver not available.\n" % driverName
else:
    print  "%s driver IS available.\n" % driverName

Set working directory to 'data'



In [ ]:

    
os.chdir('./data')

Interactive input system



In [ ]:

    
layername = raw_input("Name of Layer: ")

pointnumber = raw_input("How many points do you want to insert? ")

pointcoordinates = []
for number in range(1, (int(pointnumber)+1)):
    x = raw_input(("What is the Latitude (WGS 84) of Point %s ? " % str(number)))
    y = raw_input(("What is the Longitude (WGS 84) of Point %s ? " % str(number)))
    pointcoordinates += [(float(x), float(y))]

# e.g.:
# pointcoordinates =[(4.897070, 52.377956), (5.104480, 52.092876)]

Create shape file from input



In [ ]:

    
# Set filename
fn = layername + ".shp"

ds = drv.CreateDataSource(fn)

# Set spatial reference
spatialReference = osr.SpatialReference()
spatialReference.ImportFromEPSG(4326)

## Create Layer
layer=ds.CreateLayer(layername, spatialReference, ogr.wkbPoint)

# Get layer Definition
layerDefinition = layer.GetLayerDefn()

for pointcoord in pointcoordinates:
    ## Create a point
    point = ogr.Geometry(ogr.wkbPoint)

    ## SetPoint(self, int point, double x, double y, double z = 0)
    point.SetPoint(0, pointcoord[0], pointcoord[1]) 

    ## Feature is defined from properties of the layer:e.g:
    feature = ogr.Feature(layerDefinition)

    ## Lets add the points to the feature
    feature.SetGeometry(point)

    ## Lets store the feature in a layer
    layer.CreateFeature(feature)

ds.Destroy()

Convert shapefile to KML with bash



In [ ]:

    
bashcommand = 'ogr2ogr -f KML -t_srs crs:84 points.kml points.shp'
os.system(bashcommand)