7 min read

Python Geospatial Development

If you want to follow through the examples in this article, make sure you have the following Python libraries installed on your computer:

Reading and writing geo-spatial data

In this section, we will look at some examples of tasks you might want to perform that involve reading and writing geo-spatial data in both vector and raster format.

Task: Calculate the bounding box for each country in the world

In this slightly contrived example, we will make use of a Shapefile to calculate the minimum and maximum latitude/longitude values for each country in the world. This “bounding box” can be used, among other things, to generate a map of a particular country. For example, the bounding box for Turkey would look like this:

Python Geospatial Development

Start by downloading the World Borders Dataset from:

 

http://thematicmapping.org/downloads/world_borders.php

Decompress the .zip archive and place the various files that make up the Shapefile (the .dbf, .prj, .shp, and .shx files) together in a suitable directory.

We next need to create a Python program that can read the borders of each country. Fortunately, using OGR to read through the contents of a Shapefile is trivial:

import osgeo.ogr
shapefile = osgeo.ogr.Open(“TM_WORLD_BORDERS-0.3.shp”)
layer = shapefile.GetLayer(0)
for i in range(layer.GetFeatureCount()):
feature = layer.GetFeature(i)



The feature consists of a geometry and a set of fields. For this data, the geometry is a polygon that defines the outline of the country, while the fields contain various pieces of information about the country. According to the Readme.txt file, the fields in this Shapefile include the ISO-3166 three-letter code for the country (in a field named ISO3) as well as the name for the country (in a field named NAME). This allows us to obtain the country code and name like this:


countryCode = feature.GetField(“ISO3”)

countryName = feature.GetField(“NAME”)




We can also obtain the country’s border polygon using:


geometry = feature.GetGeometryRef()




There are all sorts of things we can do with this geometry, but in this case we want to obtain the bounding box or envelope for the polygon:


minLong,maxLong,minLat,maxLat = geometry.GetEnvelope()




Let’s put all this together into a complete working program:


# calcBoundingBoxes.py

import osgeo.ogr

shapefile = osgeo.ogr.Open(“TM_WORLD_BORDERS-0.3.shp”)

layer = shapefile.GetLayer(0)

countries = [] # List of (code,name,minLat,maxLat,

# minLong,maxLong) tuples.

for i in range(layer.GetFeatureCount()):

feature = layer.GetFeature(i)

countryCode = feature.GetField(“ISO3”)

countryName = feature.GetField(“NAME”)

geometry = feature.GetGeometryRef()

minLong,maxLong,minLat,maxLat = geometry.GetEnvelope()

countries.append((countryName, countryCode,

minLat, maxLat, minLong, maxLong))

countries.sort()

for name,code,minLat,maxLat,minLong,maxLong in countries:

print “%s (%s) lat=%0.4f..%0.4f, long=%0.4f..%0.4f”

% (name, code,minLat, maxLat,minLong, maxLong)




Running this program produces the following output:


% python calcBoundingBoxes.py

Afghanistan (AFG) lat=29.4061..38.4721, long=60.5042..74.9157

Albania (ALB) lat=39.6447..42.6619, long=19.2825..21.0542

Algeria (DZA) lat=18.9764..37.0914, long=-8.6672..11.9865

…




Task: Save the country bounding boxes into a Shapefile


While the previous example simply printed out the latitude and longitude values, it might be more useful to draw the bounding boxes onto a map. To do this, we have to convert the bounding boxes into polygons, and save these polygons into a Shapefile.


Creating a Shapefile involves the following steps:


    

  1. Define the spatial reference used by the Shapefile’s data. In this case, we’ll use the WGS84 datum and unprojected geographic coordinates (that is, latitude and longitude values). This is how you would define this spatial reference using OGR:

    import osgeo.osr
    
spatialReference = osgeo.osr.SpatialReference()
    
spatialReference.SetWellKnownGeogCS(‘WGS84’)
    

    

    We can now create the Shapefile itself using this spatial reference:
    

    import osgeo.ogr
    
driver = osgeo.ogr.GetDriverByName(“ESRI Shapefile”)
    
dstFile = driver.CreateDataSource(“boundingBoxes.shp”))
    
dstLayer = dstFile.CreateLayer(“layer”, spatialReference)
    

    

    2. 

  2. After creating the Shapefile, you next define the various fields that will hold the metadata for each feature. In this case, let’s add two fields to store the country name and its ISO-3166 code:

    fieldDef = osgeo.ogr.FieldDefn(“COUNTRY”, osgeo.ogr.OFTString)
    
fieldDef.SetWidth(50)
    
dstLayer.CreateField(fieldDef)
    

    fieldDef = osgeo.ogr.FieldDefn(“CODE”, osgeo.ogr.OFTString)
    
fieldDef.SetWidth(3)
    
dstLayer.CreateField(fieldDef)
    

    

    3. 

  3. We now need to create the geometry for each feature—in this case, a polygon defining the country’s bounding box. A polygon consists of one or more linear rings; the first linear ring defines the exterior of the polygon, while additional rings define “holes” inside the polygon. In this case, we want a simple polygon with a square exterior and no holes:

    linearRing = osgeo.ogr.Geometry(osgeo.ogr.wkbLinearRing)
    
linearRing.AddPoint(minLong, minLat)
    
linearRing.AddPoint(maxLong, minLat)
    
