Geospatial Development By Example with Python

For this book, we suggest using Python 2.7; this version of Python is fully compatible with the libraries and packages that we will use in the examples and also has precompiled binary files available on the Internet for Windows users. We will keep all the examples as compatible as possible with Python 3.4 so that it would be easy to port future upgrades.

Windows users may find compatibility problems with the 64-bit packages, so we recommend the 32-bit version of Python for them.

For Linux users, we will show the installation procedures for Ubuntu Linux distribution and use package managers, so you don't need to worry about versions and requirements; the package managers will do this for you.

The libraries that you will install are written in Python and other languages, the most common being C and C++. These libraries can abstract classes, methods, and functions to Python objects or have an extra layer that makes the connection; when this happens, we say that the library has Python bindings.

sudo apt-get install python

sudo apt-get install python-pip

A Python package is a directory containing one or more Python files (that is, modules) plus one __init__.py file (this can be just an empty file). This file tells Python Interpreter that the directory is a package.

When writing Python code, we can import packages and modules and use them in our code. The Python community does this a lot; many packages use other packages and so on, forming an intricate network of requirements and dependencies.

In order to facilitate the installation of packages and all the requirements for it to run, Python has a package manager called pip.

Pip looks for packages in a central repository (or on a user-defined place), downloads it, then downloads its dependencies, and installs them. Some packages also use libraries in other languages, such as C. In these cases, these libraries need to be compiled during the installation. Ubuntu users don't have problem with this because many compilers are already installed on the system, but this won't work on Windows by default.

In this topic, we will go through the installation process of every package used in the book.

NumPy is a package for scientific computing with Python. It handles multidimensional arrays of operations in a very efficient way. NumPy is required by OpenCV to run and will be used by many raster operations that we will perform in the examples. NumPy is also an efficient data container and will be our tool to calculate massive image data.

GDAL (Geospatial Data Abstraction Library) is composed of two packages that come together: OGR handles geospatial vector file formats, including coordinate system transformations and vector operations. GDAL is the raster part of the library, and in version 1.11, it comes packed with 139 drivers that can read, and some even create rasters. GDAL also comes packed with functions for raster transformations and calculations such as resizing, clipping, reprojecting, and so on.

In the following tables, there's an excerpt of the list of GDAL and OGR drivers with the most common formats that you may find:

The following table describes the OGR drivers:

Mapnik is a map rendering package. It is a free toolkit to develop mapping applications. It produces high-quality maps and is used on many applications, including OpenStreetMaps.

Shapely is a package for the manipulation and analysis of two dimensional geometries. It can perform operations such as union and subtraction of geometries. It also can perform tests and comparisons, such as when a geometry intersects other geometries.

Some packages do not require compilation steps. For Windows users, these are easier to install because they can be obtained and installed directly with pip with a single command.

c:\Python27\scripts\pip install django tabulate requests xmltodict psycopg2

sudo pip install django tabulate requests xmltodict psycopg2

Collecting django
  Downloading Django-1.9-py2.py3-none-any.whl (6.6MB)
    100% |################################| 6.6MB 43kB/s
Installing collected packages: django
Successfully installed django-1.9

IDEs are fancy text editors with tools and inspections regarding programming languages. You can surely use any text editor or IDE of your preference; none of the tasks in this book depends on the IDE, but an IDE will facilitate our work a lot because the suggested configuration will help you avoid mistakes and save time on typing, running, and debugging your code. The IDE checks the code for you and detects underlying errors; it even guesses what you are typing and completes the statements for you, runs the code with a simple command, and if there are exceptions, it provides links to the place where the exception occurred. For Windows or Linux, go to http://www.jetbrains.com/pycharm/ and click on the big orange button Get Pycharm Now. On the next page, select the free community edition.

Perform the following steps:

Now that we have all we need installed, we will go through our first example. In this example, we will test the installation and then see a glimpse of OGR's capabilities.

To do this, we will open a vector file containing the boundaries of all the countries in the world and make a list of country names. The objective of this simple example is to present the logic behind OGR objects and functions and give an understanding of how geospatial files are represented. Here's how:

Congratulations! You successfully opened a Shapefile and counted the number of features inside it. Now, let's understand what this code does.

The first line imports the ogr package. From this point on, all the functions are available as ogr.FunctionName(). Note that ogr doesn't follow the Python naming conventions for functions.

The line after the comment opens the OGR datasource (this opens the shapefile containing the data) and assigns the object to the datasource variable. Note that the path, even on Windows, uses a forward slash (/) and not a backslash.

