Python：如何从加拿大的shapefile创建一个区域分布地图？

Question

我的目标是用Python创建一个加拿大的分级统计地图。假设我有一个包含每个加拿大省/地区值的字典:

myvalues={'Alberta': 1.0,
 'British Columbia': 2.0,
 'Manitoba': 3.0,
 'New Brunswick': 4.0,
 'Newfoundland and Labrador': 5.0,
 'Northwest Territories': 6.0,
 'Nova Scotia': 7.0,
 'Nunavut': 8.0,
 'Ontario': 9.0,
 'Prince Edward Island': 10.0,
 'Quebec': 11.0,
 'Saskatchewan': 12.0,
 'Yukon': 13.0}

现在我想根据myvalues中相应的值，使用连续的彩图（例如红色阴影）对每个省份进行着色。 该怎么做？

到目前为止，我只能够在matplotlib中绘制加拿大省/地区，但它们的形状以唯一的颜色出现，我不知道如何根据myvalues中的数字更改它们的颜色（也许我需要玩弄patches，但我不知道如何）。

这是您可以找到shapefile的位置：http://www.filedropper.com/canadm1_1

这是我的迄今为止的代码：

import shapefile
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
#   -- input --
sf = shapefile.Reader("myfolder\CAN_adm1.shp")
recs    = sf.records()
shapes  = sf.shapes()
Nshp    = len(shapes)
cns     = []
for nshp in xrange(Nshp):
    cns.append(recs[nshp][1])
cns = array(cns)
cm    = get_cmap('Dark2')
cccol = cm(1.*arange(Nshp)/Nshp)
#   -- plot --
fig     = plt.figure()
ax      = fig.add_subplot(111)
for nshp in xrange(Nshp):
    ptchs   = []
    pts     = array(shapes[nshp].points)
    prt     = shapes[nshp].parts
    par     = list(prt) + [pts.shape[0]]
    for pij in xrange(len(prt)):
     ptchs.append(Polygon(pts[par[pij]:par[pij+1]]))
    ax.add_collection(PatchCollection(ptchs,facecolor=None,edgecolor='k', linewidths=.5))
ax.set_xlim(-160,-40)
ax.set_ylim(40,90)

目前我得到的图片是:

编辑

我认为解决方案必须在以下行中:

cm    = get_cmap('OrRd')
cccol = cm(1.*arange(Nshp)/Nshp)

上述脚本创建了一个cccol数组，实际上它的形状如下：

array([[ 1.        ,  0.96862745,  0.9254902 ,  1.        ],
       [ 0.99766244,  0.93356402,  0.84133796,  1.        ],
       [ 0.99520185,  0.89227221,  0.74749713,  1.        ],
       [ 0.99274125,  0.84306037,  0.64415227,  1.        ],
       [ 0.99215686,  0.78754327,  0.5740254 ,  1.        ],
       [ 0.99186467,  0.71989237,  0.50508269,  1.        ],
       [ 0.98940408,  0.60670514,  0.39927722,  1.        ],
       [ 0.97304114,  0.50618995,  0.32915034,  1.        ],
       [ 0.94105344,  0.40776625,  0.28732027,  1.        ],
       [ 0.88521339,  0.28115341,  0.19344868,  1.        ],
       [ 0.8220992 ,  0.16018455,  0.10345252,  1.        ],
       [ 0.73351789,  0.04207613,  0.02717416,  1.        ],
       [ 0.61959248,  0.        ,  0.        ,  1.        ]])

我不知道为什么它有4列，但我想如果我能将此数组的值与values字典中指定的那些值链接起来，我就可以解决这个问题。 有什么想法吗？

编辑2

我已经找到了“诀窍”的方法，就在cccol = cm()中。 为了将其与省份相关联，我尝试分配 cccol = cm(myvalues.values(i) for i in myvalues.keys())

这样（至少在我的脑海中）每个颜色都是基于相关键分配的，没有任何错位。 问题是我遇到了错误：

TypeError：Cannot cast array data from dtype('O') to dtype('int32') according to the rule 'safe'。

如何解决这个问题？

- FaCoffee

---

