OSGeo (GDAL/OGR) exercise Geoscripting.

In this assignment you’ll interpret the previously build ArcGIS model/script in a complete open source python environment that makes use of the OSGeo (GDAL/OGR) module.

For this assignment you need the OS-GEO linux vitual machine, which is installed on the lab PC’s. The installation includes a python version which has all modules needed pre-installed.

To load the OSGeo modules in python (http://gdal.org/python/):

In [2]:
import os
from osgeo import ogr
from osgeo import osr

As in the ArcPy assignment a point should be created first, which can next be used as the centre of a buffer operation.

In [3]:
wkt = "POINT (173914.00 441864.00)"
pt = ogr.CreateGeometryFromWkt(wkt)
print pt.GetX(), pt.GetY()
173914.0 441864.0

This will construct an in memory geometry of a point.

An empty dataset (shape-file) should be created where the constructed buffers can be stored, so they can also be used to visualize the results later on. A driver and file should be specified.

In [13]:
driver = ogr.GetDriverByName('ESRI Shapefile')
outputSf = 'test.shp'

# Remove output shapefile if it already exists
if os.path.exists(outputSf):
    driver.DeleteDataSource(outputSf)

# Create output shapefile 
datasource = driver.CreateDataSource(outputSf)

And set the coordinate reference system, the reference of the National Dutch Grid (rijksdriehoekstelsel) is known in the Geodetic Parameter Dataset as EPSG28992 (http://www.epsg.org/)

In [14]:
proj = osr.SpatialReference()
proj.ImportFromEPSG(28992)
Out[14]:
0

Use these parameters to construct a new layer for the to be constructed buffers

In [15]:
layer = datasource.CreateLayer('test', geom_type=ogr.wkbPolygon, srs = proj)
feature = ogr.Feature(layer.GetLayerDefn())

As all necessary conditions are set, and the final buffer operation can be performed. Features can be buffered with a given radius e.g. 1000 meters in size.

Search in the OSGeo documentation how to construct a buffer around the in memory point-geometry http://gdal.org/python/ .

Hint: You need the OGR module to perform vector based spatial operations. Make sure you work with geometries

In [16]:
bufferDistance = 5000
poly = pt.Buffer(bufferDistance)

Each new geometry of a constructed buffer should be added to the shape-file.

In [17]:
# Wkt of the constructed buffer polygon
polygon = ogr.CreateGeometryFromWkt(poly.ExportToWkt())
feature.SetGeometry(polygon)

# Add the features to the layer/shapefile
layer.CreateFeature(feature)
Out[17]:
0

Make it a habit to clean-up the memory after finishing loops/scripts:

In [18]:
# Clean-up the added buffer polygon for next loop
polygon.Destroy()

feature.Destroy()
datasource.Destroy()

The result of the assignment can be visualized in QGIS

The results can be visualized using Cartopy and Matplotlib

In [19]:
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature

import matplotlib.pyplot as plt
import matplotlib.cm as cm
%matplotlib inline 

Plot the buffered location on the OSM basemap

To be able to use epsg codes for projection settings, cartopy depends on pyepsg, this module is not installed by default.
I can be obtained from https://pypi.python.org/pypi/pyepsg/0.1.0.
Unpack and open the folder in the command shell to install:
python setup.py install

In [20]:
ax = plt.axes(projection=ccrs.epsg(28992))

extent = [(pt.GetX()-10000), (pt.GetX()+10000), (pt.GetY()-10000), (pt.GetY()+10000)]
#extent = [(pt.GetX()-1000), (pt.GetY()-1000), (pt.GetX()+1000), (pt.GetY()+1000)]
#extent = [51.5, 5.0, 52.0, 6.0]
print extent
#ax.set_extent(extent, crs=ccrs.PlateCarree())
ax.set_extent(extent, crs=ccrs.epsg(28992))

#Add the OSM basemap (this makes the drawing very slow)
bg = cimgt.OSM()
ax.add_image(bg, 12, alpha=0.9)

#Add the buffered object
#fname = 'test.shp'
shape_feature = ShapelyFeature(Reader(outputSf).geometries(),
                                ccrs.epsg(28992), facecolor='red')
ax.add_feature(shape_feature)
[163914.0, 183914.0, 431864.0, 451864.0]
Out[20]:
<cartopy.mpl.feature_artist.FeatureArtist at 0xaccbe08c>
In [21]:
plt.show()

Reproject the buffer object to lon/lat

In [22]:
from osgeo import ogr, osr
import os

driver = ogr.GetDriverByName('ESRI Shapefile')

# input SpatialReference
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(28992)

# output SpatialReference
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(4326)

# create the CoordinateTransformation
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

# get the input layer
inDataSet = driver.Open(r'test.shp')
inLayer = inDataSet.GetLayer()

# create the output layer
outputShapefile = r'test_4326.shp'
if os.path.exists(outputShapefile):
    driver.DeleteDataSource(outputShapefile)
outDataSet = driver.CreateDataSource(outputShapefile)
outLayer = outDataSet.CreateLayer("test_4326", geom_type=ogr.wkbMultiPolygon)

# add fields
inLayerDefn = inLayer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
    fieldDefn = inLayerDefn.GetFieldDefn(i)
    outLayer.CreateField(fieldDefn)

# get the output layer's feature definition
outLayerDefn = outLayer.GetLayerDefn()

# loop through the input features
inFeature = inLayer.GetNextFeature()
while inFeature:
    # get the input geometry
    geom = inFeature.GetGeometryRef()
    # reproject the geometry
    geom.Transform(coordTrans)
    # create a new feature
    outFeature = ogr.Feature(outLayerDefn)
    # set the geometry and attribute
    outFeature.SetGeometry(geom)
    for i in range(0, outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
    # add the feature to the shapefile
    outLayer.CreateFeature(outFeature)
    # destroy the features and get the next input feature
    outFeature.Destroy()
    inFeature.Destroy()
    inFeature = inLayer.GetNextFeature()

# close the shapefiles
inDataSet.Destroy()
outDataSet.Destroy()

Draw the reprojected buffer on a basemap

In [3]:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
In [13]:
%matplotlib inline 

fig = plt.figure(figsize=(6,4))

#Custom adjust of the subplots
plt.subplots_adjust(left=0.05,right=0.95,top=0.90,bottom=0.05,wspace=0.15,hspace=0.05)
ax = plt.subplot(111)

#Let's create a basemap around Netherlands & Belgium
m = Basemap(resolution='i',projection='merc', llcrnrlat=51.5,urcrnrlat=52.5,llcrnrlon=5.0,urcrnrlon=6.0,lat_ts=51.5)
m.drawcountries(linewidth=0.5,color='gray')
m.drawcoastlines(linewidth=0.5,color='gray')
m.drawmapboundary(fill_color='gray')

#m.drawparallels(np.arange(51.,5.,1.),labels=[1,0,0,0],color='gray',dashes=[3,1],linewidth=0.4) # draw parallels
#m.drawmeridians(np.arange(1.,9.,1.),labels=[0,0,0,1],color='gray',dashes=[3,1],linewidth=0.4) # draw meridians

#m.fillcontinents(color='black',lake_color='gray')

m.readshapefile('test_4326', 'buffer', drawbounds=True, linewidth=2, color='b')


plt.title('Benelux in Mercator Projection')
plt.show()