Plotting Maps: Visualizing Haiti Earthquake Crisis Data (2)

Haiti_shape.py

import pandas as pd
import numpy as np
from pandas import DataFrame

data = pd.read_csv('ch08/Haiti.csv')

data = data[(data.LATITUDE > 18) & (data.LATITUDE < 20) &
		(data.LONGITUDE > -75) & (data.LONGITUDE < -70)
		& data.CATEGORY.notnull()]

def to_cat_list(catstr):
	stripped = (x.strip() for x in catstr.split(','))
	return [x for x in stripped if x]

def get_all_categories(cat_series):
	cat_sets = (set(to_cat_list(x)) for x in cat_series) 
	return sorted(set.union(*cat_sets))

def get_english(cat):
	code, names = cat.split('.') 
	if '|' in names:
		names = names.split(' | ')[1] 
	return code, names.strip()

all_cats = get_all_categories(data.CATEGORY)
english_mapping = dict(get_english(x) for x in all_cats)

def get_code(seq):
	return [x.split('.')[0] for x in seq if x]

all_codes = get_code(all_cats)
code_index = pd.Index(np.unique(all_codes))
dummy_frame = DataFrame(np.zeros((len(data), len(code_index))),
						index=data.index, columns=code_index)

for row, cat in zip(data.index, data.CATEGORY): 
	codes = get_code(to_cat_list(cat)) 
	dummy_frame.ix[row, codes] = 1

data = data.join(dummy_frame.add_prefix('category_'))

from mpl_toolkits.basemap import Basemap 
import matplotlib.pyplot as plt

def basic_haiti_map(ax=None, lllat=17.25, urlat=20.25, lllon=-75, urlon=-71):
	# create polar stereographic Basemap instance. 
	m = Basemap(ax=ax, projection='stere', 
				lon_0=(urlon + lllon) / 2, 
				lat_0=(urlat + lllat) / 2,
				llcrnrlat=lllat, urcrnrlat=urlat, 
				llcrnrlon=lllon, urcrnrlon=urlon, 
				resolution='f')
	# draw coastlines, state and country boundaries, edge of map. 
	#m.drawcoastlines()
	m.drawstates()
	m.drawcountries()
	return m

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 10)) 
fig.subplots_adjust(hspace=0.05, wspace=0.05)

lllat=18.4; urlat=18.7; lllon=-72.5; urlon=-72.1

m = basic_haiti_map(ax, lllat=lllat, urlat=urlat,
					lllon=lllon, urlon=urlon) 

cat_data = data[data['category_2a'] == 1]

# compute map proj coordinates.
x, y = m(cat_data.LONGITUDE.values, cat_data.LATITUDE.values)

m.plot(x, y, 'k.', alpha=0.5)
ax.set_title('Food shortages reported in Port-au-Prince')

shapefile_path = 'ch08/PortAuPrince_Roads/PortAuPrince_Roads'
m.readshapefile(shapefile_path, 'roads', drawbounds=True, linewidth=0.5)

plt.show()

This is an example code for "Plotting Maps: Visualizing Haiti Earthquake Crisis Data" on a book, Python for Data Analysis. Page 246-247

Haiti.csv and Shapefile can be downloaded from https://github.com/pydata/pydata-book/tree/master/ch08

The resulting figure can be seen at http://i.imgur.com/492oM4w.png