一开始也是网上搜了很多方法,比起直接在网络上搜,直接在微信上搜方法(可返回相关公众号博主的一些方法)也是很有用的;

我尝试了两种方法,其中第一种是因为salem包没有安装成功(最终成功的方法是询问同学salem其它支撑包的版本,做以替换才可以),这个方法基本参考以下代码:

Python利用shp剪裁nc数据(11)

修改自己的参数后:

###批量nc裁剪
import os
os.environ['USE_PYGEOS']='0'
import geopandas as gpd
import numpy as np
import xarray as xr
import rioxarray
from shapely.geometry import mapping
import matplotlib.pyplot as plt

os.chdir('E:/seasonyield_predict/datapre/climate/pre')
filenames=glob.glob('*.nc')
#读取shp文件
shp_dir='E:/seasonyield_predict/datapre/region/clip/clip_dbc.shp'
db_shp=gpd.read_file(shp_dir)
outpath='E:/seasonyield_predict/datapre/climate/pre_cliped_nc/'

for file in filenames:
    data=xr.open_dataset(file)
    data=data.prep.loc[:,:,:]
    data=data.transpose('time','lon','lat')
    for i in range(len(data)):
        datab=data[i].transpose('lon','lat')
        datab.rio.write_crs("EPSG:4326", inplace=True)
        datab.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True)
        clipped=datab.rio.clip(db_shp.geometry.apply(mapping),db_shp.crs,drop=False)
        out=outpath+str(data[i].time.values)[0:10]+'_pre.nc'
        clipped.to_netcdf(out,mode='w', format="NETCDF4")
        print(out)

 

第二种方法用到salem包,主要参考下述代码:

一文学会python批量裁剪tif nc并可视化

修改自己的参数后:

import os
os.environ['USE_PYGEOS']='0'
import geopandas as gpd
import xarray as xr
import matplotlib.pyplot as plt
import salem
import glob
import numpy as np
from osgeo import gdal,osr,ogr
# from datetime import datetime,timedelta

os.chdir('E:/seasonyield_predict/datapre/PAR/par_nc')
filenames=glob.glob('*.nc')
#读取shp文件
shp_dir='E:/seasonyield_predict/datapre/region/clip4326/db.shp'
db_shp=gpd.read_file(shp_dir)
outpath='E:/seasonyield_predict/datapre/PAR/par_cliped_nc/'
minx = db_shp.bounds['minx'].min()
miny = db_shp.bounds['miny'].min()
maxx = db_shp.bounds['maxx'].max()
maxy = db_shp.bounds['maxy'].max()
lon_range = [minx, maxx] # 经度范围
lat_range = [miny, maxy] # 纬度范围 

for file in filenames:
    data=xr.open_dataset(filenames[0])
    data=data.scpdsi.loc['1984-01-01':'2013-12-31',:,:]
    for i in range(len(data)):
        maskregion=data[i].salem.roi(shape=db_shp)
        maskregion=maskregion.salem.subset(corners=((lon_range[0],lat_range[0]),(lon_range[1],lat_range[1])))
        out=outpath+str(data[i].time.values)[0:10]+'.nc'
        maskregion.to_netcdf(out,mode='w', format="NETCDF4")
        print(out)

其中遇到的问题:

1. 两种方法运行后,原来的nc文件范围是全国范围,用rio.clip和salem.rio两种方法裁剪后,显示出来的图像确实只有想要的矢量区范围,但是查看属性的时候它的实际范围还是全国,这样最后再转为tif的时候,还是全国的一个tif;询问了同学后,才明白可以在salem方法中指定范围,问题解决!

2.salem方法中,指定范围后裁剪显示错误,然后询问了同学,是矢量文件投影的问题,把原来的投影坐标系换成4326坐标系,解决问题!

关于坐标系的转换,参考:

坐标系不一致?GIS最全重投影方法

Logo

腾讯云面向开发者汇聚海量精品云计算使用和开发经验,营造开放的云计算技术生态圈。

更多推荐