# Plotting Geographical Data using Basemap

#### Matplotlib for Python Developers

\$26.99

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.

Basemap is a Matplotlib toolkit, a collection of application-specific functions that extends Matplotlib functionalities, and its complete documentation is available at http://matplotlib.sourceforge.net/basemap/doc/html/index.html.

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.

# pyplot module import
import matplotlib.pyplot as plt
# basemap import
from mpl_toolkits.basemap import Basemap
# Numpy import
import numpy as np

These are the usual imports along with the basemap module.

# Lambert Conformal map of USA lower 48 states
m = Basemap(llcrnrlon=-119, llcrnrlat=22, urcrnrlon=-64,
urcrnrlat=49, projection='lcc', lat_1=33, lat_2=45,
lon_0=-95, resolution='h', area_thresh=10000)

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

# draw the coastlines of continental area
m.drawcoastlines()
# draw country boundaries
m.drawcountries(linewidth=2)
# draw states boundaries (America only)
m.drawstates()

We start adding features to the map. In this case, we have just added:

• The coast lines
• The country borders (with a bigger line style)
• The state borders inside the country (they are only available for America)
# fill the background (the oceans)
m.drawmapboundary(fill_color='aqua')
# fill the continental area
# we color the lakes like the oceans
m.fillcontinents(color='coral',lake_color='aqua')

We give some colors to our map. We color the ocean with aqua color and the interior of the continents are coral (but lakes have the same color of the ocean).

# draw parallels and meridians
m.drawparallels(np.arange(25,65,20),labels=[1,0,0,0])
m.drawmeridians(np.arange(-120,-40,20),labels=[0,0,0,1])

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.

After adding a title to it, the result is as shown:

## Using satellite background

Basemap can also use a terrain image as a map background.

m = Basemap(llcrnrlon=-119, llcrnrlat=22, urcrnrlon=-64,
urcrnrlat=49, projection='lcc', lat_1=33, lat_2=45,
lon_0=-95, resolution='h', area_thresh=10000)

We are using the same map as before.

# display blue marble image (from NASA) as map background
m.bluemarble()

We add the satellite images taken from the NASA images library.

# draw the coastlines of continental area
m.drawcoastlines()
# draw country boundaries
m.drawcountries(linewidth=2)
# draw states boundaries (America only)
m.drawstates()

Then, we draw the Basemap features over it, as done before, and this results in a very pretty image, as shown in the following screenshot:

## Plot data over a map

We got to know Matplotlib as a tool that can plot datasets easily. It would be really nice if we can mix this with Basemap projections. Well, of course it's possible, and this is how it is done.

We need some geographical data to plot over a map, so we take a group of cities, and we'll plot them on a map. The cities along with their coordinates (taken from Wikipedia) are:

City

Latitude

Longitude

London

51° 30′ 28″ N

0° 7′ 41″ W

New York

40° 43′ 0″ N

74° 0′ 0″ W

40° 24′ 0″ N

3° 41′ 0″ W

Cairo

30° 3′ 28.8″ N

31° 13′ 44.4″ E

Moscow

55° 45′ 6″ N

37° 37′ 4″ E

Delhi

28° 36′ 36″ N

77° 13′ 48″ E

Dakar

14° 41′ 34″ N

17° 26′ 48″ W

From the previous table, we can prepare these three lists:

# Cities names and coordinates
cities = ['London', 'New York', 'Madrid', 'Cairo', 'Moscow',
'Delhi', 'Dakar']
lat = [51.507778, 40.716667, 40.4, 30.058, 55.751667,
28.61, 14.692778]
lon = [-0.128056, -74, -3.683333, 31.229, 37.617778,
77.23, -17.446667]

where we have recorded the names and coordinates of different cities.

# orthogonal projection of the Earth
m = Basemap(projection='ortho', lat_0=45, lon_0=10)

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 map
m.drawmapboundary()
# draw the coasts borders and fill the continents
m.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 coordinates
x, 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 coordinates
plt.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).

# the map, a Miller Cylindrical projection
m = Basemap(projection='mill',
llcrnrlon=-180. ,llcrnrlat=-60,
urcrnrlon=180. ,urcrnrlat=80.)

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.

# read the shapefile archive
s = m.readshapefile('<location of shapefile>/copper', 'copper')

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.

# prepare map coordinate lists for copper smelters locations
x, y = zip(*m.copper)

We prepare a list of coordinates (in the map domain) for the copper smelters locations; zip() receives the m.copper array unpacked (each sublist is passed as a separate parameter to zip()).

# draw coast lines and fill the continents
m.drawcoastlines()
m.fillcontinents()

We draw the coast lines and fill the continents

# draw a blue dot at smelters location
plt.plot(x, y, 'b.')

We can then draw a blue dot at the smelters' locations.

When we run the example, we can see a map with dots (in blue) that represent the places where the smelters are located:

# Summary

In this article, we have seen how to plot geographical data using Basemap.

If you have read this article you may be interested to view :

### Books to Consider

Oracle Data Integrator Essentials [Video]
\$ 72.25
Microsoft Data Protection Manager 2010
\$ 35.99