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列,但我想如果我能以某种方式将这个数组的值链接到valuesdict中指定的值,我可以解决问题.有任何想法吗？

编辑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'.

如何解决这个问题？

Answer 1

这并不直接回答你的问题,但希望能够解决你的问题.你看过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没有一种简单的方法来绘制等值线颜色的图例.(此处报道的问题)