In this tutorials, the following functions will be covered:
First, import things that is needed for the following steps.
In [1]:
import geopandas as gpd # for reading and manupulating shapefile
import matplotlib.pyplot as plt # for making figure
import colouringmap.mapping_polygon as mpoly # for making maps
# magic line for matlotlib figure to be shown inline in jupyter cell
%matplotlib inline
In [2]:
from palettable.colorbrewer.qualitative import Dark2_7 # to get the colormap for more custom manupulation in the last step
Second, read the file, and take a look on the attribute table of the shapefile.
In [3]:
grid_res = gpd.read_file('data/community_results.shp')
grid_res.head()
Out[3]:
Then, start playing with the shapefile. The following show how to prepare the figure before drawing the map on the figure.
The first line in the following cell is just a standard way to create a matplotlib figure, along with an "ax".
The second line prepare the "ax" for mapping things. Setting map_context to the context of the shapefile is to
make sure the shapefile is within the figure, so it will be shown within the figure.
The figure is just a matplotlib figure & ax, so you can also set the map context manually by something like this:
ax.set_xlim([minx, maxx])
ax.set_ylim([miny, maxy])
In [4]:
fig,ax = plt.subplots(figsize=(7,7))
ax = mpoly.prepare_map(ax, map_context=grid_res)
The following show how to create a map with just the shape of the polygons with a same colour.
In [5]:
fig,ax = plt.subplots(figsize=(7,7))
ax = mpoly.prepare_map(ax, map_context=grid_res)
ax = mpoly.map_shape(grid_res, ax, lw=.1, alpha=.7)
So, lets try to map the 'com' column to colours.
use mpoly.map_category function, throw in the geopandas gdf, the column name, and the ax.
In [6]:
fig,ax = plt.subplots(figsize=(7,7))
ax = mpoly.prepare_map(ax, map_context=grid_res)
ax = mpoly.map_category(grid_res, 'com', ax)
The resulting map is not so good.
The map return that there are too many categories.
Normally, there are two ways to coupe with this problem:
In this case, the number of category is way too high, this can be observe by the legend. There are no colormap that can support so much number of categories.
So, in the following, I try to reduce the number of category to 6 major cats, plus one category as "other".
The first part in the following cell is to find the "major" categories, which is determined by the appearance frequency of the "com".
The second part is to create a new column in the attribute table, to store the major cats data.
In [7]:
## find major categories
coms = grid_res.com.tolist() # get the cats column from the attribute table
comset = list(set(coms)) # get the unique cat id
comcount = [ coms.count(c) for c in comset ] # count the appearance of the cats
comset2 = [ (c,n) for n,c in sorted(zip(comcount,comset), reverse=True) ] # sort the cats by the frequency
major = [c[0] for c in comset2[:6]] # 6 major cats
print major
## create new column of categories with major/other
collist = []
for c in coms:
if c in major:
collist.append(c) # the cat id if it is in the major cats
else:
collist.append(-1) # or -1 for other cats
print len(collist)==len(grid_res) # just a check
grid_res['com2'] = collist # put the new column into the attribute table
Then, map the major categories, with some modification (lw, ec, alpha) for better looking on the edge of the polygon.
In [8]:
fig,ax = plt.subplots(figsize=(7,7))
ax = mpoly.prepare_map(ax, map_context=grid_res)
ax = mpoly.map_category(grid_res, 'com2', ax, lw=.1, ec='k', alpha=.7)
Actually, this can also be done by setting the cat_order paramter to the desire category list like this. (Experimental)
In [9]:
fig,ax = plt.subplots(figsize=(7,7))
ax = mpoly.prepare_map(ax, map_context=grid_res)
ax = mpoly.map_category(grid_res, 'com2', ax, lw=.1, ec='k', alpha=.7, cat_order=[0,1,3,11,14,18])
The result is not perfect yet. The colour of the "other" cat cannot be controlled. And, its colour looks just as it is as important as the other cats.
So, lets manually assign colours to the categories, and set the "other" cat to something less attractive.
In [10]:
colors = Dark2_7.mpl_colors # get this colormap from palletable
collist2 = []
for c in grid_res.com.tolist():
if c in major:
collist2.append(colors[major.index(c)]) # if the cat is a major, then use the colormap
else:
collist2.append('lightgrey') # silver colour for those "other" cat
len(collist2)==len(grid_res)
grid_res['colour'] = collist2 # assign the new list with color to the attribute table
In [11]:
grid_res.head() # lets take a look
Out[11]:
This time, map the polygon directly with the colour column. Note that the function is "map_colour()".
In [12]:
fig,ax = plt.subplots(figsize=(7,7))
ax = mpoly.prepare_map(ax, map_context=grid_res)
ax = mpoly.map_colour(grid_res, 'colour', ax, lw=.1, ec='k', alpha=.7)
The map looks better now.
But, it would be better if we also map the administrative boundaries to the map, to show which colour belong to which area. So lets add the boundary shapefile.
The following read the administrative boundary of the Tokyo 23 special wards.
In [13]:
borders = gpd.read_file('data/tokyo_special_ward.shp')
borders.head()
Out[13]:
Before we map the boundaries to the map, lets check if the two shapefile have the same projection.
In [14]:
print borders.crs==grid_res.crs # check if the two shapefile have the same projection
print borders.crs # check the two projections
print grid_res.crs
Turns out they are not same, so lets do some projection.
In [15]:
borders = borders.to_crs(grid_res.crs) # convert the borders projection to the same as the grid_res
print borders.crs==grid_res.crs # now check again if the two shapefile have the same projection
So they are now in the same projection.
Now, lets add the administrative boundaries to the map.
Because the boundaries are set to black colour (ec='k'), so the grid edge colour is changed to white (ec='w')
In [16]:
fig,ax = plt.subplots(figsize=(7,7))
ax = mpoly.prepare_map(ax, map_context=grid_res)
ax = mpoly.map_colour(grid_res, 'colour', ax, lw=.1, ec='w', alpha=.7)
ax = mpoly.add_border(borders, ax, ec='k', alpha=.4)
The map is better now.
Sometimes, we may need to add the labels of the area name to the map. So, lets try to find the name of each administrative.
In [17]:
print borders['NAME_2']
The column that record the names seems to be in the "NAME_2" column.
So, lets try to add it to the previous map.
In [23]:
fig,ax = plt.subplots(figsize=(7,7))
ax = mpoly.prepare_map(ax, map_context=grid_res)
ax = mpoly.map_colour(grid_res, 'colour', ax, lw=.1, ec='w', alpha=.7)
ax = mpoly.add_border(borders, ax, ec='k', alpha=.4)
ax = mpoly.add_label(borders, ax, 'NAME_2', font_colour='k', font_size=10)
To be honest, this map can be helpful for observation, but it is not very nice looking.
So, maybe use the previous map for publication, and this for discussion. XD
That is all in this tutorial.
In [ ]: