In this post I’ll discuss how to plot data points on a shapefile. In a previous post I discussed how to install basemap using pip, the package manager for Python. Since basemap is an extension of matplotlib, we have a lot of familiar plotting functions and options at our disposal. Of particular importance is the ability to use projection data in plotting the shapefile, and plotting the data points.

import numpy as np
import scipy.stats
from mpl_toolkits.basemap import Basemap

After we import what we need, we can create a map object, m, with some projection. Here, we’ve used the Miller cylindrical projection, which is better for stuff nearer the equator than the poles.

Next, we call the readshapefile() method of the map object, m, to incorporate the data from the shapefile. Then, after generating some uniformly distributed data points, we call the scatter() method of the map object to plot the data points.

# the map, a Miller Cylindrical projection
m = Basemap(projection='mill', llcrnrlon=-104. ,llcrnrlat=28, urcrnrlon=-96. ,urcrnrlat=34.)

# incorporate the shapefile
s = m.readshapefile('./Texas_Counties', 'texascounties')

# generate some points
x = scipy.stats.uniform(-102,2).rvs(3)
y = scipy.stats.uniform(30,2).rvs(3)

# superimpose the points on the shapefile
m.scatter( x, y, marker='o', facecolor='none', edgecolor='b', alpha=0.5, latlon=True )

The latlon=True argument of the scatter() function is particularly important. This informs scatter() that the data points should be interpreted as latitude and longitude. Beyond that, the rest of the scatter() arguments are familiar from the matplotlib library.