如何在Python中将Cartopy地图上用户绘制的多边形转换为二进制掩码？

我来帮你搞定这个问题！你要实现的核心是把用户在Cartopy地图上画的多边形坐标，和格陵兰模拟数据的经纬度网格做匹配，判断每个网格点是否在多边形内部，从而生成带地理信息的二进制掩码（内部为1，外部为0）。下面一步步拆解实现：

1. 先捕获用户绘制的多边形坐标 #

首先得通过Matplotlib的事件处理，捕获用户的鼠标操作：左键添加多边形顶点，右键闭合多边形并触发掩码生成。这里要注意把屏幕坐标转换成真实的经纬度坐标，Cartopy的地理轴一定要用transData.inverted()来做坐标转换，不然点的位置会完全不对。

给你一段可直接复用的代码片段：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
from matplotlib.path import Path

# 全局变量存储多边形点和闭合状态
poly_points = []
closed = False

def handle_mouse_click(event):
    global poly_points, closed
    if event.button == 1 and not closed:  # 左键添加顶点
        # 把屏幕坐标转成经纬度
        lon, lat = event.inaxes.transData.inverted().transform((event.x, event.y))
        poly_points.append((lon, lat))
        # 实时绘制多边形预览
        if len(poly_points) > 1:
            x, y = zip(*poly_points)
            event.inaxes.plot(x, y, 'r-', linewidth=2, transform=ccrs.PlateCarree())
            event.canvas.draw()
    elif event.button == 3 and len(poly_points) >= 3 and not closed:  # 右键闭合多边形
        # 闭合多边形，把第一个点再加一遍
        poly_points.append(poly_points[0])
        x, y = zip(*poly_points)
        event.inaxes.plot(x, y, 'r-', linewidth=2, transform=ccrs.PlateCarree())
        event.canvas.draw()
        closed = True
        # 触发掩码生成函数
        generate_binary_mask(poly_points)

# 初始化格陵兰区域的Cartopy地图
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()}, figsize=(10, 8))
ax.coastlines(resolution='50m', color='black')
ax.set_extent([-70, -10, 50, 85])  # 格陵兰的大致经纬度范围
ax.gridlines(draw_labels=True)

# 绑定鼠标点击事件
cid = fig.canvas.mpl_connect('button_press_event', handle_mouse_click)
plt.show()

2. 生成匹配数据网格的二进制掩码 #

拿到闭合的多边形坐标后，就可以和你的格陵兰模拟数据的经纬度网格做匹配了。核心是用matplotlib.path.Path的contains_points方法，批量判断每个网格点是否在多边形内部。

如果你的模拟数据已经有lon_grid和lat_grid（二维网格），直接用下面的函数就行；如果是一维经纬度，先用np.meshgrid转成二维：

def generate_binary_mask(poly_points):
    # 这里替换成你自己的格陵兰模拟数据的经纬度
    # 示例：生成一个格陵兰区域的测试网格，实际用你自己的lon、lat数组
    lon = np.linspace(-70, -10, 200)
    lat = np.linspace(50, 85, 150)
    lon_grid, lat_grid = np.meshgrid(lon, lat)

    # 去掉多边形最后一个闭合点，创建Path对象
    poly_path = Path(poly_points[:-1])

    # 把二维网格转成Nx2的点数组，方便批量判断
    grid_points = np.column_stack((lon_grid.ravel(), lat_grid.ravel()))

    # 批量判断每个点是否在多边形内部，返回布尔数组
    is_inside = poly_path.contains_points(grid_points)

    # 把布尔数组转成和网格同形状的二进制掩码（1=内部，0=外部）
    mask = is_inside.reshape(lon_grid.shape).astype(int)

    # 现在mask就是和你的经纬度网格完全对应的二进制掩码了
    print(f"生成的掩码形状：{mask.shape}，和网格形状{lon_grid.shape}一致")

    # 可以把掩码和经纬度一起保存，方便后续和模拟数据结合
    np.savez('greenland_poly_mask.npz', lon=lon, lat=lat, mask=mask)

3. 几个要注意的细节 #

• 坐标系统必须一致：确保你的多边形坐标和模拟数据的坐标是同一个投影系统！比如都是WGS84经纬度（ccrs.PlateCarree()），如果你的模拟数据用的是极地投影（比如 stereographic），要先把多边形的经纬度转成数据的投影坐标再做判断，示例代码如下：

  # 假设你的数据用的是极地立体投影
  data_proj = ccrs.Stereographic(central_longitude=-45, central_latitude=70)
  # 把多边形的经纬度转成数据投影坐标
  poly_proj_coords = data_proj.transform_points(ccrs.PlateCarree(), *zip(*poly_points[:-1])).T[:2].T
  poly_path = Path(poly_proj_coords)
  # 同时你的数据网格也要是该投影下的坐标，或者把网格转成经纬度再判断
  

• 大网格的效率优化：如果你的网格特别大（比如2000x2000以上），contains_points可能有点慢，推荐用shapely库的矢量判断，速度会快很多：

  from shapely.geometry import Polygon
  from shapely.vectorized import contains
  
  # 创建shapely的Polygon对象
  poly = Polygon(poly_points[:-1])
  # 矢量判断每个网格点是否在内部
  mask = contains(poly, lon_grid, lat_grid).astype(int)
  

  用之前记得装shapely：pip install shapely
• Tkinter集成的注意点：如果你的程序是用Tkinter做GUI的，要确保Matplotlib用的是Tk后端，在初始化的时候加：

  import matplotlib
  matplotlib.use('TkAgg')
  

备注：内容来源于stack exchange，提问作者Nevpzo