无法在 cartopy 自定义形状边界和立体投影中添加网格线

Fail to add gridlines in cartopy custom shape boundary and stereographic projection

提问人:Jellywzx 提问时间:3/6/2023 最后编辑:python_userJellywzx 更新时间:3/22/2023 访问量:151

问:

我正在绘制自定义形状边界图,重点关注南大洋太平洋部分(160E~180~60W,-60S~-90S)。添加网格线时,区域(180~60W)缺少网格线。

Example output showing the problem

我尝试过许多解决方案(主要是矩形投影),但失败了。

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()
蟒蛇 Cartopy 网格线

评论


答:

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()

spole

评论

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()

map plotting result