python - 通过 shapefile 剪切 NetCDF 文件
问题描述
我有一个大型的全局 .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
推荐阅读
- javascript - ESLINT:在评论块空格/制表符内进行 Linting
- apache-kafka - Spring 集成中的 MutableMessageBuilderFactory
- python - 模块 face_recognition 将实时网络摄像头保存到 TXT 文件
- javascript - Array.concat 反映连接后修改输入嵌套数组时结果数组的变化
- mysql - 如何在sequelize中生成具有多个或条件的查询
- r - 查找函数调用的物种特定系数的有效方法
- asp.net - 从 HTML 输入转换为字符串
- javascript - 鼠标悬停事件的JS问题以显示额外的按钮
- matlab - 如何在一个图中绘制多条线并保存它们
- r - 我想使用 bibliometrix 在 R 中转换一个 excel 文件