Python Libraries for Geospatial Development

Exclusive offer: get 50% off this eBook here
Python Geospatial Development - Second Edition

Python Geospatial Development - Second Edition — Save 50%

Learn to build sophisticated mapping applications from scratch using Python tools for geospatial development with this book and ebook

$29.00    $14.50
by Erik Westra | June 2013 | Open Source

This article by Erik Westra the author of Python Geospatial Development - Second Edition, examines a number of libraries and other tools which can be used for geospatial development in Python.

More specifically, we will cover:

  • Python libraries for reading and writing geospatial data

  • Python libraries for dealing with map projections

  • Libraries for analyzing and manipulating geospatial data directly within your Python programs

  • Tools for visualizing geospatial data

Note that there are two types of geospatial tools which are not discussed in this article: geospatial databases and geospatial web toolkits. Both of these will be examined in detail later in this book.

(For more resources related to this topic, see here.)

Reading and writing geospatial data

While you could in theory write your own parser to read a particular geospatial data format, it is much easier to use an existing Python library to do this. We will look at two popular libraries for reading and writing geospatial data: GDAL and OGR.

GDAL/OGR

Unfortunately, the naming of these two libraries is rather confusing. Geospatial Data Abstraction Library ( GDAL), was originally just a library for working with raster geospatial data, while the separate OGR library was intended to work with vector data. However, the two libraries are now partially merged, and are generally downloaded and installed together under the combined name of "GDAL". To avoid confusion, we will call this combined library GDAL/OGR and use "GDAL" to refer to just the raster translation library.

A default installation of GDAL supports reading 116 different raster file formats, and writing to 58 different formats. OGR by default supports reading 56 different vector file formats, and writing to 30 formats. This makes GDAL/OGR one of the most powerful geospatial data translators available, and certainly the most useful freely-available library for reading and writing geospatial data.

GDAL design

GDAL uses the following data model for describing raster geospatial data:

Let's take a look at the various parts of this model:

  • A dataset holds all the raster data, in the form of a collection of raster "bands", along with information that is common to all these bands. A dataset normally represents the contents of a single file.

  • A raster band represents a band, channel, or layer within the image. For example, RGB image data would normally have separate bands for the red, green, and blue components of the image.

  • The raster size specifies the overall width and height of the image, in pixels.

  • The georeferencing transform converts from (x, y) raster coordinates into georeferenced coordinates—that is, coordinates on the surface of the earth. There are two types of georeferencing transforms supported by GDAL: affine transformations and ground control points.

    • An affine transformation is a mathematical formula allowing the following operations to be applied to the raster data:

      More than one of these operations can be applied at once; this allows you to perform sophisticated transforms such as rotations.

      Affine transformations are sometimes referred to as linear transformations.

    • Ground Control Points ( GCPs) relate one or more positions within the raster to their equivalent georeferenced coordinates, as shown in the following figure:

      Note that GDAL does not translate coordinates using GCPs— that is left up to the application, and generally involves complex mathematical functions to perform the transformation.

  • The coordinate system describes the georeferenced coordinates produced the georeferencing transform. The coordinate system includes the projection and datum, as well as the units and scale used by the raster data.

  • The metadata contains additional information about the dataset as a whole.

Each raster band contains the following (among other things):

  • The band raster size: This is the size (number of pixels across and number of lines high) for the data within the band. This may be the same as the raster size for the overall dataset, in which case the dataset is at full resolution, or the band's data may need to be scaled to match the dataset.

  • Some band metadata providing extra information specific to this band.

  • A color table describing how pixel values are translated into colors.

  • The raster data itself.

GDAL provides a number of drivers which allow you to read (and sometimes write) various types of raster geospatial data. When reading a file, GDAL selects a suitable driver automatically based on the type of data; when writing, you first select the driver and then tell the driver to create the new dataset you want to write to.

GDAL example code

A Digital Elevation Model ( DEM) file contains height values. In the following example program, we use GDAL to calculate the average of the height values contained in a sample DEM file. In this case, we use a DEM file downloaded from the GLOBE elevation dataset:

