文章目录
- 前言
- 实例数据
- 代码部分
- 需要的库
- 加载文件
- 写入地理信息
- 裁剪NC
- 结果
- 完整代码奉上
前言
Hello大家好呀,最近正好需要用到多个SHP去裁剪NC,按照我以前的两种办法(办法1和办法2)操作的话,我自己都会破防,在偶然的情况下发现了rioxarray这个库。这个库可以将坐标信息直接写入到NC之中,从而直接使用shp裁剪,比原先的方法更为简单明了。
实例数据
本次测试数据采用2m气温数据,下载链接:点我直接下载;
NC数据打开后如下图所示,该数据为逐月的三维数据。
SHP数据我就直接采用NOAA公开的海岸线数据吧,主要是涉及到中国的SHP极大可能性会被平台封文章。下载链接:点我直接下载。数据如下图所示,我用的文件是GSHHS_i_L1.shp
代码部分
需要的库
本次我们需要用到以下四个库,这几个库都可以直接用pip安装比较方便。
import xarray as xr
import rioxarray
import geopandas as gpd
from shapely.geometry import mapping
加载文件
ds = xr.open_dataset(r'C:\Users\Administrator\Downloads/air.2m.mon.mean.nc') # 加载NC文件
# ds['lon'] = ds['lon'] -180
shp = gpd.read_file(r"D:\BaiduNetdiskDownload\gshhg-shp-2.3.7\GSHHS_shp\i/GSHHS_i_L1.shp") # 加载shp文件
写入地理信息
ds.rio.write_crs("epsg:4326", inplace=True) # 向NC文件中写入WGS84坐标系信息
ds.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True) # 指定NC文件的xy识别名称
裁剪NC
ds = ds.rio.clip(shp.geometry.apply(mapping),shp.crs,drop=False)
结果
结果是出来了,但是我们会发现,数据只有0-180度之间的,而-180-0之间的数据没有了,造成这个问题的原因就是NCEP数据的lon是从0-360的,而这个库识别不出来。因此在裁剪之前我们需要仔细的观察一下数据的lon是否是从-180开始的。如果是从-180开始的就不用管了,如果是从0开始的就需要将加载文件那一步中的ds[‘lon’] = ds[‘lon’] -180取消注释。
经过lon处理后的结果就是正常的啦!!
完整代码奉上
import xarray as xr
import rioxarray
import geopandas as gpd
from shapely.geometry import mappingds = xr.open_dataset(r'C:\Users\Administrator\Downloads/air.2m.mon.mean.nc')
# ds['lon'] = ds['lon'] - 180 # 根据实际NC文件的经度范围确定是否开启这一行
shp = gpd.read_file(r"D:\BaiduNetdiskDownload\gshhg-shp-2.3.7\GSHHS_shp\i/GSHHS_i_L1.shp")ds.rio.write_crs("epsg:4326", inplace=True)
ds.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True)ds = ds.rio.clip(shp.geometry.apply(mapping),shp.crs,drop=False)ds.to_netcdf(r'D:\CSDN/air.2m.mon.mean_clip.nc') # 保存NC文件