基于Matplotlib与Cartopy正射投影绘图的两个技术问题咨询

问题描述 #

我希望绘制出类似示例的地图，但需移除绘图数据范围外的海岸线、边界线和网格线。第二个问题：能否将x、y坐标轴沿绘图数据的轮廓设置，使经纬度标签位于该轮廓旁而非方形边框旁？以下是我使用的代码：

import xarray as xr
from matplotlib import pyplot as plt
import cartopy.crs as ccrs
import cartopy

ds = xr.open_dataset(r'D:\Python\qq_z1000.nc')
ds2 = ds['qq'].mean(dim='time')
minn = ds2.min()
maxx = ds2.max()

fig = plt.figure(figsize=(8, 8), dpi=300, num="False")
central_lon, central_lat = 12.5, 47.5
extent = [-22.5, 45, 25, 65]
ax = plt.axes(projection=ccrs.Orthographic(central_lon, central_lat))
ax.set_extent(extent)

gl = ax.gridlines(draw_labels=True, linestyle='--', color='grey', dms=True, 
                  x_inline=False, y_inline=False, linewidth=0.3)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
gl.xlocator = plt.FixedLocator(range(-180, 181, 10))

ax.coastlines(resolution='50m')
ax.add_feature(cartopy.feature.BORDERS, edgecolor='black')
im = ax.pcolormesh(ds2.longitude, ds2.latitude, ds2, transform=ccrs.PlateCarree(), 
                   cmap='jet', vmin=minn, vmax=maxx)

cbar = fig.colorbar(im, orientation='horizontal', pad=0.05, shrink=0.8)

#plt.savefig('my_map2.jpg', dpi=300, bbox_inches='tight')
plt.show()

我尝试修改im = ax.pcolormesh(ds2.longitude, ds2.latitude, ds2, transform=ccrs.PlateCarree(), cmap='jet', vmin=minn, vmax=maxx)这行代码的多种写法，但均未解决问题。

问题1：移除数据范围外的海岸线、边界线和网格线 #

核心是通过数据的经纬度范围裁剪地理要素，只保留数据覆盖区域内的内容：

1. 提取数据有效范围：从ds2中获取经纬度的最大、最小值，确定数据覆盖的地理边界。
2. 裁剪海岸线与边界：遍历内置的shapefile要素，只添加落在数据范围内的地理图形。
3. 限制网格线范围：将网格线的定位器设置为数据范围内的经纬度刻度，避免绘制超出区域的网格。

对应的代码修改部分：

# 获取数据的经纬度范围
lon_min, lon_max = ds2.longitude.min().item(), ds2.longitude.max().item()
lat_min, lat_max = ds2.latitude.min().item(), ds2.latitude.max().item()

# 替换原网格线设置：只显示数据范围内的网格
gl = ax.gridlines(draw_labels=True, linestyle='--', color='grey', dms=True, 
                  x_inline=False, y_inline=False, linewidth=0.3)