from osgeo import gdal,gdalconst import struct dataset = gdal.Open("data/e10g") band = dataset.GetRasterBand(1) fmt = "<" + ("h" * band.XSize) totHeight = 0 for y in range(band.YSize): scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, band.DataType) values = struct.unpack(fmt, scanline) for value in values: if value == -500: # Special height value for the sea -> ignore. continue totHeight = totHeight + value average = totHeight / (band.XSize * band.YSize) print "Average height =", average

As you can see, this program obtains the single raster band from the DEM file, and then reads through it one scanline at a time. We then use the struct standard Python library module to read the individual height values out of the scanline. Because the GLOBE dataset uses a special height value of -500 to represent the ocean, we exclude these values from our calculations. Finally, we use the remaining height values to calculate the average height, in meters, over the entire DEM data file.

OGR design

OGR uses the following model for working with vector-based geospatial data:

Let's take a look at this design in more detail:

  • The data source represents the file you are working with—though it doesn't have to be a file. It could just as easily be a URL or some other source of data.

  • The data source has one or more layers , representing sets of related data. For example, a single data source representing a country may contain a "terrain" layer, a "contour lines" layer, a "roads" later, and a "city boundaries" layer. Other data sources may consist of just one layer. Each layer has a spatial reference and a list of features.

  • The spatial reference specifies the projection and datum used by the layer's data.

  • A feature corresponds to some significant element within the layer. For example, a feature might represent a state, a city, a road, an island, and so on. Each feature has a list of attributes and a geometry.

  • The attributes provide additional meta-information about the feature. For example, an attribute might provide the name for a city's feature, its population, or the feature's unique ID used to retrieve additional information about the feature from an external database.

  • Finally, the geometry describes the physical shape or location of the feature. Geometries are recursive data structures that can themselves contain sub-geometries—for example, a "country" feature might consist of a geometry that encompasses several islands, each represented by a subgeometry within the main "country" geometry.

    The geometry design within OGR is based on the Open Geospatial Consortium's "Simple Features" model for representing geospatial geometries. For more information, see http://www.opengeospatial.org/standards/sfa .

Like GDAL, OGR also provides a number of drivers which allow you to read (and sometimes write) various types of vector-based geospatial data. When reading a file, OGR selects a suitable driver automatically; when writing, you first select the driver and then tell the driver to create the new data source to write to.

OGR example code

The following example program uses OGR to read through the contents of a shapefile, printing out the value of the NAME attribute for each feature along with the geometry type:

from osgeo import ogr shapefile = ogr.Open("TM_WORLD_BORDERS-0.3.shp") layer = shapefile.GetLayer(0) for i in range(layer.GetFeatureCount()): feature = layer.GetFeature(i) name = feature.GetField("NAME") geometry = feature.GetGeometryRef() print i, name, geometry.GetGeometryName()

Documentation

GDAL and OGR are well documented, but with a catch for Python programmers. The GDAL/OGR library and associated command-line tools are all written in C and C++. Bindings are available which allow access from a variety of other languages, including Python, but the documentation is all written for the C++ version of the libraries. This can make reading the documentation rather challenging—not only are all the method signatures written in C++, but the Python bindings have changed many of the method and class names to make them more "pythonic".

Fortunately, the Python libraries are largely self-documenting, thanks to all the docstrings embedded in the Python bindings themselves. This means you can explore the documentation using tools such as Python's built-in pydoc utility, which can be run from the command line like this:

% pydoc -g osgeo

This will open up a GUI window allowing you to read the documentation using a web browser. Alternatively, if you want to find out about a single method or class, you can use Python's built-in help() command from the Python command line, like this:

>>> import osgeo.ogr >>> help(osgeo.ogr.DataSource.CopyLayer)

Not all the methods are documented, so you may need to refer to the C++ docs on the GDAL website for more information, and some of the docstrings are copied directly from the C++ documentation—but in general the documentation for GDAL/OGR is excellent, and should allow you to quickly come up to speed using this library.

Availability

