If you want to follow through the examples in this article, make sure you have the following Python libraries installed on your computer:
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.
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:
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
…
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:
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)
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)
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:
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.
dstFile.Destroy()
# 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()
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:
I remember deciding to pursue my first IT certification, the CompTIA A+. I had signed…
Key takeaways The transformer architecture has proved to be revolutionary in outperforming the classical RNN…
Once we learn how to deploy an Ubuntu server, how to manage users, and how…
Key-takeaways: Clean code isn’t just a nice thing to have or a luxury in software projects; it's a necessity. If we…
While developing a web application, or setting dynamic pages and meta tags we need to deal with…
Software architecture is one of the most discussed topics in the software industry today, and…