用 shapefile 或 geopandas 绘制蒙面的南极洲

Posted 2023-03-23

【中文标题】用 shapefile 或 geopandas 绘制蒙面的南极洲

【英文标题】：Plotting a masked Antarctica with a shapefile or geopandas

【发布时间】：2018-10-26 23:15:17

【问题描述】：

我正在尝试绘制南极洲周围的数据，同时掩盖大陆。当我使用basemap 并且它可以选择使用map.fillcontinents() 轻松掩盖大陆时，basemap 考虑的大陆包括我不想掩盖的冰架。

我尝试使用我在 Internet 上找到的代码中的 geopandas。这行得通，除了海岸线在我认为是南极洲多边形的开始/结束时产生了一条不希望的线：

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.collections import PatchCollection
import geopandas as gpd
import shapely
from descartes import PolygonPatch

lats = np.arange(-90,-59,1)
lons = np.arange(0,361,1)
X, Y = np.meshgrid(lons, lats)
data = np.random.rand(len(lats),len(lons))

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

fig=plt.figure(dpi=150)
ax = fig.add_subplot(111)

m = Basemap(projection='spstere',boundinglat=-60,lon_0=180,resolution='i',round=True)  
xi, yi = m(X,Y)

cf = m.contourf(xi,yi,data)

patches = []
selection = world[world.name == 'Antarctica']
for poly in selection.geometry:
    if poly.geom_type == 'Polygon':
        mpoly = shapely.ops.transform(m, poly)
        patches.append(PolygonPatch(mpoly))
    elif poly.geom_type == 'MultiPolygon':
        for subpoly in poly:
            mpoly = shapely.ops.transform(m, poly)
            patches.append(PolygonPatch(mpoly))
    else:
        print(poly, 'blah')
ax.add_collection(PatchCollection(patches, match_original=True,color='w',edgecolor='k'))

当我尝试使用其他 shapefile（例如 land 可从Natural Earth Data 免费下载的 shapefile 时，会出现同一行。所以我在 QGIS 中编辑了这个 shapefile 来移除南极洲的边界。现在的问题是我不知道如何屏蔽 inside shapefile 的所有内容（也找不到如何做到这一点）。我还尝试通过设置linewidth=0 将前面的代码与geopandas 结合起来，并在顶部添加我创建的shapefile。问题是它们并不完全相同：

关于如何使用 shapefile 或使用 geopandas 但不使用线进行屏蔽的任何建议？

编辑：将 Thomas Khün 以前的 answer 与我编辑的 shapefile 一起使用会产生一个很好地掩蔽的南极洲/大陆，但海岸线超出了地图的圆形边缘：

我上传了here我使用的编辑过的shapefile，但它是Natural Earth Data 50m land shapefile没有那行。

【问题讨论】：

你为什么不使用你编辑的 shapefile 来掩盖南极洲？你能用你编辑的 shapefile 显示你用来绘制轮廓的代码吗？

@ThomasKühn 我不知道如何使用我编辑的 shapefile 来掩盖南极洲。我只是使用basemap's readshapefile:m.readshapefile('myshpfile_location', 'myshpfile_name',color='k',linewidth=0.6) 绘制 shapefile 的轮廓。它没有选项（或找不到），例如facecolor。

您能看看this 或this 是否对您有帮助吗？如果没有，请提供您的 shapefile，我可以尝试为您写一个答案。

@ThomasKühn 我编辑了帖子以显示当我尝试使用您的第一个答案时会发生什么。感谢您的意见！

【参考方案1】：

这里是一个如何实现你想要的例子。我基本上遵循Basemap example 如何处理shapefiles 并添加了一点shapely magic 以将轮廓限制在地图边界内。请注意，我首先尝试从ax.patches 中提取地图轮廓，但不知何故不起作用，因此我定义了一个半径为boundinglat 的圆，并使用底图坐标转换功能对其进行了转换。

import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.collections import PatchCollection
from matplotlib.patches import Polygon

import shapely
from shapely.geometry import Polygon as sPolygon

boundinglat = -40 
lats = np.arange(-90,boundinglat+1,1)
lons = np.arange(0,361,1)
X, Y = np.meshgrid(lons, lats)
data = np.random.rand(len(lats),len(lons))

fig, ax = plt.subplots(nrows=1, ncols=1, dpi=150)

m = Basemap(
    ax = ax,
    projection='spstere',boundinglat=boundinglat,lon_0=180,
    resolution='i',round=True
)  
xi, yi = m(X,Y)

cf = m.contourf(xi,yi,data)

#adjust the path to the shapefile here:
result = m.readshapefile(
    'shapefiles/AntarcticaWGS84_contorno', 'antarctica',
    zorder = 10, color = 'k', drawbounds = False)

#defining the outline of the map as shapely Polygon:
rim = [np.linspace(0,360,100),np.ones(100)*boundinglat,]
outline = sPolygon(np.asarray(m(rim[0],rim[1])).T)

#following Basemap tutorial for shapefiles
patches = []
for info, shape in zip(m.antarctica_info, m.antarctica):
    #instead of a matplotlib Polygon, create first a shapely Polygon
    poly = sPolygon(shape)
    #check if the Polygon, or parts of it are inside the map:
    if poly.intersects(outline):
        #if yes, cut and insert
        intersect = poly.intersection(outline)        
        verts = np.array(intersect.exterior.coords.xy)
        patches.append(Polygon(verts.T, True))

ax.add_collection(PatchCollection(
    patches, facecolor= 'w', edgecolor='k', linewidths=1., zorder=2
))

plt.show()

结果如下：

希望这会有所帮助。