GDAL/OGR runs on modern Unix machines, including Linux and Mac OS X, as well as most versions of Microsoft Windows. The main website for GDAL can be found at:

http://gdal.org

The main website for OGR is at:

http://gdal.org/ogr

To download GDAL/OGR, follow the Downloads link on the main GDAL website. Windows users may find the FWTools package useful, as it provides a wide range of geospatial software for win32 machines, including GDAL/OGR and its Python bindings. FWTools can be found at:

http://fwtools.maptools.org

For those running Mac OS X, prebuilt binaries can be obtained from:

http://www.kyngchaos.com/software/frameworks

Make sure that you install GDAL Version 1.9 or later, as you will need this version to work through the examples in this book.

Being an open source package, the complete source code for GDAL/OGR is available from the website, so you can compile it yourself. Most people, however, will simply want to use a prebuilt binary version.

Dealing with projections

One of the challenges of working with geospatial data is that geodetic locations (points on the Earth's surface) are mapped into a two-dimensional Cartesian plane using a cartographic projection. Whenever you have some geospatial data, you need to know which projection that data uses. You also need to know the datum (model of the Earth's shape) assumed by the data.

A common challenge when dealing with geospatial data is that you have to convert data from one projection/datum to another. Fortunately, there is a Python library pyproj which makes this task easy.

pyproj

pyproj is a Python "wrapper" around another library called PROJ.4. "PROJ.4" is an abbreviation for Version 4 of the PROJ library. PROJ was originally written by the US Geological Survey for dealing with map projections, and has been widely used in geospatial software for many years. The pyproj library makes it possible to access the functionality of PROJ.4 from within your Python programs.

Design

The pyproj library consists of the following pieces:

pyproj consists of just two classes: Proj and Geod. Proj converts from longitude and latitude values to native map (x, y) coordinates, and vice versa. Geod performs various Great Circle distance and angle calculations. Both are built on top of the PROJ.4 library. Let's take a closer look at these two classes.

Proj

Proj is a cartographic transformation class, allowing you to convert geographic coordinates (that is, latitude and longitude values) into cartographic coordinates (x, y values, by default in meters) and vice versa.

When you create a new Proj instance, you specify the projection, datum, and other values used to describe how the projection is to be done. For example, to use the Transverse Mercator projection and the WGS84 ellipsoid, you would do the following:

projection = pyproj.Proj(proj='tmerc', ellps='WGS84')

Once you have created a Proj instance, you can use it to convert a latitude and longitude to an (x, y) coordinate using the given projection. You can also use it to do an inverse projection—that is, converting from an (x, y) coordinate back into a latitude and longitude value again.

The helpful transform() function can be used to directly convert coordinates from one projection to another. You simply provide the starting coordinates, the Proj object that describes the starting coordinates' projection, and the desired ending projection. This can be very useful when converting coordinates, either singly or en masse.

Geod

Geod is a geodetic computation class, which allows you to perform various Great Circle calculations. We looked at Great Circle calculations earlier, when considering how to accurately calculate the distance between two points on the Earth's surface. The Geod class, however, can do more than this:

  • The fwd() method takes a starting point, an azimuth (angular direction) and a distance, and returns the ending point and the back azimuth (the angle from the end point back to the start point again):

     

  • The inv() method takes two coordinates and returns the forward and back azimuth as well as the distance between them:

     

  • The npts() method calculates the coordinates of a number of points spaced equidistantly along a geodesic line running from the start to the end point:

     

When you create a new Geod object, you specify the ellipsoid to use when performing the geodetic calculations. The ellipsoid can be selected from a number of predefined ellipsoids, or you can enter the parameters for the ellipsoid (equatorial radius, polar radius, and so on) directly.

Example code

The following example starts with a location specified using UTM zone 17 coordinates. Using two Proj objects to define the UTM Zone 17 and lat/long projections, it translates this location's coordinates into latitude and longitude values:

