问题描述
我有一个很大的全局.nc文件数据集,我试图将它们裁剪到较小的区域。我将此区域存储为.shp文件。
我尝试使用来自Qgis的gdal,但是需要通过转换每个变量来做到这一点,我必须为所有文件一个一个地选择每个变量和相同的形状,并且要通过400个文件,每个变量似乎不是最好的主意。还会返回分隔的.tiff文件,而不是我要的.nc文件。
我有这个小脚本,但是没有满足我的要求
import glob
import subprocess
import os
ImageList = sorted(glob.glob('*.nc'))
print('number of images to process: ',len(ImageList))
Shapefile = 'NHAF-250m.shp'
# Create output directory
OutDir = './Clipped_Rasters/'
if not os.path.exists(OutDir):
os.makedirs(OutDir)
for Image in ImageList:
print('Processing ' + Image)
Outimage = OutDir + Image.replace('.nc','_BurnedArea_Clipped.tif') # Defines Output Image
# Clip image
subprocess.call('gdalwarp -q -cutline /Users/path/to/file/NHAF-250-vector/ -tr 0.25 0.25 -of GTiff NETCDF:'+Image+":burned_area "+Outimage,shell=True)
print('Done.' + '\n')
print('All images processed.')
提前谢谢
解决方法
我建议使用xarray
处理netcdf数据,并使用geopandas
+ rasterio
处理Shapefile。
import geopandas
import xarray
import rasterio
import glob
shapefile = 'NHAF-250m.shp'
sf = geopandas.read_file(shapefile)
shape_mask = rasterio.features.geometry_mask(sf.iloc[0],out_shape=(len(ndvi.y),len(ndvi.x)),transform=ndvi.geobox.transform,invert=True)
shape_mask = xarray.DataArray(shape_masj,dims=("y","x"))
file_list = sorted(glob.glob('*.nc'))
for file in file_list:
nc_file = xarray.open_dataset(file)
# Then apply the mask
masked_netcdf_file = nc_file.where(shape_mask == True,drop=True)
# store again as netcdf or do what every you want with the masked array