将多个GeoTIFF格式的栅格时间序列转换为NetCDF格式

Question

我有多个*.tif格式的栅格时间序列文件，希望将它们合并成一个NetCDF文件。数据类型为uint16。

我可以使用gdal_translate命令将每个图像转换为netcdf格式:

 gdal_translate -of netcdf -co FORMAT=NC4 20150520_0164.tif foo.nc

然后使用NCO进行一些脚本编写，从文件名中提取日期并进行连接，但我想知道是否可以使用Python更有效地完成这项工作，使用xarray及其新的rasterio后端。

我可以轻松读取文件：

import glob
import xarray as xr
f = glob.glob('*.tif')
da = xr.open_rasterio(f[0]) 
da

返回

<xarray.DataArray (band: 1, y: 5490, x: 5490)>
[30140100 values with dtype=uint16]
Coordinates:
  * band     (band) int64 1
  * y        (y) float64 5e+05 5e+05 5e+05 5e+05 5e+05 4.999e+05 4.999e+05 ...
  * x        (x) float64 8e+05 8e+05 8e+05 8e+05 8.001e+05 8.001e+05 ...
Attributes:
    crs:      +init=epsg:32620

我可以将其中之一写入NetCDF：

ds.to_netcdf('foo.nc')

但理想情况下，我希望能够使用类似于xr.open_mfdataset的方法，写入时间值（从文件名中提取），然后将整个聚合写入netCDF。并且让dask处理内存问题。

请问xarray和dask能够做到这样的吗？

- Rich Signell

Answer 1

Xarray应该能够为您完成连接步骤。我稍微修改了您的示例如下。将文件名解析成有用的内容是由您自己决定的。

import glob
import pandas as pd
import xarray as xr

def time_index_from_filenames(filenames):
    '''helper function to create a pandas DatetimeIndex
       Filename example: 20150520_0164.tif'''
    return pd.DatetimeIndex([pd.Timestamp(f[:8]) for f in filenames])

filenames = glob.glob('*.tif')
time = xr.Variable('time', time_index_from_filenames(filenames))
chunks = {'x': 5490, 'y': 5490, 'band': 1}
da = xr.concat([xr.open_rasterio(f, chunks=chunks) for f in filenames], dim=time)

- jhamman

---

{btsdaf} - Rich Signell

{btsdaf} - jhamman

嗨@RichSignell，你介意和我分享你的脚本吗？我有1981年以来的月降水数据，想将其转换为带有时间维度的单个NetCDF文件。 - user97103