import pyproj UTM_X = 565718.5235 UTM_Y = 3980998.9244 srcProj = pyproj.Proj(proj="utm", zone="11", ellps="clrk66", units="m") dstProj = pyproj.Proj(proj="longlat", ellps="WGS84", datum="WGS84") long,lat = pyproj.transform(srcProj, dstProj, UTM_X, UTM_Y) print "UTM zone 11 coordinate (%0.4f, %0.4f) = %0.4f, %0.4f" \ % (UTM_X, UTM_Y, lat, long)

Continuing on with this example, let's take the calculated lat/long values and, using a Geod object, calculate another point 10 kilometers northeast of that location:

angle = 315 # 315 degrees = northeast. distance = 10000 geod = pyproj.Geod(ellps="WGS84") long2,lat2,invAngle = geod.fwd(long, lat, angle, distance) print "%0.4f, %0.4f is 10km northeast of %0.4f, %0.4f" \ % (lat2, long2, lat, long)

Documentation

The documentation available on the pyproj website, and in the docs directory provided with the source code, is excellent as far as it goes. It describes how to use the various classes and methods, what they do and what parameters are required. However, the documentation is rather sparse when it comes to the parameters used when creating a new Proj object. As the documentation says:

A Proj class instance is initialized with proj map projection control parameter key/value pairs. The key/value pairs can either be passed in a dictionary, or as keyword arguments, or as a proj4 string (compatible with the proj command).

The documentation does provide a link to a website listing a number of standard map projections and their associated parameters, but understanding what these parameters mean generally requires you to delve into the PROJ documentation itself. The documentation for PROJ is dense and confusing, even more so because the main manual is written for PROJ Version 3, with addendums for later versions. Attempting to make sense of all this can be quite challenging.

Fortunately, in most cases you won't need to refer to the PROJ documentation at all. When working with geospatial data using GDAL or OGR, you can easily extract the projection as a "proj4 string" which can be passed directly to the Proj initializer. If you want to hardwire the projection, you can generally choose a projection and ellipsoid using the proj="..." and ellps="..." parameters, respectively. If you want to do more than this, though, you will need to refer to the PROJ documentation for more details.

To find out more about PROJ, and to read the original documentation, you can find everything you need at: http://trac.osgeo.org/proj

Python Geospatial Development - Second Edition Learn to build sophisticated mapping applications from scratch using Python tools for geospatial development with this book and ebook
Published: May 2013
eBook Price: $29.00
Book Price: $49.99
See more
Select your format and quantity:

Availability

Prebuild versions of pyproj are available for MS Windows, with source code distributions for other platforms. The main web page for pyproj can be found at:

http://code.google.com/p/pyproj

How you go about installing it depends on which operating system you are running.

Make sure that you install Version 4.8.0 or later of the PROJ framework, and Version 1.9.2 or later of the pyproj library. These versions are required to follow the examples in this book.

  • MS Windows

    For computers running MS Windows, installation is easy: just go to the downloads page at the website mentioned earlier and and choose the appropriate installer for your version of Python. The installer includes everything you need, including the PROJ framework.

  • Linux

    For computers running Linux, you have to download and install the PROJ framework separately, before installing pyproj. For Linux machines, you can generally obtain PROJ.4 as an RPM or source tarball which you can then compile yourself. Once this has been done, you can download the pyproj source code from the above website, and compile and install it in the usual way:

    python setup.py build python setup.py install

  • Macintosh

    If your computer runs Mac OS X, you will also have to download and install PROJ separately. You can install a compiled version of the PROJ framework either as part of a "GDAL Complete" installation, or by just installing the PROJ framework by itself. Either are available at:

    http://www.kyngchaos.com/software/frameworks

Once you have installed PROJ.4, you will have to download and build your own copy of the pyproj library. Before you can compile pyproj, you will need to have Apple's developer tools installed. Doing this is a two-step process:

  1. Download and install the latest version of XCode. XCode is available for free from the App store, or if you are running an older version of OS X you can download it from:

    https://developer.apple.com/xcode

  2. Run XCode, and choose the Preferences command. Within the Downloads tab, click on the Install button beside the Command Line Tools item:

