首页 文章

如何将Matplotlib Basemap的maskoceans()应用于多边形补丁

提问于
浏览
4

我想制作一个新西兰的地块,根据一些数据对每个地区进行颜色编码 . 然而,我用来制作图像的形状文件延伸到土地边缘之外并进入海洋 . 这意味着当我为多边形着色时,它们最终也会对海洋的部分进行着色 . 不是想要的回应!

我使用的代码是:

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(),并且相信这可能是解决这个问题的最好方法,但是我不明白如何将它应用到我的特定情况(例如,如何访问lat,lons?什么是数据数组?这个案例?)

或者有没有办法让新西兰的 Map 边界成为一个硬边界,所以只打印内部与多边形补丁重叠?

enter image description here

1 回答

  • 2

    你需要一些多边形来掩盖多余的区域 . 在这里获取掩码文件(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()
    

    enter image description here

相关问题