使用 shapefile 或 geopandas 绘制蒙面南极洲
Plotting a masked Antarctica with a shapefile or geopandas
我试图在掩盖大陆的同时绘制南极洲周围的数据。当我使用 basemap
并且它有一个选项可以使用 map.fillcontinents()
轻松掩盖大陆时,basemap
考虑的大陆包括我不想掩盖的冰架。
我尝试使用我在 Internet 上找到的代码中的 geopandas
。这是有效的,除了海岸线在我假设的南极洲多边形的 beginning/end 中产生了一条不需要的线:
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 之前的 和我编辑的 shapefile 生成了一个很好的蒙版 Antarctica/continents,但是海岸线超出了地图的圆形边缘:
我上传了 here the edited shapefile I used, but it's the Natural Earth Data 50m land shapefile 没有线。
这是一个如何实现您想要的示例。我基本上按照 Basemap example how to deal with shapefiles
and added a bit of 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()
结果如下所示:
希望对您有所帮助。
我试图在掩盖大陆的同时绘制南极洲周围的数据。当我使用 basemap
并且它有一个选项可以使用 map.fillcontinents()
轻松掩盖大陆时,basemap
考虑的大陆包括我不想掩盖的冰架。
我尝试使用我在 Internet 上找到的代码中的 geopandas
。这是有效的,除了海岸线在我假设的南极洲多边形的 beginning/end 中产生了一条不需要的线:
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 之前的
我上传了 here the edited shapefile I used, but it's the Natural Earth Data 50m land shapefile 没有线。
这是一个如何实现您想要的示例。我基本上按照 Basemap example how to deal with shapefiles
and added a bit of 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()
结果如下所示:
希望对您有所帮助。