Last week, I needed to draw some maps and display data on these maps also. I have spent some hours to discover how to complete this task. Thus, in this post, I will summarize step by step to accomplish the initial goal.

## Basic drawing

2. Secondly, we need to choose the library to draw. BaseMap is a common lib for drawing shapefile but it is overkilled in my situation i.e. I do not need to use its fancy features. Thus, pyshp is my selection. It is lightweight but still powerful enough to satisfy my requirement. Moreover, installation and usage are as easy as ABC :-) . Last but not least, I have used Python 2.7 and matplotlib 1.4.2 to draw but the conversion to Python 3.x is not hard.
3. Let us draw :-D
# import libraries
import shapefile as shp
import matplotlib.pyplot as plt

plt.figure()
for shape in sf.shapeRecords():

# end index of each components of map
l = shape.shape.parts

len_l = len(l)  # how many parts of countries i.e. land and islands
x = [i[0] for i in shape.shape.points[:]] # list of latitude
y = [i[1] for i in shape.shape.points[:]] # list of longitude
l.append(len(x)) # ensure the closure of the last component
for k in xrange(len_l):
# draw each component of map.
# l[k] to l[k + 1] is the range of points that make this component
plt.plot(x[l[k]:l[k + 1]],y[l[k]:l[k + 1]], 'k-')

# display
plt.show()


The final result is not bad.

## Display histogram

Using the same setup as above, we will display the histogram in 2D map to illustrate the density. There are some tricks to complete the task so let us view the code first.

import shapefile as shp
import matplotlib.pyplot as plt
import matplotlib as mpl
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np

def drange(start, stop, step):
result = list()
r = start
i = 0.0
while r < stop:
result.append(r)
i += 1.0
r = round(start + i * step, 3)
result.append(r)
return result

f = open(fname, 'r')
f.close()
x = list()
y = list()
for line in lines:
comp = line.strip('\n').split(',')
x.append(float(comp[1]))
y.append(float(comp[0]))
return x, y

plt.figure()
for shape in sf.shapeRecords():
l = shape.shape.parts
len_l = len(l)  # how many parts of countries i.e. land and
x = [i[0] for i in shape.shape.points[:]]
y = [i[1] for i in shape.shape.points[:]]
l.append(len(x))
for k in xrange(len_l):
plt.plot(x[l[k]:l[k + 1]],y[l[k]:l[k + 1]], 'k-')

lower_x = 103.6
lower_y = 1.2
upper_x = 104.1
upper_y = 1.5
# plot the frame of city
x = [lower_x, lower_x, upper_x, upper_x, lower_x]
y = [lower_y, upper_y, upper_y, lower_y, lower_y]
plt.plot(x, y)

fname = '1.35_103.9_1.4_103.95'
comp = fname.split('_')
loc = [float(i) for i in comp]

# plot the line
area_x = [loc[1], loc[1], loc[3], loc[3], loc[1]]
area_y = [loc[0], loc[2], loc[2], loc[0], loc[0]]
# plot a box in map
plt.plot(area_x, area_y, 'b--')

# histogram
xedges = drange(lower_x, upper_x, 0.01) # 0.01 is the size of bin
yedges = drange(lower_y, upper_y, 0.01)
H, xedges, yedges = np.histogram2d(x, y, bins=(xedges, yedges))

im = plt.imshow(H.T, cmap='OrRd', interpolation='bilinear', origin='low', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])

# add label to axis and do not scale unit of x, y axises
plt.ylabel('Latitude', fontsize=20)
plt.xlabel('Longitude', fontsize=20)
ax = plt.gca()
ax.ticklabel_format(useOffset=False)
divider = make_axes_locatable(ax)                          # set size of color bar
cax = divider.append_axes("right", size="5%", pad=0.05)    # set size of color bar
plt.colorbar(im, cax=cax)                                  # set size of color bar

# finally, show the plot
plt.show()


The file user_density.txt could be downloaded via this link. The data inside this file is generated randomly. Each line of file contains 2 information latitude and longitude of each point and both are separated via comma. Function readLoc is used to load its data.

There are 2 important notes for the above code

• I do not use the normal order in numpy manual which is (y, x)
H, xedges, yedges = np.histogram2d(x, y, bins=(xedges, yedges))

• We use the transpose matrix of H.T when using imshow. Option bilinear is used to make it look like heatmap.
im = plt.imshow(H.T, cmap='OrRd', interpolation='bilinear', origin='low', extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]])


The final result is quite beautiful