提问人:Jellywzx 提问时间:3/6/2023 最后编辑:python_userJellywzx 更新时间:3/22/2023 访问量:151
无法在 cartopy 自定义形状边界和立体投影中添加网格线
Fail to add gridlines in cartopy custom shape boundary and stereographic projection
问:
我正在绘制自定义形状边界图,重点关注南大洋太平洋部分(160E~180~60W,-60S~-90S)。添加网格线时,区域(180~60W)缺少网格线。
我尝试过许多解决方案(主要是矩形投影),但失败了。
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import cartopy
import cartopy.crs as ccrs
import numpy as np
proj = ccrs.Stereographic(central_longitude=228, central_latitude=-70)
fig = plt.figure(figsize=(10,10))
ax = plt.axes([0,0,1,1],projection=proj)
latmin = -90
latmax = -60
lonmin = 150
lonmax = 305
lats = np.linspace(latmax, latmin, latmax - latmin + 1)
lons = np.linspace(lonmin, lonmax, lonmax - lonmin + 1)
vertices = [(lon, latmin) for lon in range(lonmin, lonmax + 1, 1)] + \
[(lon, latmax) for lon in range(lonmax, lonmin - 1, -1)]
boundary = mpath.Path(vertices)
ax.set_boundary(boundary, transform=ccrs.PlateCarree())
ax.set_extent([170, 320, -90, -60], ccrs.PlateCarree())
ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='k')
ax.gridlines(draw_labels=True)
plt.show()
答:
2赞
swatchai
3/6/2023
#1
Gridlines
在这种特殊情况下,绘图是有问题的。该图涉及特殊的边界,即跨越国际日期变更线的纬线——这些可能是问题的原因。我希望@Rutger_Kassies的答案是你需要的。但是,如果您仍然需要原始绘图范围,这里有另一种方法。
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import cartopy
import cartopy.crs as ccrs
import numpy as np
proj = ccrs.Stereographic(central_longitude=228, central_latitude=-70)
fig = plt.figure(figsize=(10/2, 10/2))
ax = plt.axes([0,0,1,1], projection=proj)
latmin = -90
latmax = -60
lonmin = 150
lonmax = 305
lats = np.linspace(latmax, latmin, latmax - latmin + 1)
lons = np.linspace(lonmin, lonmax, lonmax - lonmin + 1)
vertices = [(lon, latmin) for lon in range(lonmin, lonmax + 1, 1)] + \
[(lon, latmax) for lon in range(lonmax, lonmin - 1, -1)]
boundary = mpath.Path(vertices)
ax.set_boundary(boundary, transform=ccrs.PlateCarree())
ax.set_extent([170, 320, -90, -60], ccrs.PlateCarree())
ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='k')
#ax.add_feature(cartopy.feature.OCEAN, zorder=1, edgecolor='none', alpha=0.3)
# Build gridlines in 2 steps
# First set of gridlines: meridians only
# proper `xlim` is needed to get all the x-labels' visibility set correctly
gl1 = ax.gridlines(draw_labels=True, crs=ccrs.PlateCarree(), \
xlocs=range(-180, 359, 20), \
ylocs=[], \
xlim=[-2847619, 3537282]) # get these from ax.get_xbound() of previous run
# Second set of gridlines: parallels of latitude only.
# The lines are terminated at the meridian of the dateline!
gl2 = ax.gridlines(draw_labels=True, crs=ccrs.PlateCarree(), \
y_inline=True, \
ylocs=range(-80, -50, 10), \
xlocs=[])
# Draw the missing parts of the parallel lines at 70, 80 deg_S
# Colors are set to `red` and `green` intentionally
lons = range(-180, -50, 5)
ax.plot(lons, -80*np.ones(len(lons)), transform=ccrs.Geodetic(), lw=0.5, color="red", zorder=30)
ax.plot(lons, -70*np.ones(len(lons)), transform=ccrs.Geodetic(), lw=0.5, color="green", zorder=30)
# Print the bounds of the plot
print("Bounds:", ax.get_xbound(), ax.get_ybound())
plt.show()
评论
0赞
swatchai
4/28/2023
另请查看此链接:github.com/SciTools/cartopy/issues/697#issuecomment-157025400 查看 pelson 的评论以获取另一种方法。
0赞
Stanislav Ivanov
3/22/2023
#2
在回答中,https://stackoverflow.com/a/65690841/4265407 描述了另一种边界投影方法,该方法允许网格正确绘制。
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import cartopy
import cartopy.crs as ccrs
proj = ccrs.Stereographic(central_longitude=228, central_latitude=-70)
fig = plt.figure(figsize=(15,8))
ax = plt.axes([0,0,1,1],projection=proj)
ylim = [-90,-60]
xlim = [150,305]
rect = mpath.Path([[xlim[0], ylim[0]],
[xlim[1], ylim[0]],
[xlim[1], ylim[1]],
[xlim[0], ylim[1]],
[xlim[0], ylim[0]],
]).interpolated(40)
proj_to_data = ccrs.PlateCarree()._as_mpl_transform(ax) - ax.transData
rect_in_target = proj_to_data.transform_path(rect)
ax.set_boundary(rect_in_target)
ax.set_extent([150,305,-90,-58], ccrs.PlateCarree())
ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='k')
ax.gridlines(draw_labels=True, x_inline=False, y_inline=False,
xlocs=range(-180, 180, 10))
plt.show()
评论