您提供的shapefile无效，因为它只包含*.shp文件，而没有附带的文件（例如CAN_adm1.shx、CAN_adm1.dbf等）。您能否提供一个包含所有文件的ZIP文件链接？ - Jeff G

@JeffG 这是压缩文件夹的链接：http://www.filedropper.com/canadm1_1。非常抱歉，我会在问题本身中修改链接。谢谢。 - FaCoffee

Answer 1

这并没有直接回答你的问题，但希望能解决你的问题。你看过GeoPandas吗？它提供了一个简单的API来处理和绘制shapefiles。你可以用几行代码复制你的代码，包括绘制choropleth。

import geopandas as gpd
canada = gpd.read_file('CAN_adm1.shp')
canada.plot('myvalues', cmap='OrRd')

该示例假设您的Shapefile上有一个包含要绘制的值的每个省的属性，该属性称为"myvalues"。如果值未存储在Shapefile中，则可以使用canada.merge将您的values地图合并到GeoDataframe上。

请注意：目前，GeoPandas没有一种简单的方法来绘制分级颜色图例。（已报告的问题见此处）

- Jeff G

---

geopandas 在这里是一个不错的解决方案。这是我通过将 myvalues 与 canada 合并而得到的图像，正如 Jeff 所建议的：http://i.stack.imgur.com/KiAry.png - lanery

Answer 2

请求: 请将您的values字典改名为其他名称。该名称使得撰写此答案更加困难 :)

我还没有测试过，但可以尝试以下操作:

color_numbers = values.values()
    # assumes the provinces are listed in the same order in values as 
    # they are in the shape file
for nshp in xrange(Nshp):
    ptchs   = []
    # ... code omitted ...
    the_facecolor = [(color_numbers[nshp]-1)/(Nshp-1), 0, 0];   #1..13 -> 0..1, then add G=B=0.
        # change the computation if the values in the values dictionary are no longer 1..13
    ax.add_collection(PatchCollection(ptchs, facecolor=the_facecolor, edgecolor='k', linewidths=.5))

你正在获得的输出中有所有的蓝色补丁，或者说[0,0,1]。由于该行不在cccol中，我不认为cccol是问题所在。此外，你添加的代码在创建后从未实际引用cccol！（请添加你开始使用的代码示例的链接！ :) ）无论如何，设置facecolor应该会有所帮助。将值条目转换为0..1范围，然后制作[R,G,B]颜色条目，应该会给你红色的阴影。

- cxw

---

刚刚将 values 改为了 myvalues。希望没问题 :) - FaCoffee

有一个小缺陷，我认为：(color_numbers [nshp] - 1)将始终返回负的RGB。 我是正确的吗？ 我尝试删除-1，结果得到完全黑色的地图。 怎么回事？ - FaCoffee

你当前示例中的 color_numbers 值范围从 1.0 到 13.0 - 它们是 myvalues 冒号后面的部分。 ...-1 将这些数字缩放到 0-12，/(Nshp-1) 将其更改为 0-1。如果您使用的具体数字发生变化，则需要更改该公式。 - cxw

Answer 3

您提到了关于cccol是列表嵌套列表的疑惑。实际上，它是由RGBA元组（红、绿、蓝、alpha透明度）组成的列表。这些表示从橙色到红色的13种“等距”颜色。

在您的情况下，您不需要等距颜色，而是需要与myvalues相对应的颜色。请按照以下步骤进行操作：

cmap = matplotlib.cm.get_cmap('OrRd')
norm = matplotlib.colors.Normalize(min(myvalues.values()), max(myvalues.values()))
color_producer = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap)

现在，color_producer有一个to_rgba方法，它从myvalues中获取值并将其转换为正确的颜色。 Normalize将myvalues的最小值和最大值范围设置为Red-Orange colormap的极端颜色。

现在，当您创建每个省份的PatchCollection时，可以将其facecolor设置为color_producer返回的RGBA元组。

# Change the province name passed as you iterate through provinces.
rgba = color_producer.to_rgba(myvalues['Manitoba'])
PatchCollection(ptchs, facecolor=rgba, edgecolor='k', linewidths=.5)

- Jeff G