.shp文件与.nc文件的数据处理

代码如下：

import geopandas as gpd
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.font_manager import FontProperties
from tqdm import tqdm

# ================== 配置部分 ==================
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.weight'] = 'bold'
cn_font = FontProperties(fname='AaKuangPaiShouShu-2.ttf', size=10)
plt.rcParams['axes.unicode_minus'] = False

# ================== 文件路径 ==================
shp_path = r'ZG_S.shp'
nc_path = r'compound_events_2010.nc'


# ================ 数据预处理部分 ================
def load_polygons(path):
    """加载并返回闭合多边形列表（带坐标系处理）"""
    gdf = gpd.read_file(path)

    polygons = []
    for geom in gdf.geometry:
        if geom.geom_type == 'Polygon':
            coords = list(geom.exterior.coords)
            if coords[0] != coords[-1]:
                coords.append(coords[0])
            polygons.append(coords)
        elif geom.geom_type == 'MultiPolygon':
            for poly in geom.geoms:
                coords = list(poly.exterior.coords)
                if coords[0] != coords[-1]:
                    coords.append(coords[0])
                polygons.append(coords)
    return polygons


# 加载多边形（带坐标系转换）
polygons = load_polygons(shp_path)
# polygons = load_polygons(shp_path)[0:2]  # 取前两个测试
print(f"成功加载 {len(polygons)} 个多边形")

# ================ 数据处理部分 ================
# 读取NetCDF数据
ds = xr.open_dataset(nc_path)
dem = ds['compound_event'].squeeze()

# 维度校验与调整（确保lat在前）
if dem.dims != ('time', 'lat', 'lon'):
    dem = dem.transpose('time', 'lat', 'lon')
print(f"数据维度确认: {dem.dims}")

# 数据聚合（保留浮点型）
dem_data = dem.sum(axis=0)
# dem_data = dem_data.isel(lat=slice(0, 4), lon=slice(0, 3))    # 取一部分做测试

print(f"聚合后数据形状: {dem_data.shape}")

# 网格生成（匹配lat/lon顺序）
lon_grid, lat_grid = np.meshgrid(dem.lon.values, dem.lat.values, indexing='xy')
# lon_grid, lat_grid = np.meshgrid(dem.lon.values[0:3], dem.lat.values[0:4], indexing='xy')  # 测试
print(f"网格维度: lon_grid={lon_grid.shape}, lat_grid={lat_grid.shape}")

# ============== 掩膜生成核心修正 ==============
mask = np.zeros(lat_grid.shape, dtype=bool)  # 注意维度顺序

for poly in tqdm(polygons, desc="处理多边形"):
    if len(poly) < 3:
        continue

    poly_array = np.array(poly)
    # 坐标范围筛选（扩大0.1度容差）
    lon_min, lon_max = poly_array[:, 0].min() - 0.1, poly_array[:, 0].max() + 0.1
    lat_min, lat_max = poly_array[:, 1].min() - 0.1, poly_array[:, 1].max() + 0.1

    # 生成有效区域切片
    valid_lon = (lon_grid >= lon_min) & (lon_grid <= lon_max)
    valid_lat = (lat_grid >= lat_min) & (lat_grid <= lat_max)
    valid_area = valid_lon & valid_lat

    if not np.any(valid_area):
        continue

    # 提取有效点（注意xy顺序）
    points = np.column_stack([
        lon_grid[valid_area].ravel(),
        lat_grid[valid_area].ravel()
    ])

    # 创建路径并判断
    path = Path(poly_array[:, :2])
    poly_mask = path.contains_points(points)

    # 精准更新掩膜（修复形状不匹配问题）
    temp_mask = np.zeros_like(valid_area, dtype=bool)
    temp_mask[valid_area] = poly_mask
    mask |= temp_mask  # 合并掩膜

# ============== 掩膜应用优化 ==============
# 维度对齐（直接使用xarray广播）
dem_masked = dem_data.where(mask)  # xarray自动对齐维度
print(f"有效区域占比: {np.mean(mask) * 100:.2f}%")

# ================ 可视化部分 ================
fig, ax = plt.subplots(figsize=(14 / 2.54, 8 / 2.54), dpi=600)

# 精确绘图（注意lat/lon顺序）
mesh = ax.pcolormesh(
    dem_masked.lon,
    dem_masked.lat,
    dem_masked,
    cmap='YlGn',
    shading='auto',
    vmin=0,
    vmax=dem_masked.max(),
)

# 叠加多边形边界
for poly in polygons:
    x, y = np.array(poly).T
    ax.plot(x, y, '--k', linewidth=0.5, alpha=0.8)

# 颜色条设置
cbar = fig.colorbar(mesh, ax=ax, extend='both', pad=0.02)
cbar.set_label('Compound Event Frequency', weight='bold')

# 坐标轴优化
ax.set_xlabel('Longitude', fontweight='bold')
ax.set_ylabel('Latitude', fontweight='bold')
ax.set_aspect('auto')

# plt.savefig('中国1.tif', bbox_inches='tight', dpi=600)
plt.show()
plt.close()

print("处理完成，结果已保存")

效果如下：