Python: как создать карту choropleth из шейп файла Канады?

Моя цель здесь - создать choropleth map Канады в 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, но я не знаю, как).

Здесь вы можете найти шейп файл: 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)

Вышеупомянутый script создает массив 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 dict, я могу решить эту проблему. Любые идеи?

РЕДАКТИРОВАТЬ 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'.

Как обойти это?

Ответ 1

Это напрямую не отвечает на ваш вопрос, но, надеюсь, решает вашу проблему точно так же. Вы посмотрели GeoPandas? Он предоставляет простой API для работы с формами шейпов и их построения. Вы можете копировать свой код, включая прорисовку choropleth, всего несколькими строчками:

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

В этом примере предполагается, что ваш шейп файл имеет атрибут в каждой провинции, который содержит значения, которые вы хотите построить, и атрибут называется "myvalues". Если значения не сохранены в шейп файле, вы можете использовать canada.merge для объединения вашей карты values в GeoDataframe.

Одно предостережение: в это время GeoPandas не имеет простого способа построить легенду для цветов choropleth. (вопрос, указанный здесь)

Ответ 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 должна помочь, насколько я знаю. Преобразуя запись values в диапазон 0..1, затем делая цветные записи [R, G, B], должны давать вам оттенки красного цвета.

Ответ 3

Вы упомянули путаницу о cccol, являющемся списком списков. Это список кортежей RGBA (красная, зеленая, синяя, альфа-прозрачность). Они представляют 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 в крайние цвета красно-оранжевой цветовой карты.

Теперь, когда вы создаете каждую провинцию PatchCollection, вы можете установить ее facecolor в кортеж RGBA, возвращенный color_producer:

# 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)