Как использовать набор обрезанных путей для многоугольника базовой карты

Я хочу использовать imshow (например) для отображения некоторых данных внутри границ страны (для примера я выбрал США). Простой пример ниже иллюстрирует то, что я хочу:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import RegularPolygon

data = np.arange(100).reshape(10, 10)
fig = plt.figure()
ax = fig.add_subplot(111)
im = ax.imshow(data)
poly = RegularPolygon([ 0.5,  0.5], 6, 0.4, fc='none', 
                      ec='k', transform=ax.transAxes)
im.set_clip_path(poly)
ax.add_patch(poly)
ax.axis('off')
plt.show()

Результат:

введите описание изображения здесь

Теперь я хочу сделать это, но вместо простого многоугольника я хочу использовать сложную форму США. Я создал несколько примеров данных, содержащихся в массиве «Z», как показано в приведенном ниже коде. Именно эти данные я хочу отобразить с помощью цветовой карты, но только в пределах материковой части США.

Пока что пробовал следующее. Я получаю файл формы здесь, содержащийся в файле "nationp010g.shp.tar. gz "и я использую модуль Basemap в Python для построения США. Обратите внимание, что это единственный найденный мной метод, который дает мне возможность получить многоугольник в нужной мне области. Если есть альтернативные методы, меня тоже заинтересуют. Затем я создаю многоугольник с именем mainpoly, который является почти тем многоугольником, который я хочу покрасить в синий цвет:

введите описание изображения здесь

Обратите внимание, как было окрашено только одно тело, все остальные непересекающиеся многоугольники остались белыми:

введите описание изображения здесь

Так что область, окрашенная в синий цвет, - это почти то, что я хочу, обратите внимание, что рядом с Канадой есть нежелательные границы, потому что граница на самом деле проходит через некоторые озера, но это небольшая проблема. Настоящая проблема в том, почему мои данные imshow не отображаются внутри США? Сравнивая мой первый и второй примеры кода, я не понимаю, почему я не получаю обрезанное imshow во втором примере, как в первом. Любая помощь будет оценена в понимании того, что мне не хватает.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap as Basemap
from matplotlib.patches import Polygon

# Lambert Conformal map of lower 48 states.
m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,
            projection='lcc',lat_1=33,lat_2=45,lon_0=-95)


shp_info = m.readshapefile('nationp010g/nationp010g', 'borders', drawbounds=True) # draw     country boundaries.

for nshape,seg in enumerate(m.borders):
    if nshape == 1873: #This nshape denotes the large continental body of the USA, which we want
        mainseg = seg
        mainpoly =  Polygon(mainseg,facecolor='blue',edgecolor='k')



nx, ny = 10, 10
lons, lats = m.makegrid(nx, ny) # get lat/lons of ny by nx evenly space grid.
x, y = m(lons, lats) # compute map proj coordinates.

Z = np.zeros((nx,ny))
Z[:] = np.NAN

for i in np.arange(len(x)):
    for j in np.arange(len(y)):
        Z[i,j] = x[0,i] 

ax = plt.gca()
im = ax.imshow(Z, cmap = plt.get_cmap('coolwarm') )
im.set_clip_path(mainpoly)
ax.add_patch(mainpoly)
plt.show()

Обновлять

Я понимаю, что линия

ax.add_patch(mainpoly)

не добавляет даже многоугольника к графику. Я не правильно использую? Насколько мне известно, mainpoly правильно рассчитывалась с использованием метода Polygon (). Я проверил, что ввод координат разумный:

plt.plot(mainseg[:,0], mainseg[:,1] ,'.') 

который дает

введите описание изображения здесь


person Dipole    schedule 06.09.2014    source источник
comment
Мне также любопытно, почему я опускаю голоса. Пожалуйста, скажите мне, чтобы я мог стать лучше!   -  person Dipole    schedule 06.09.2014
comment
Что ты пробовал? Что ты спрашиваешь? Действительно ли вопрос: «Как мне упростить путь?» ? Если да, то почему важно, чтобы путь лежал в США? Вы смотрели какие-нибудь библиотеки геометрии? Предположительно, вы хотите закрыть Чесапикский залив, возможно, пролив Лонг-Айленд, но не Мексиканский залив. Вопрос супер открытый, не требует больших исследовательских усилий и читается как «пожалуйста, сделайте мою работу за меня, дайте мне код !! 1!» циничным / раздражительным / пресыщенным членам SO.   -  person tacaswell    schedule 06.09.2014
comment
Хорошо, я обновлю дополнительную информацию, спасибо.   -  person Dipole    schedule 06.09.2014
comment
См. Мой обновленный вопрос.   -  person Dipole    schedule 06.09.2014
comment
Так намного лучше.   -  person tacaswell    schedule 06.09.2014
comment
Спасибо, я буду обновлять больше, если добьюсь большего прогресса, и я ценю эту помощь!   -  person Dipole    schedule 06.09.2014
comment
Кроме того, часть проблемы может заключаться в том, что синий цвет лица на патче имеет более высокий z-oredr, чем ваше изображение, и, следовательно, скрывает его.   -  person tacaswell    schedule 06.09.2014
comment
Я попробовал это, установив для цвета значение «none», все равно не дало результата. Я не совсем уверен, как продолжить отладку экземпляра mainpoly - я читал документацию matplotlib.patches, но ее трудно понять. Прямо сейчас я просто экспериментирую со всеми методами, которые есть у моего объекта Polygon, возможно, это мне что-то скажет.   -  person Dipole    schedule 07.09.2014


