python - 如何将 Matplotlib Basemap 的 maskoceans() 应用于多边形面片

Question

我想制作一张新西兰的地 block ，每个区域根据一些数据进行颜色编码。但是，我用来生成图像的 shapefile 超出了陆地边缘并进入了海洋。这意味着当我为多边形着色时，它们最终也会为海洋的一部分着色。不是想要的响应!

我使用的代码是:

from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon

plt.figure(figsize=(15,15))
lonmax = 180
lonmin = 165
latmax = -33
latmin = -48
map = Basemap(llcrnrlon=lonmin,llcrnrlat=latmin,urcrnrlon=lonmax,urcrnrlat=latmax, resolution = 'i')

map.drawmapboundary(fill_color='white')
map.fillcontinents(color='white',lake_color='white')
map.drawcoastlines()

map.readshapefile('../data/raw/statsnzregional-council-2016-generalised-version-SHP/regional-council-2016-generalised-version', 'regional_council')

ax = plt.gca() # get current axes instance
cm = matplotlib.cm.get_cmap('viridis')
norm = matplotlib.colors.Normalize(vmin=0.05, vmax=0.35)

for info, shape in zip(map.regional_council_info, map.regional_council):
    poly = Polygon(shape, facecolor=cm(norm(region_percent[info['REGC2016_N']])),edgecolor='k')
    ax.add_patch(poly)

plt.show()

产生下图。此图像非常接近我想要的，但我希望着色停止在陆地边界，而不是同时着色海洋。

我研究了 Basemap 的 maskoceans() 并相信这可能是解决这个问题的最佳方法，但我不明白如何将它应用到我的特定情况(例如，如何访问纬度、经度？什么是在这种情况下是数据数组？)

或者有没有办法使新西兰的 map 边界成为硬边界，以便仅打印与多边形补丁的内部重叠？

Answer 1

您需要一些多边形来遮盖多余的区域。在此处获取掩码文件(nz_mask_w.shp、nz_mask_e.shp): https://github.com/swatchai/cartopy_asean_proj/tree/master/shape_files这是代码:

import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import shapefile  # used to read my shapefiles

fig = plt.figure(figsize=(15,15))
lonmax = 179.95
lonmin = 165
latmax = -33
latmin = -48
map = Basemap(llcrnrlon=lonmin, \
              llcrnrlat=latmin, \
              urcrnrlon=lonmax, \
              urcrnrlat=latmax, \
              resolution = 'i')

# this is map theme (change to your need)
map.readshapefile('data/statsnzregional/regional-council-2016-generalised-version', \
                  name='regional_council')

ax = plt.gca()  # get current axes instance

#cm = matplotlib.cm.get_cmap('viridis')
#norm = matplotlib.colors.Normalize(vmin=0.05, vmax=0.35)

for info, shape in zip(map.regional_council_info, map.regional_council):
    poly = Polygon(shape, \
                facecolor=cm(norm(int(info['REGC2016']))/100.), \
                edgecolor='k', \
                zorder=1)  
    # original:- facecolor=cm(norm(info['REGC2016']))
    ax.add_patch(poly)


# mask out the excess areas (use files in data sub folder)
sf = shapefile.Reader("data/nz_mask_w")
ss = sf.shapes()
poly1 = Polygon(ss[0].points)

sf = shapefile.Reader("data/nz_mask_e")
ss = sf.shapes()
poly2 = Polygon(ss[0].points)

ax.add_collection( PatchCollection([poly1,poly2], \
                                   zorder=12, \
                                   facecolor='lightblue', \
                                   edgecolor='lightblue' ) )

map.drawcoastlines(color='blue', linewidth=0.3)
plt.show()