The next line gets the first layer of the data source by its index (0). Some data sources can have many layers, but this is not the case of a Shapefile, which has only one layer. So, when working with Shapefiles, we always know that the layer of interest is layer 0.

In the final line, the print statement prints the number of features returned by layer.GetFeatureCount(). Here, we will use Python's string formatting, where the curly braces are replaced by the argument passed to format().

Now, perform the following steps:

The layer is iterable, but first, we need to ensure that we are at the beginning of the layer list with layer.ResetReading() (this is one of OGR's "gotcha" points).

The feature.GetFieldAsString(4) method returns the value of field 4 as a Python string. There are two ways of knowing whether the country names are in field 4:

Your complete code should look similar to the following:

import ogr


# Open the shapefile and get the first layer.
datasource = ogr.Open("../data/world_borders_simple.shp")
layer = datasource.GetLayerByIndex(0)
print("Number of features: {}".format(layer.GetFeatureCount()))

# Inspect the fields available in the layer.
feature_definition = layer.GetLayerDefn()
for field_index in range(feature_definition.GetFieldCount()):
    field_definition = feature_definition.GetFieldDefn(field_index)
    print("\t{}\t{}\t{}".format(field_index,
                                field_definition.GetTypeName(),
                                field_definition.GetName()))

# Print a list of country names.
layer.ResetReading()
for feature in layer:
    print(feature.GetFieldAsString(4)) 

Now, the objective is to know how much area is occupied by each country. However, the coordinates of country borders are expressed in latitude and longitude, and we can't calculate areas in this coordinate system. We want the area to be in the metric system, so first we need to convert the spatial reference system of the geometries.

Let's also take a step further in the programming techniques and start using functions to avoid the repetition of code. Perform the following steps:

That's wonderful! You just created a piece of very useful and reusable code. You now have a function that can open any shapefile, print the number of features, and return the ogr datasource. From now on, you can reuse this function in any of your projects.

You are already familiar with how this code works, but there are a few novelties here that deserve an explanation. The def statement defines a function with the def function_name(arguments): syntax.

Remember when I told you that OGR doesn't follow Python's naming convention? Well, the convention is that function names should all be in lowercase with an underscore between words. A good hint for names is to follow the verb_noun rule.

Right after the function's definition, you can see a description between triple quotes; this is a docstring, and it is used to document the code. It's optional but very useful for you to know what the function does.

Now, let's get back to our code. The second important thing to point out is the return statement. This makes the function return the values of the variables listed after the statement—in this case, the datasource.

Now, the next step is to create a function that performs the transformation. In OGR, the transformation is made in the feature's geometry, so we need to iterate over the features, get the geometry, and transform its coordinates. We will do this using the following steps:

Very well, unity conversion is another very important procedure in geoprocessing, and you just implemented it in your calculate_areas function.

However, having a list of numbers as the output is not very useful to us. So, it's time to combine everything that we did so far in order to extract valuable information.

You programmed three functions so far; now, let's add another one to our list by converting the code that generated a list of country names to a function and add this function to world_areas.py, as follows:

def get_country_names(datasource):
    """Returns a list of country names."""
    layer = datasource.GetLayerByIndex(0)
    country_names = []
    layer.ResetReading()
    for feature in layer:
        country_names.append(feature.GetFieldAsString(4))
    return country_names

Now, we have four functions, which are:

All these functions return iterables, with each item sharing the same index on all of them, thus making it easy to combine the information.

So, let's take advantage of this feature to sort the countries by area size and return a list of the five biggest countries and their areas. For this, add another function, as follows:

def get_biggest_countries(countries, areas, elements=5):
    """Returns a list of n countries sorted by area size."""    
    countries_list = [list(country) 
                      for country in zip(areas, countries)]

    sorted_countries = sorted(countries_list, 
                              key=itemgetter(0), reverse=True) 
    return sorted_countries[:5]

In the first line, the two lists are zipped together, producing a list of country-area pairs. Then, we used the Python list's sorted method, but as we don't want the lists to be sorted by both values, we will define the key for sorting. Finally, the list is sliced, returning only the desired number of values.

In this chapter, we had a brief introduction to the libraries and packages that we will use in this book. By installing these libraries, you also learned the general procedure of how to search and install Python packages. You can use this procedure in other cases whenever you feel the need for other libraries in your applications.

Then, we wrote code that made use of the OGR library to open a shapefile and perform area calculation and sorting. These simple procedures showed a little bit of the internal organization of OGR, how it handles geographic data, and how it is possible to extract information from them. In the next chapter, we will use some of the techniques learned here to read data and process vector points.