当前位置：首页 > news >正文

利用Shp裁剪nc数据

news 来源：原创 2025/5/19 7:52:20

保证：

1、nc数据为等经纬度投影（不是等经纬度先转换为等经纬度）

2、Shp文件必须只有最外围（如果要剪裁京津冀，需要先得到该地区最外围的shp，不能包含内部的区域的shp），且shp必须为polygon，不能是polyline

剪裁法1：

import xarray as xr
import geopandas as gpd
from shapely.geometry import mapping
import rioxarray# 1. 创建插值后的 DataArray，使用纬度和经度作为坐标
da_interp = xr.DataArray(z_target_grid,dims=("lat", "lon"),coords={"lat": our_lats[:, 0], "lon": our_lons[0, :]},name="SWDOWN2"
)# 2. 设置空间维度和坐标参考系
da_interp.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True)
da_interp.rio.write_crs("EPSG:4326", inplace=True)# 3. 强制 shapefile 转为 EPSG:4326
shp = gpd.read_file(r'D:\ZZZZorder_jobs\250508\AAA\shp_Nor_China\test6.shp')
shp = shp.to_crs("EPSG:4326")
clipped = da_interp.rio.clip(shp.geometry.apply(mapping), shp.crs, drop=True)

绘图看看

import matplotlib.pyplot as plt# 创建子图
fig, axs = plt.subplots(1, 2, figsize=(14, 6))# 剪裁前的数据绘图
da_interp.plot(ax=axs[0], cmap='viridis')
axs[0].set_title("Before Clipping")# 剪裁后的数据绘图
clipped.plot(ax=axs[1], cmap='viridis')
axs[1].set_title("After Clipping")plt.suptitle("Comparison of Data Before and After Clipping")
plt.tight_layout()
plt.show()

如下图：

剪裁法2

import rasterio.features
import numpy as np
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import mapping# 读取 shapefile 并确认投影
shp = gpd.read_file(r'D:\ZZZZorder_jobs\250508\AAA\shp_Nor_China\test6.shp',crs="EPSG:4326")# 构建掩膜
mask = rasterio.features.geometry_mask([mapping(geom) for geom in shp.geometry],transform=da_interp.rio.transform(),out_shape=(da_interp.sizes['lat'], da_interp.sizes['lon']),invert=True  # 保留 geometry 内的区域为 True
)# 应用掩膜
masked = da_interp.where(mask)

绘图看看：

fig, axs = plt.subplots(1, 2, figsize=(18, 6))# 原始图
da_interp.plot(ax=axs[0], cmap='viridis')
axs[0].set_title('Before Masking')# 掩膜后图
masked.plot(ax=axs[1], cmap='viridis')
axs[1].set_title('After Masking by Shape')plt.suptitle('Comparison of Data Before and After Masking', fontsize=16)
plt.tight_layout()
plt.show()