python空间数据处理-part2:nc文件批量裁剪导出为tif

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

我尝试了两种方法，其中第一种是因为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最全重投影方法