Ответы (1)


Я тоже так долго думал об этой проблеме.
И я обнаружил, что язык NCL имеет функцию маскировки данных за пределами некоторой границы.
Вот пример:

http://i5.tietuku.com/bdb1a6c007b82645.png

Контурная диаграмма отображается только в пределах границы Китая. Щелкните здесь, чтобы просмотреть код.

Я знаю, что у python есть пакет PyNCL, который поддерживает весь код NCL в среде Python.
Но я действительно хочу построить такую ​​фигуру, используя базовую карту. Если вы в этом разобрались, опубликуйте сообщение в Интернете. Узнаю в первый раз.

Спасибо!

Добавить 2016-01-16

В каком-то смысле я понял это.
Это моя идея и код, и он вдохновлен тем вопросом, который я задал сегодня.

Мой метод:
1. Превратите шейп-файл интересующей области (например, США) в shapely.polygon.
2. Протестируйте каждую точку значения внутри / вне многоугольника.
3. Если точка значения равна вне области исследования, замаскируйте его как np.nan

Введение * многоугольник xxx был городом в Китае в формате шейп-файла ESRI. * fiona, здесь использовалась красивая упаковка.

# generate the shapely.polygon
shape = fiona.open("xxx.shp")
pol = shape.next()
geom = shape(pol['geometry'])
poly_data = pol["geometry"]["coordinates"][0]
poly = Polygon(poly_data)

Это выглядит так:

http://i4.tietuku.com/2012307faec02634.png

### test the value point 
### generate the grid network which represented by the grid midpoints.
lon_med  = np.linspace((xi[0:2].mean()),(xi[-2:].mean()),len(x_grid))
lat_med  = np.linspace((yi[0:2].mean()),(yi[-2:].mean()),len(y_grid))

value_test_mean = dsu.mean(axis = 0)
value_mask =  np.zeros(len(lon_med)*len(lat_med)).reshape(len(lat_med),len(lon_med))
for i in range(0,len(lat_med),1):
    for j in range(0,len(lon_med),1):
        points = np.array([lon_med[j],lat_med[i]])
        mask = np.array([poly.contains(Point(points[0], points[1]))])
        if mask == False:
            value_mask[i,j] = np.nan
        if mask == True:
            value_mask[i,j] = value_test_mean[i,j]


# Mask the np.nan value 
Z_mask = np.ma.masked_where(np.isnan(so2_mask),so2_mask)

# plot!
fig=plt.figure(figsize=(6,4))
ax=plt.subplot()

map = Basemap(llcrnrlon=x_map1,llcrnrlat=y_map1,urcrnrlon=x_map2,urcrnrlat=y_map2)
map.drawparallels(np.arange(y_map1+0.1035,y_map2,0.2),labels=  [1,0,0,1],size=14,linewidth=0,color= '#FFFFFF')
lon_grid  = np.linspace(x_map1,x_map2,len(x_grid))
lat_grid  = np.linspace(y_map1,y_map2,len(y_grid))
xx,yy = np.meshgrid(lon_grid,lat_grid)
pcol =plt.pcolor(xx,yy,Z_mask,cmap = plt.cm.Spectral_r ,alpha =0.75,zorder =2)

результат

http://i4.tietuku.com/c6620c5b6730a5f0.png

http://i4.tietuku.com/a22ad484fee627b9.png

исходный результат

http://i4.tietuku.com/011584fbc36222c9.png

person Han Zhengzu    schedule 11.12.2015
comment
Спасибо - я признаю, что давно не смотрел на это, но я дам вам знать, если что-то появится. - person Dipole; 11.12.2015
comment
Большой. Я тоже попробую. - person Han Zhengzu; 11.12.2015