Plotting Geographical Data using Basemap

Packt Publishing

Python developers who want to learn Matplotlib need look no further. This book covers it all with a practical approach including lots of code and images. Take this chance to learn 2D plotting through real-world examples.

Toolkits are not present in the default Matplotlib installation (in fact, they also have a different namespace, mpl_toolkits), so we have to install Basemap separately. We can download it from http://sourceforge.net/projects/matplotlib/, under the matplotlib-toolkits menu of the download section, and then install it following the instructions in the documentation link mentioned previously.

Basemap is useful for scientists such as oceanographers and meteorologists, but other users may also find it interesting. For example, we could parse the Apache log and draw a point on a map using GeoIP localization for each connection.

We use the 0.99.3 version of Basemap for our examples.

First example

Let's start playing with the library. It contains a lot of things that are very specific, so we're going to just give an introduction to the basic functions of Basemap.

Here, we initialize a Basemap object, and we can see it has several parameters depending upon the projection chosen.

Let's see what a projection is: In order to represent the curved surface of the Earth on a two-dimensional map, a map projection is needed.

This conversion cannot be done without distortion. Therefore, there are many map projections available in Basemap, each with its own advantages and disadvantages. Specifically, a projection can be:

equal-area (the area of features is preserved)

conformal (the shape of features is preserved)

No projection can be both (equal-area and conformal) at the same time.

In this example, we have used a Lambert Conformal map. This projection requires additional parameters to work with. In this case, they are lat_1, lat_2, and lon_0.

Along with the projection, we have to provide the information about the portion of the Earth surface that the map projection will describe. This is done with the help of the following arguments:

Argument

Description

llcrnrlon

Longitude of lower-left corner of the desired map domain

llcrnrlat

Latitude of lower-left corner of the desired map domain

urcrnrlon

Longitude of upper-right corner of the desired map domain

urcrnrlat

Latitude of upper-right corner of the desired map domain

The last two arguments are:

Argument

Description

resolution

Specifies what the resolution is of the features added to the map (such as coast lines, borders, and so on), here we have chosen high resolution (h), but crude, low, and intermediate are also available.

area_thresh

Specifies what the minimum size is for a feature to be plotted. In this case, only features bigger than 10,000 square kilometer

We draw a 20 degrees graticule of parallels and meridians for the map. Note how the labels argument controls the positions where the graticules are labeled. labels is an array having four elements:

[left, right, top, bottom]

These elements define the label of the parallels and the meridian when they intersect the borders of the plot. In this case, parallels are labeled when they intersect the left border and meridians are labeled at the bottom.

We now prepare a map using an orthogonal projection that displays the Earth in the way a satellite would see it. The additional arguments, lat_0 and lon_0, represent the points at the center of the projection (what the satellite looks down at).

# draw the borders of the mapm.drawmapboundary()# draw the coasts borders and fill the continentsm.drawcoastlines()m.fillcontinents()

We then draw the map's border (the edge of the map projection region) and the coastal lines, and then fill the continents.

# map city coordinates to map coordinatesx, y = m(lon, lat)

Here we convert the latitude and longitudes of the different cities into map domain coordinates—in particular, note that the resulting lists are values in meters on the map.

Calling a Basemap instance with arrays of longitudes and latitudes returns those locations in the native map projection coordinates.

# draw a red dot at cities coordinatesplt.plot(x, y, 'ro')

Now that we have the cities' locations in the map coordinates, we can plot a red dot at their positions.

# for each city,for city, xc, yc in zip(cities, x, y):# draw the city name in a yellow (shaded) box plt.text(xc+250000, yc-150000, city, bbox=dict(facecolor='yellow', alpha=0.5))

We also want to display the name of the city next to the point in the map. In order to do this, we use the text() function to write the name of the city (inside a nice yellow box, a bit translucent because of the alpha channel) next to the points position. Note the big numbers that are used to adapt the text's position. They need a little bit of hand tweaking and remember that they are in meters.

The following image is created as a result of executing the preceding code:

Plotting shapefiles with Basemap

Through the DATA.gov portal, the US government is releasing a huge quantity of high quality datasets that are free to use and analyze. Some of the datasets contain geographical information in a particular format: Shapefile.

A Shapefile, which commonly refers to a collection of files, is a popular geospatial data format for Geographical Information Systems (GIS).

Shapefiles store geometrical primitives such as points, lines, and polygons (the shapes) to represent a geographical feature in a dataset. Each item can also have attributes and information associated to it, which are used to describe what it represents.

We will use the dataset available at the URL http://www.data.gov/details/16. It represents the locations of the copper smelters in the world (it also contains several other attributes and characteristics about the smelters, but we are not going to use them here).

We use a map from a Miller cylindrical projection. We limit the latitude (while keeping the world-wide longitude) because the excluded areas don't have smelters, and so we have more space for the zones where they are present.

Reading a shapefile is as simple as calling the readshapefile() function and passing the shapefile location. The additional argument (in this case, copper) is the name of the map attribute that will be created to hold the shapefiles' vertices and features. m.copper will contain the smelters locations in map domain coordinates, while s contains only general information about the Shapefile.

Alerts & Offers

Series & Level

We understand your time is important. Uniquely amongst the major publishers, we seek to develop and publish the broadest range of learning and information products on each technology. Every Packt product delivers a specific learning pathway, broadly defined by the Series type. This structured approach enables you to select the pathway which best suits your knowledge level, learning style and task objectives.

Learning

As a new user, these step-by-step tutorial guides will give you all the practical skills necessary to become competent and efficient.

Beginner's Guide

Friendly, informal tutorials that provide a practical introduction using examples, activities, and challenges.

Essentials

Fast paced, concentrated introductions showing the quickest way to put the tool to work in the real world.

Cookbook

A collection of practical self-contained recipes that all users of the technology will find useful for building more powerful and reliable systems.

Blueprints

Guides you through the most common types of project you'll encounter, giving you end-to-end guidance on how to build your specific solution quickly and reliably.

Mastering

Take your skills to the next level with advanced tutorials that will give you confidence to master the tool's most powerful features.

Starting

Accessible to readers adopting the topic, these titles get you into the tool or technology so that you can become an effective user.

Progressing

Building on core skills you already have, these titles share solutions and expertise so you become a highly productive power user.