This installs the command-line tools you will need to compile pyproj.

Once you have the developer tools installed, download the source code to pyproj from the website mentioned earlier. Then open up a Terminal window and cd into the main source code directory, then type the following commands:

python setup.py build sudo python.setup.py install

The sudo command allows pyproj to install itself inside your Python installation's site-packages directory. You'll be asked to enter your password before this is done.

Once this has finished, you can check that it worked by running the Python interpreter and typing the following command:

This will open up a GUI window allowing you

The Python prompt should reappear without any error messages being shown.

Analyzing and manipulating geospatial data

Because geospatial data works with geometrical features such as points, lines, and polygons, you often need to perform various calculations using these geometrical features. Fortunately, there are some very powerful tools for doing exactly this. For reasons we will describe shortly, the library of choice for performing this type of computational geometry in Python is Shapely .

Shapely

Shapely is a Python package for the manipulation and analysis of two-dimensional geospatial geometries. Shapely is based on the GEOS library, which implements a wide range of geospatial data manipulations in C++. GEOS is itself based on a library called the Java Topology Suite, which provides the same functionality for Java programmers. Shapely provides a Pythonic interface to GEOS which makes it easy to use these manipulations directly from your Python programs.

Design

The Shapely library is organized as follows:

All of Shapely's functionality is built on top of GEOS. Indeed, Shapely requires GEOS to be installed before it can run.

Shapely itself consists of eight major classes, representing different types of geometrical shapes:

  • The Point class represents a single point in space. Points can be two-dimensional (x, y), or three-dimensional (x, y, z).

  • The LineString class represents a sequence of points joined together to form a line. LineStrings can be simple (no crossing line segments) or complex (where two line segments within the LineString cross).

  • The LinearRing class represents a line string which finishes at the starting point. The line segments within a LinearRing cannot cross or touch.

  • The Polygon class represents a filled area, optionally with one or more "holes" inside it.

  • The MultiPoint class represents a collection of Points.

  • The MultiLineString class represents a collection of LineStrings.

  • The MultiPolygon class represents a collection of Polygons.

  • The GeometryCollection class represents a collection of any combination of Points, LineStrings, LinearRings, and Polygons.

As well as being able to represent these various types of geometries, Shapely provides a number of methods and attributes for manipulating and analyzing these geometries. For example, the LineString class provides a length attribute that equals the length of all the line segments that make up the LineString, and a crosses() method that returns true if two LineStrings cross. Other methods allow you to calculate the intersection of two polygons, dilate or erode geometries, simplify a geometry, calculate the distance between two geometries, and build a polygon that encloses all the points within a given list of geometries (called the convex_hull attribute).

Note that Shapely is a spatial manipulation library rather than a geospatial manipulation library. It has no concept of geographical coordinates. Instead, it assumes that the geospatial data has been projected onto a two-dimensional Cartesian plane before it is manipulated, and the results can then be converted back into geographic coordinates if desired.

Example code

The following program creates two Shapely geometry objects, a circle and a square, and calculates their intersection:

The intersection will be a polygon in the shape of a quarter circle , as indicated by the dark grey portion of the preceding image:

import shapely.geometry pt = shapely.geometry.Point(0, 0) circle = pt.buffer(1.0) square = shapely.geometry.Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]) intersect = circle.intersection(square) for x,y in intersect.exterior.coords: print x,y

Notice how the circle is constructed by taking a Point geometry and using the buffer() method to create a Polygon representing the outline of a circle.

Documentation

Shapely comes with excellent documentation , with detailed descriptions, extended code samples, and many illustrations that clearly show how the various classes, methods, and attributes work.

The Shapely documentation is entirely self-contained; there is no need to refer to the GEOS documentation, or to the Java Topology Suite it is based on, unless you particularly want to see how things are done in these libraries. The only exception is that you may need to refer to the GEOS documentation if you are compiling GEOS from source and are having problems getting it to work.

Availability

Shapely will run on all major operating systems, including MS Windows, Mac OS X, and Linux. Shapely's main website can be found at:

