Python实现法国境内经纬度标记及NetCDF生成问题求助

问题背景 #

已创建北半球0.5°分辨率的全0二维经纬度数组，需要将法国境内的坐标点值改为1，最终生成包含经度、纬度和二进制值的NetCDF文件，用于后续与清单数据相乘得到法国区域清单。但当前代码处理法国多边形时失效，仅用4组坐标测试时正常。

代码问题分析 #

1. 几何对象类型错误：france['geometry']返回的是GeoSeries（带索引的几何集合），而非单个几何对象。调用contains(point)时会返回布尔数组，而非单个布尔值，导致判断逻辑失效。
2. 双重循环效率极低：遍历全量经纬度点的双重循环不仅速度慢，还容易因索引对应错误导致结果异常。
3. 未显式处理MultiPolygon：法国的几何形状是MultiPolygon（包含本土及海外领地），直接用GeoSeries判断会出现逻辑偏差。

修正方案 #

步骤1：提取单个几何对象 #

通过iloc[0]从GeoSeries中取出法国的实际几何对象（MultiPolygon/Polygon）。

步骤2：批量生成并判断点 #

用np.meshgrid生成所有经纬度点的坐标对，再通过shapely的vectorized.contains批量判断点是否在几何内，替代低效的双重循环。

步骤3：生成NetCDF文件 #

使用xarray库将坐标和掩码数组整合为Dataset，直接导出为NetCDF格式，操作更简洁。

完整修正代码 #

import numpy as np
import geopandas as gpd
import xarray as xr
from shapely.vectorized import contains

# 生成经纬度数组
lon = np.arange(-179.75, 180.25, 0.5)
lat = np.arange(0.25, 89.25, 0.5)
# 创建全0掩码数组（注意维度：lat在前，lon在后，符合NetCDF常规存储）
tab = np.zeros((len(lat), len(lon)))

# 读取国家边界 shp 文件
shapefile = gpd.read_file("./TM_WORLD_BORDERS-0.3.shp")

# 筛选法国的几何对象
france_geo = shapefile[shapefile['NAME'] == 'France'].geometry.iloc[0]

# 生成所有经纬度点的网格
lon_grid, lat_grid = np.meshgrid(lon, lat)
# 批量判断每个点是否在法国境内，返回布尔数组
mask = contains(france_geo, lon_grid.flatten(), lat_grid.flatten())
# 将布尔数组转为整数并重塑为原网格形状
tab = mask.reshape(tab.shape).astype(int)

# 生成NetCDF文件
ds = xr.Dataset(
    data_vars={'mask': (['lat', 'lon'], tab)},
    coords={'lon': lon, 'lat': lat}
)
# 导出为NetCDF
ds.to_netcdf('france_mask.nc')

# 验证结果（可选）
print("掩码数组中值为1的点数量：", np.sum(tab))

关键说明 #

• 维度顺序调整：将原代码中(len(lon), len(lat))改为(len(lat), len(lon))，符合NetCDF数据的常规存储规范（纬度为第一维度），避免后续与清单数据相乘时出现维度不匹配问题。
• 批量判断优化：shapely.vectorized.contains是向量化操作，比双重循环效率提升数十倍，尤其适合大分辨率网格。
• MultiPolygon兼容：contains方法天然支持MultiPolygon，会自动判断点是否属于其中任意一个子多边形。

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