# Plotting Geographical Data using Basemap

Your details (so we can tell your friend who this is from) *
This question is for testing whether you are a human visitor and to prevent automated spam submissions.
Q
L
f
r
Q
h
by Sandro Tosi | October 2009 | Open Source

This article by Sandro Tosi is dedicated to Basemap, a Matplotlib toolkit to draw geographical data.

We can use Matplotlib to draw on geographical map projections using the Basemap external toolkit. Basemap provides an efficient way to draw Matplotlib plots over real world maps.

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 importimport matplotlib.pyplot as plt# basemap importfrom mpl_toolkits.basemap import Basemap# Numpy importimport numpy as np`

These are the usual imports along with the basemap module.

`# Lambert Conformal map of USA lower 48 statesm = 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

## Matplotlib for Python Developers

 Build remarkable publication-quality plots the easy way
Published: November 2009
eBook Price: £16.99
Book Price: £27.99
See more
`# draw the coastlines of continental aream.drawcoastlines()# draw country boundariesm.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 oceansm.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 meridiansm.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 backgroundm.bluemarble()`

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

`# draw the coastlines of continental aream.drawcoastlines()# draw country boundariesm.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:

## Matplotlib for Python Developers

 Build remarkable publication-quality plots the easy way
Published: November 2009
eBook Price: £16.99
Book Price: £27.99
See more

## 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 coordinatescities = ['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 Earthm = 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 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).

`# the map, a Miller Cylindrical projectionm = 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 archives = 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 locationsx, 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 continentsm.drawcoastlines()m.fillcontinents()`

We draw the coast lines and fill the continents

`# draw a blue dot at smelters locationplt.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.

## Sandro Tosi

Sandro Tosi is a Debian Developer, Open Source evangelist, and Python enthusiast. After completing a B.Sc. in Computer Science from the University of Firenze, he worked as a consultant for an energy multinational as System Analyst and EAI Architect, and now works as System Engineer for one of the biggest and most innovative Italian Internet companies.

## Books From Packt

well this was very dissapointing by
As this same sort of example floats around the web, what is the point? clearly i donot know who was first but this kind of info is duplicated what you should do is show how to project an square image of e.g. surface temperature
how to draw links between two cities by
Hi Nice information on basemap :) I would like to know how to draw links or add edges between two cities e.g. New York and New Delhi.