Python：如何根据加拿大的 shapefile 创建等值线图？

Question

我的目标是在 Python 中创建一个 choropleth map 的加拿大。假设我有一本字典，其中包含引用每个加拿大人的值 province/territory:

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 中绘制加拿大 provinces/territory，但它们的形状以独特的颜色出现，我不知道如何根据 [=] 中的数字更改它15=]（也许我需要玩 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 列，但我想如果我能以某种方式 link 这个数组的值到 values 字典中指定的值，我可以解决问题。有什么想法吗？

编辑 2

我发现 "trick" 在 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'。

如何解决这个问题？

Answer 1

请求：请将您的 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！（请将 link 添加到您开始的代码示例中！:)）

无论如何，据我所知，设置 facecolor 应该会有帮助。将 values 条目转换为范围 0..1，然后制作 [R,G,B] 颜色条目，应该会给你红色阴影。

Answer 2

您提到了关于 cccol 是列表列表的困惑。它是 RGBA 元组列表（红色、绿色、蓝色、alpha 透明度）。这些代表从橙色到红色的 13 "equally spaced" 种颜色。

在您的情况下，您不需要等间距的颜色，而是与 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 的最小和最大范围设置为红橙色颜色图的极端颜色。

现在当你创建每个省份的 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)

Answer 3

这并没有直接回答您的问题，但希望能同样解决您的问题。你看过GeoPandas了吗？它提供了一个简单的 API 用于处理和绘制 shapefile。您可以复制您的代码，包括绘制等值线，只需几行：

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

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

一个警告：目前 GeoPandas 没有一种简单的方法来绘制等值线颜色的图例。 (issue reported here)

Python：如何根据加拿大的 shapefile 创建等值线图？

Python: how to create a choropleth map out of a shapefile of Canada?

python

patch

matplotlib

shapefile

choropleth