Analyzing the data

In the previous section, we obtained an osgeo.ogr.Geometry object representing each country's outline. While there are a number of things we can do with this geometry object directly, in this case we'll take the outline and copy it into Shapely so that we can take advantage of Shapely's geospatial analysis capabilities. To do this, we have to export the geometry object out of OGR and import it as a Shapely object. For this, we'll use the WKT format. Still in the Python interpreter, let's grab a single feature's geometry and convert it into a Shapely object:

>>> import shapely.wkt
>>> feature = layer.GetFeature(0)
>>> geometry = feature.GetGeometryRef()
>>> wkt = geometry.ExportToWkt()
>>> outline = shapely.wkt.loads(wkt)

Because we loaded feature number 0, we retrieved the outline for Antigua and Barbuda, which would look like the following if we displayed it on a map:

The outline variable holds the outline of this country in the form of a Shapely MultiPolygon object. We can now use this object to analyze the geometry. Here are a few useful things we can do with a Shapely geometry:

Let's display the latitude and longitude for our feature's centroid:

>>> print outline.centroid.x, outline.centroid.y
-61.791127517 17.2801365868

Because Shapely doesn't know which coordinate system the polygon is in, it uses the more generic x and y attributes for a point, rather than talking about latitude and longitude values. Remember that latitude corresponds to a position in the north-south direction, which is the y value, while longitude is a position in the east-west direction, which is the x value.

We can also display the outline's bounding box:

>>> print outline.bounds
(-61.891113, 16.989719, -61.666389, 17.724998)

In this case, the returned values are the minimum longitude and latitude and the maximum longitude and latitude (that is, min_x, min_y, max_x, max_y).

There's a lot more we can do with Shapely, of course, but this is enough to prove that the Shapely library is working, and that we can read geospatial data from a shapefile and convert it into a Shapely geometry object for analysis.

This is as far as we want to go with using the Python shell directly—the shell is great for quick experiments like this, but it quickly gets tedious having to retype lines (or use the command history) when you make a typo. For anything more serious, you will want to write a Python program. In the final section of this chapter, we'll do exactly that: create a Python program that builds on what we have learned to solve a useful geospatial analysis problem.