linearRing.AddPoint(maxLong, maxLat)
    
linearRing.AddPoint(minLong, maxLat)
    
linearRing.AddPoint(minLong, minLat)
    

    polygon = osgeo.ogr.Geometry(osgeo.ogr.wkbPolygon)
    
polygon.AddGeometry(linearRing)
    

    

    You may have noticed that the coordinate (minLong, minLat)was added to the linear ring twice. This is because we are defining line segments rather than just points—the first call to AddPoint()defines the starting point, and each subsequent call to AddPoint()adds a new line segment to the linear ring. In this case, we start in the lower-left corner and move counter-clockwise around the bounding box until we reach the lower-left corner again:
    

    Python Geospatial Development
    

     
    

    Once we have the polygon, we can use it to create a feature:
    

    feature = osgeo.ogr.Feature(dstLayer.GetLayerDefn())
    
feature.SetGeometry(polygon)
    
feature.SetField(“COUNTRY”, countryName)
    
feature.SetField(“CODE”, countryCode)
    
dstLayer.CreateFeature(feature)
    
feature.Destroy()
    

    

    Notice how we use the setField() method to store the feature’s metadata. We also have to call the Destroy() method to close the feature once we have finished with it; this ensures that the feature is saved into the Shapefile.
    4. 

  4. Finally, we call the Destroy() method to close the output Shapefile:

    dstFile.Destroy()
    

    

    5. 

  5. Putting all this together, and combining it with the code from the previous recipe to calculate the bounding boxes for each country in the World Borders Dataset Shapefile, we end up with the following complete program:

    # boundingBoxesToShapefile.py
    
import os, os.path, shutil
    
import osgeo.ogr
    
import osgeo.osr
    
# Open the source shapefile.
    
srcFile = osgeo.ogr.Open(“TM_WORLD_BORDERS-0.3.shp”)
    
srcLayer = srcFile.GetLayer(0)
    
# Open the output shapefile.
    
if os.path.exists(“bounding-boxes”):
    
shutil.rmtree(“bounding-boxes”)
    
os.mkdir(“bounding-boxes”)
    
spatialReference = osgeo.osr.SpatialReference()
    
spatialReference.SetWellKnownGeogCS(‘WGS84’)
    
driver = osgeo.ogr.GetDriverByName(“ESRI Shapefile”)
    
dstPath = os.path.join(“bounding-boxes”, “boundingBoxes.shp”)
    
dstFile = driver.CreateDataSource(dstPath)
    
dstLayer = dstFile.CreateLayer(“layer”, spatialReference)
    
fieldDef = osgeo.ogr.FieldDefn(“COUNTRY”, osgeo.ogr.OFTString)
    
fieldDef.SetWidth(50)
    
dstLayer.CreateField(fieldDef)
    
fieldDef = osgeo.ogr.FieldDefn(“CODE”, osgeo.ogr.OFTString)
    
fieldDef.SetWidth(3)
    
dstLayer.CreateField(fieldDef)
    
# Read the country features from the source shapefile.
    
for i in range(srcLayer.GetFeatureCount()):
    
feature = srcLayer.GetFeature(i)
    
countryCode = feature.GetField(“ISO3”)
    
countryName = feature.GetField(“NAME”)
    
geometry = feature.GetGeometryRef()
    
minLong,maxLong,minLat,maxLat = geometry.GetEnvelope()
    
# Save the bounding box as a feature in the output
    
# shapefile.
    
linearRing = osgeo.ogr.Geometry(osgeo.ogr.wkbLinearRing)
    
linearRing.AddPoint(minLong, minLat)
    
linearRing.AddPoint(maxLong, minLat)
    
linearRing.AddPoint(maxLong, maxLat)
    
linearRing.AddPoint(minLong, maxLat)
    
linearRing.AddPoint(minLong, minLat)
    
polygon = osgeo.ogr.Geometry(osgeo.ogr.wkbPolygon)
    
polygon.AddGeometry(linearRing)
    
feature = osgeo.ogr.Feature(dstLayer.GetLayerDefn())
    
feature.SetGeometry(polygon)
    
feature.SetField(“COUNTRY”, countryName)
    
feature.SetField(“CODE”, countryCode)
    
dstLayer.CreateFeature(feature)
    
feature.Destroy()
    
# All done.
    
srcFile.Destroy()
    
dstFile.Destroy()
    

    

    6. 



The only unexpected twist in this program is the use of a sub-directory called bounding-boxes to store the output Shapefile. Because a Shapefile is actually made up of multiple files on disk (a .dbf file, a .prj file, a .shp file, and a .shx file), it is easier to place these together in a sub-directory. We use the Python Standard Library module shutil to delete the previous contents of this directory, and then os.mkdir() to create it again.


If you aren’t storing the TM_WORLD_BORDERS-0.3.shp Shapefile in the same directory as the script itself, you will need to add the directory where the Shapefile is stored to your osgeo.ogr.Open() call. You can also store the boundingBoxes.shp Shapefile in a different directory if you prefer, by changing the path where this Shapefile is created.


Running this program creates the bounding box Shapefile, which we can then draw onto a map. For example, here is the outline of Thailand along with a bounding box taken from the boundingBoxes.shp Shapefile:


Python Geospatial Development


 




 





 





 
Packt
  






LEAVE A REPLY 





Please enter your comment!





Please enter your name here