Book Image

Python Geospatial Analysis Essentials

By : Erik Westra
Book Image

Python Geospatial Analysis Essentials

By: Erik Westra

Overview of this book

Table of Contents (13 chapters)
Python Geospatial Analysis Essentials
Credits
About the Author
About the Reviewers
www.PacktPub.com
Preface
Index

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:

  • We can calculate the centroid, which is the center-most point in the geometry.

  • We can calculate the bounding box for the geometry. This is a rectangle defining the northern, southern, eastern, and western edges of the polygon.

  • We can calculate the intersection between two geometries.

  • We can calculate the difference between two geometries.

    Note

    We could also calculate values such as the length and area of each polygon. However, because the World Borders Dataset uses what are called unprojected coordinates, the resulting length and area values would be measured in degrees rather than meters or miles. This means that the calculated lengths and areas wouldn't be very useful. We will look at the nature of map projections in the following chapter and find a way to get around this problem so we can calculate meaningful length and area values for polygons. But that's too complex for us to tackle right now.

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.