http://pypi.python.org/pypi/Shapely

The website has everything you need, including the documentation and downloads for the Shapely library, in both source code form and prebuilt binaries for MS Windows.

If you are installing Shapely on a Windows computer, the prebuilt binaries include the GEOS library built-in. Otherwise, you will be responsible for installing GEOS before you can use Shapely.

Make sure that you install Shapely Version 1.2 or later; you will need this version to work through the examples in this book.

The GEOS library's website is at:

http://trac.osgeo.org/geos

To install GEOS in a Unix-based computer, you can either download the source code from the GEOS website and compile it yourself, or you can install a suitable RPM or APT package which includes GEOS. If you are running Mac OS X, you can either try to download and build GEOS yourself, or you can install the prebuild GEOS framework, which is available from the following website:

http://www.kyngchaos.com/software/frameworks

If you've installed the "GDAL Complete" package from the above website, you'll already have GEOS installed on your Mac OS X computer.

After installing GEOS, you need to download, compile, and install the Shapely library. This can be slightly tricky on a Mac OS X computer, so you may find the following blog post useful:

http://tumblr.pauladamsmith.com/post/17663153373

Visualizing geospatial data

It's very hard, if not impossible, to understand geospatial data unless it is turned into a visual form—that is, until it is rendered as an image of some sort. Converting geospatial data into images requires a suitable toolkit. While there are several such toolkits available, we will look at one in particular: Mapnik.

Mapnik

Mapnik is a freely-available toolkit for building mapping applications. It takes geospatial data from a PostGIS database, shapefile, or any other format supported by GDAL/OGR, and turns it into clearly-rendered, good-looking images.

There are a lot of complex issues involved in rendering maps well, and Mapnik does a good job of allowing the application developer to control the rendering process. Rules control which features should appear on the map, while "symbolizers" control the visual appearance of these features.

Mapnik allows developers to create XML stylesheets that control the map-creation process. Just as with CSS stylesheets, Mapnik's stylesheets allow you complete control over the way geospatial data is rendered. Alternatively, you can create your styles by hand if you prefer.

Mapnik itself is written in C++, though bindings are included which allow access to almost all of the Mapnik functionality via Python. Because these bindings are included in the main code base rather than being added by a third party developer, support for Python is built right into Mapnik. This makes Python eminently suited to developing Mapnik-based applications.