gl.top_labels = False
gl.right_labels = False
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
# 只设置数据范围内的经纬度刻度
gl.xlocator = plt.FixedLocator(range(int(lon_min)//10*10, int(lon_max)//10*10+10, 10))
gl.ylocator = plt.FixedLocator(range(int(lat_min)//10*10, int(lat_max)//10*10+10, 10))

# 替换原海岸线和边界添加方式：只加载数据范围内的要素
from cartopy.io.shapereader import Reader

# 裁剪海岸线
coast_shp = cartopy.io.shapereader.natural_earth(resolution='50m', category='physical', name='coastline')
for geom in Reader(coast_shp).geometries():
    b = geom.bounds
    if b[0] <= lon_max and b[2] >= lon_min and b[1] <= lat_max and b[3] >= lat_min:
        ax.add_geometries([geom], ccrs.PlateCarree(), edgecolor='black', facecolor='none')

# 裁剪边界线
borders_shp = cartopy.io.shapereader.natural_earth(resolution='50m', category='cultural', name='admin_0_boundary_lines_land')
for geom in Reader(borders_shp).geometries():
    b = geom.bounds
    if b[0] <= lon_max and b[2] >= lon_min and b[1] <= lat_max and b[3] >= lat_min:
        ax.add_geometries([geom], ccrs.PlateCarree(), edgecolor='black', facecolor='none')

问题2：沿数据轮廓设置坐标轴并放置经纬度标签 #

Orthographic投影下无法直接将坐标轴绑定到数据轮廓，可通过以下两种方式实现近似效果：

方式1：用数据轮廓遮挡方形边框 #

通过生成数据的凸包路径，添加一个覆盖背景的补丁，隐藏原有的方形边框，只显示数据区域：

import numpy as np
from scipy.spatial import ConvexHull
from matplotlib.path import Path
from matplotlib.patches import PathPatch

# 生成数据掩码，提取有效经纬度点
mask = ~ds2.isnull()
lon_valid = ds2.longitude[mask.any(dim='latitude')].values
lat_valid = ds2.latitude[mask.any(dim='longitude')].values

# 创建数据凸包路径
points = np.column_stack((lon_valid, lat_valid))
hull = ConvexHull(points)
hull_path = Path(points[hull.vertices])

# 转换为投影坐标并添加补丁
proj_path = ax.projection.transform_path(hull_path, ccrs.PlateCarree())
patch = PathPatch(proj_path, facecolor='white', edgecolor='black', zorder=-1)
ax.add_patch(patch)

# 隐藏原有的方形边框
ax.spines['geo'].set_visible(False)

方式2：自定义经纬度标签位置 #

关闭原网格线标签，手动在数据轮廓的关键点添加标签：

# 关闭原网格线的标签
gl.draw_labels = False

# 在数据轮廓的关键位置添加经纬度标签
# 示例：在数据范围的四个角落添加
corner_points = [(lon_min, lat_min), (lon_max, lat_min), (lon_min, lat_max), (lon_max, lat_max)]
for lon, lat in corner_points:
    # 转换为投影坐标
    x, y = ax.projection.transform_point(lon, lat, ccrs.PlateCarree())
    # 添加经度标签
    ax.text(x, y-150000, f'{lon}°', ha='center', va='top', fontsize=10)
    # 添加纬度标签
    ax.text(x+150000, y, f'{lat}°', ha='left', va='center', fontsize=10)

完整修改代码 #

import xarray as xr
from matplotlib import pyplot as plt
import cartopy.crs as ccrs
import cartopy
from cartopy.io.shapereader import Reader
import numpy as np
from scipy.spatial import ConvexHull
from matplotlib.path import Path
from matplotlib.patches import PathPatch

ds = xr.open_dataset(r'D:\Python\qq_z1000.nc')
ds2 = ds['qq'].mean(dim='time')
minn = ds2.min()
maxx = ds2.max()
mask = ~ds2.isnull()

fig = plt.figure(figsize=(8, 8), dpi=300, num="False")
central_lon, central_lat = 12.5, 47.5
extent = [-22.5, 45, 25, 65]
ax = plt.axes(projection=ccrs.Orthographic(central_lon, central_lat))
ax.set_extent(extent)

# 获取数据经纬度范围
lon_min, lon_max = ds2.longitude.min().item(), ds2.longitude.max().item()
lat_min, lat_max = ds2.latitude.min().item(), ds2.latitude.max().item()

# 设置裁剪后的网格线
gl = ax.gridlines(draw_labels=False, linestyle='--', color='grey', dms=True, 
                  x_inline=False, y_inline=False, linewidth=0.3)
gl.xlocator = plt.FixedLocator(range(int(lon_min)//10*10, int(lon_max)//10*10+10, 10))
gl.ylocator = plt.FixedLocator(range(int(lat_min)//10*10, int(lat_max)//10*10+10, 10))

# 裁剪并添加海岸线
coast_shp = cartopy.io.shapereader.natural_earth(resolution='50m', category='physical', name='coastline')
for geom in Reader(coast_shp).geometries():
    b = geom.bounds
    if b[0] <= lon_max and b[2] >= lon_min and b[1] <= lat_max and b[3] >= lat_min:
        ax.add_geometries([geom], ccrs.PlateCarree(), edgecolor='black', facecolor='none')

# 裁剪并添加边界线
borders_shp = cartopy.io.shapereader.natural_earth(resolution='50m', category='cultural', name='admin_0_boundary_lines_land')
for geom in Reader(borders_shp).geometries():
    b = geom.bounds
    if b[0] <= lon_max and b[2] >= lon_min and b[1] <= lat_max and b[3] >= lat_min:
        ax.add_geometries([geom], ccrs.PlateCarree(), edgecolor='black', facecolor='none')

# 绘制数据
im = ax.pcolormesh(ds2.longitude, ds2.latitude, ds2, transform=ccrs.PlateCarree(), 
                   cmap='jet', vmin=minn, vmax=maxx)

# 添加数据轮廓补丁，隐藏方形边框
lon_valid = ds2.longitude[mask.any(dim='latitude')].values
lat_valid = ds2.latitude[mask.any(dim='longitude')].values
points = np.column_stack((lon_valid, lat_valid))
hull = ConvexHull(points)
hull_path = Path(points[hull.vertices])
proj_path = ax.projection.transform_path(hull_path, ccrs.PlateCarree())
patch = PathPatch(proj_path, facecolor='white', edgecolor='black', zorder=-1)
ax.add_patch(patch)
ax.spines['geo'].set_visible(False)

# 自定义经纬度标签
corner_points = [(lon_min, lat_min), (lon_max, lat_min), (lon_min, lat_max), (lon_max, lat_max)]
for lon, lat in corner_points:
    x, y = ax.projection.transform_point(lon, lat, ccrs.PlateCarree())
    ax.text(x, y-150000, f'{lon}°', ha='center', va='top', fontsize=10)
    ax.text(x+150000, y, f'{lat}°', ha='left', va='center', fontsize=10)

cbar = fig.colorbar(im, orientation='horizontal', pad=0.05, shrink=0.8)

#plt.savefig('my_map2.jpg', dpi=300, bbox_inches='tight')
plt.show()

内容的提问来源于stack exchange，提问作者krzyba