Mapnik is heavily used by OpenStreetMap ( http://openstreetmap.org), EveryBlock ( http://everyblock.com), among others. Since the output of Mapnik is simply an image, it is easy to include Mapnik as part of a web-based application, or you can display the output directly in a window as part of a desktop-based application. Mapnik works equally well on the desktop and on the web.

Design

When using Mapnik, the main object you are dealing with is called the Map.

A Map object has the following parts:

When creating a Map object, you assign values for the following:

  • The overall width and height of the map, in pixels.

  • The spatial reference to use for the map.

  • The background color to draw behind the contents of the map.

You then define one or more Layers which hold the map's contents. Each Layer has the following:

  • A name.

  • A Datasource object defining where to get the data for this layer from. The Datasource can be a reference to a database, or it can be a shapefile or other GDAL/OGR data source.

  • A spatial reference to use for this layer. This can be different from the spatial reference used by the map as a whole, if appropriate.

  • A list of styles to apply to this layer. Each style is referred to by name, since the styles are actually defined elsewhere (often in an XML stylesheet).

Finally, you define one or more Styles , which tell Mapnik how to draw the various layers. Each Style has a name and of a list of Rules, which make up the main part of the style's definition. Each Rule has:

  • A minimum scale and maximum scale value (called the "scale denominator"). The Rule will only apply if the map's scale is within this range.

  • A filter expression. The Rule will only apply to those features which match this filter expression.

  • A list of Symbolizers . These define how the matching features will be drawn onto the map.

There are a number of different types of Symbolizers implemented by Mapnik:

  • LineSymbolizer is used to draw a "stroke" along a line, a linear ring, or around the outside of a polygon.

  • LinePatternSymbolizer uses the contents of an image file (specified by name) to draw the "stroke" along a line, a linear ring, or around the outside of a polygon.

  • PolygonSymbolizer is used to draw the interior of a polygon.

  • PolygonPatternSymbolizer uses the contents of an image file (again specified by name) to draw the interior of a polygon.

  • PointSymbolizer uses the contents of an image file (specified by name) to draw an image at a point.

  • TextSymbolizer draws a feature's text. The text to be drawn is taken from one of the feature's attributes, and there are numerous options to control how the text is to be drawn.

  • RasterSymbolizer is used to draw raster data taken from any GDAL dataset.

  • ShieldSymbolizer draws a textual label and a point together. This is similar to the use of a PointSymbolizer to draw the image and a TextSymbolizer to draw the label, except that it ensures that both the text and the image are drawn together.

  • BuildingSymbolizer uses a pseudo-3D effect to draw a polygon, to make it appear that the polygon is a three-dimensional building.

  • MarkersSymbolizer draws blue directional arrows or SVG markers following the direction of polygon and line geometries.

When you instantiate a Symbolizer and add it to a style (either directly in code, or via an XML stylesheet), you provide a number of parameters which define how the Symbolizer should work. For example, when using the PolygonSymbolizer, you can specify the fill color, the opacity, and a "gamma" value that helps draw adjacent polygons of the same color without the boundary being shown:

p = mapnik.PolygonSymbolizer(mapnik.Color(127, 127, 0)) p.fill_opacity = 0.8 p.gamma = 0.65

If the Rule that uses this Symbolizer matches one or more polygons, those polygons will be drawn using the given color, opacity, and gamma value.

Different rules can, of course, have different Symbolizers, as well as different filter values. For example, you might set up rules which draw countries in different colors depending on their population.

Example code

The following example program displays a simple world map using Mapnik:

import mapnik symbolizer = mapnik.PolygonSymbolizer( mapnik.Color("darkgreen")) rule = mapnik.Rule() rule.symbols.append(symbolizer) style = mapnik.Style() style.rules.append(rule) layer = mapnik.Layer("mapLayer") layer.datasource = mapnik.Shapefile(file="TM_WORLD_BORDERS-0.3.shp") layer.styles.append("mapStyle") map = mapnik.Map(800, 400) map.background = mapnik.Color("steelblue") map.append_style("mapStyle", style) map.layers.append(layer) map.zoom_all() mapnik.render_to_file(map, "map.png", "png")

If you are running Mapnik Version 2.0, you should replace the import mapnik statement in the first line of this program with import mapnik2 as mapnik.

Notice that this program creates a PolygonSymbolizer to display the country polygons, and then attaches the symbolizer to a Mapnik Rule object. The Rule then becomes part of a Mapnik Style object. We then create a Mapnik Layer object, reading the layer's map data from a shapefile data source. Finally, a Mapnik Map object is created, the layer is attached, and the resulting map is rendered to a PNG-format image file:

Documentation

Mapnik's has reasonable documentation for an open source project: there are good installation guides and some excellent tutorials, but the API documentation is often confusing. The Python documentation is derived from the C++ documentation, and concentrates on describing how the Python bindings are implemented rather than how an end user would work with Mapnik using Python—there's a lot of technical details that aren't relevant to the Python programmer, and many Python-specific descriptions are missing.

The best way to get started with Mapnik is to follow the installation instructions, and then to work your way through the supplied Python-specific tutorial. You can then check out the Learning Mapnik page on the Mapnik Wiki :

http://trac.mapnik.org/wiki/LearningMapnik

It is well worth spending some time reading through the Mapnik Wiki, even though not all of it is Python-specific. It is also a good idea to look at the Python API documentation, despite its limitations. The main page lists the various classes, which are available and a number of useful functions, many of which are documented. The classes themselves list the methods and properties (attributes) you can access, and even though many of these lack Python-specific documentation, you can generally guess what they do.

Availability

Mapnik runs on all major operating systems , including MS Windows, Mac OS X, and Linux. The main Mapnik website can be found at:

http://mapnik.org

Download links are provided for downloading the Mapnik source code, which can be readily compiled if you are running on a Unix machine, and you can also download prebuilt binaries for Windows and Mac OS X.

Make sure that you install Mapnik Version 2.0 or later; you will need to use this version as you work through the examples in this book.

Summary

In this article, we looked at a number of important libraries for developing geospatial applications using Python. We learned the following:

  • GDAL is a C++ library for reading (and sometimes writing) raster-based geospatial data.

  • OGR is a C++ library for reading (and sometimes writing) vector-based geospatial data.

  • GDAL and OGR include Python bindings that are easy to use, and support a large number of data formats.

  • The PROJ.4 library, and its Pythonic pyproj wrapper, allow you to convert between geographic coordinates (points on the Earth's surface) and cartographic coordinates (x,y coordinates on a two-dimensional plane) using any desired map projection and ellipsoid.

  • The pyproj Geod class allows you to perform various geodetic calculations based on points on the Earth's surface, a given distance, and a given angle (azimuth).

  • A geospatial data manipulation library called the Java Topology Suite was originally developed for Java. This was then rewritten in C++ under the name GEOS, and there is now a Python interface to GEOS called Shapely .

  • Shapely makes it easy to represent geospatial data in the form of Points, LineStrings, LinearRings, Polygons, MultiPoints, MultiLineStrings, MultiPolygons, and GeometryCollections.

  • As well as representing geospatial data, these classes allow you to perform a variety of geospatial calculations.

  • Mapnik is a tool for producing good-looking maps based on geospatial data.

  • Mapnik can use an XML stylesheet to control the elements that appear on the map, and how they are formatted. Styles can also be created by hand if you prefer.

  • Each Mapnik style has a list of Rules which are used to identify features to draw onto the map.

  • Each Mapnik rule has a list of Symbolizers that control how the selected features are drawn.

Resources for Article :


Further resources on this subject:


Python Geospatial Development - Second Edition Learn to build sophisticated mapping applications from scratch using Python tools for geospatial development with this book and ebook
Published: May 2013
eBook Price: $29.00
Book Price: $49.99
See more
Select your format and quantity:

About the Author :


Erik Westra

Erik Westra has been a professional software developer for over 25 years, and has worked almost exclusively in Python for the past decade. Erik's early interest in graphical user-interface design led to the development of one of the most advanced urgent courier dispatch systems used by messenger and courier companies worldwide. In recent years, Erik has been involved in the design and implementation of systems matching seekers and providers of goods and services across a range of geographical areas. This work has included the creation of real-time geocoders and map-based views of constantly changing data. Erik is based in New Zealand, and works for companies worldwide

Books From Packt


Matplotlib for Python Developers
Matplotlib for Python Developers

Expert Python Programming
Expert Python Programming

Python 3 Web Development Beginner's Guide
Python 3 Web Development Beginner's Guide

Python 2.6 Text Processing: Beginners Guide
Python 2.6 Text Processing: Beginners Guide

Python Geospatial Development
Python Geospatial Development

Python Text Processing with NLTK 2.0 Cookbook
Python Text Processing with NLTK 2.0 Cookbook

Python Testing: Beginner's Guide
Python Testing: Beginner's Guide

Python Testing Cookbook
Python Testing Cookbook


No votes yet

Post new comment

CAPTCHA
This question is for testing whether you are a human visitor and to prevent automated spam submissions.
9
U
k
8
5
d
Enter the code without spaces and pay attention to upper/lower case.
Code Download and Errata
Packt Anytime, Anywhere
Register Books
Print Upgrades
eBook Downloads
Video Support
Contact Us
Awards Voting Nominations Previous Winners
Judges Open Source CMS Hall Of Fame CMS Most Promising Open Source Project Open Source E-Commerce Applications Open Source JavaScript Library Open Source Graphics Software
Resources
Open Source CMS Hall Of Fame CMS Most Promising Open Source Project Open Source E-Commerce Applications Open Source JavaScript Library Open Source Graphics Software