将 TIFF 合成 NetCDF

引言

在科学数据处理中，我们大多数将影像存为 TIFF 格式，而忽视了NetCDF (nc)格式。其实nc文件非常利于数据传播与共享。

合成 NetCDF 格式的文件有许多好处，特别是在地理信息系统（GIS）和气象学等领域。以下是一些主要的好处：

多维数据支持：NetCDF 格式允许存储多维数据，包括时间序列、三维空间数据等。这对于存储和处理复杂的科学数据集非常有用。
元数据支持：NetCDF 具有内置的元数据支持，可以包括数据的描述、单位、坐标轴信息等。这使得数据更容易被他人理解和使用。
数据压缩：NetCDF 支持数据压缩，可以显著减小文件大小，降低存储成本，并提高数据传输效率。
数据子集和切片：NetCDF 允许按需读取数据的子集或切片，这对于处理大型数据集时非常有用，可以减少内存占用。
跨平台兼容性：NetCDF 格式是一种开放标准，跨平台兼容，可以在不同操作系统和软件环境中使用。

利用netCDF4包

接下来讲述用netCDF4包将TIF转为nc。

有读者可能会好奇，为什么选择 netCDF4 包而不是 xarray 包？xarray包明明是比较新的解决方案。

虽然 xarray 包是一个功能强大的工具，用于处理和分析多维数据，但在将 TIFF 合成 NetCDF 时，使用 netCDF4 包更为合适，原因如下：

直接控制：netCDF4 包提供了更直接的控制，允许用户以更细粒度的方式操作 NetCDF 文件的细节。这对于需要高度自定义数据存储结构的应用程序非常有用。
更低级别访问：netCDF4 包提供了更低级别的文件访问，适用于对 NetCDF 文件的底层操作。这是利用netCDF4的核心原因。因为netCDF4是向上兼容的，而用xarray写nc文件则有可能在netCDF4中打不开。
性能优化：netCDF4 包可以用于性能优化，尤其是在处理大型数据集时。用户可以选择不同的压缩策略、块大小等参数以提高数据读取和写入速度。
库的广泛应用：netCDF4 是一个广泛应用的库，许多科学软件和工具都支持该库，这使得在不同环境中交换和共享数据更为容易。

虽然 xarray 是高级工具，但在某些情况下，使用 netCDF4 包可以更好地满足特定需求，特别是在需要更多控制和性能优化时。

在这之前，假设你的数据是三维矩阵。（numpy中经常处理的结果，这时将其写为nc）

如果你的TIF是单张的，根据下面的代码把TIF转为三维矩阵：

image-20230915234458225

import os   import xarray as xr   import rioxarray as rxr   import numpy as np   from osgeo import gdal         def tif_to_3darray(folder_path):   # 获取文件夹中所有tif文件的路径       file_extension = '.tif'       file_paths = [os.path.join(folder_path, file) for file in os.listdir(folder_path) if file.endswith(file_extension)]       print([i.split('\\')[1] for i in file_paths])       dataset = gdal.Open(file_paths[1])       cols=dataset.RasterXSize#图像长度       rows=dataset.RasterYSize#图像宽度       num = len(file_paths)       array = np.empty((rows, cols, 0))              # 逐个读取tif文件并将数据添加到数据集中       for file_path in file_paths:           # 从tif文件创建xarray数据集           da = xr.open_dataset(file_path)           data = da['band_data'][0]           array = np.concatenate((array, np.expand_dims(data, axis=2)), axis=2)                     # 打印合成的三维数组的形状       print("合成的三维数组形状:", array.shape)       return array      trans5 = tif_to_3darray('v6_result')   test = rxr.open_rasterio('A201707.tif')

接下来是写入nc的核心代码，如下：

#w表示写入   data_NC.close()   # 关闭文件，当报错后一定要单独运行一下这一行，把nc文件关闭，不然他一直处于打开状态，会报错，然后找到#当前nc文件删除，重新运行   import pandas as pd   import netCDF4 as nc   data_NC = nc.Dataset('drop_data.nc', 'w', format='NETCDF4')           # 创建维度，createDimension（维度名称，维度长度）   data_NC.createDimension('lat', test.values[0].shape[0])   data_NC.createDimension('lon',test.values[0].shape[1])   data_NC.createDimension('time',trans5.shape[2])       # 创建维度变量，createVariable（变量名，值类型，维度）注意这里的维度就是上面创建维度的名称，不然会#报错   data_NC.createVariable("lat", 'f', ('lat'))   data_NC.createVariable("lon", 'f', ('lon'))   data_NC.createVariable("time", 'f', ('time'))       # 创建变量，("对应的函数值", '数据类型', ( "变量一", "变量二"))我创建的是二维数组   data_NC.createVariable("data", 'f', ("time", "lat", "lon"))       # 给维度填充数据，这里的数据是需要自己提前准备好，一般是linspace   # data_NC.variables['维度名'][:] = 维度数据       data_NC.variables['lat'][:] = test.y.values   data_NC.variables['lon'][:] = test.x.values   data_NC.variables['time'][:] = pd.date_range(start='2017-07', end='2018-06', freq='M').strftime('%Y%m').tolist()       # 给变量填充数据，和维度填充数据方式是一样的，就是存进去   for i in range(trans5.shape[2]):       data_NC.variables['data'][i, :, :] = trans5[:,:,i]       # 关闭文件，当报错后一定要单独运行一下这一行，把nc文件关闭，不然他一直处于打开状态，会报错，然后找到#当前nc文件删除，重新运行   data_NC.close()

data_NC.close() 这一行代码用于关闭之前打开的NetCDF文件，确保文件被正确关闭，以防止后续的操作发生错误。
创建维度，使用 createDimension 函数创建了三个维度：'lat'、'lon' 和 'time'。这些维度用于描述数据的形状和组织结构。
创建维度变量，使用 createVariable 函数创建了与上面创建的维度对应的维度变量：'lat'、'lon' 和 'time'。
填充维度数据，使用 data_NC.variables['维度名'][:] = 维度数据 的方式填充了维度数据。具体来说，'lat' 和 'lon' 维度的数据来自于已有tif，否则用linspace生成等间隔经纬度。而 'time' 维度的数据是一个日期范围，这里用pd.date_range生成
使用循环，通过 data_NC.variables['data'][i, :, :] = trans5[:,:,i] 的方式填充了数据变量 'data'。这个循环将三维数据存储到 'data' 变量中的每个时间。

读取结果

可以用xarray包读取写成的nc，查看是否结果正确

import xarray   import salem   file_name='drop_data.nc'   ds=xarray.open_dataset(file_name, engine='netcdf4')

ds['data'][0]

读取成功，接下来使用之前写的函数进行可视化，感兴趣的小伙伴可以看相关推文（文末）。

from utils import plot   import cartopy.crs as ccrs   import rioxarray as rxr   import matplotlib.pyplot as plt   import numpy as np   import salem   import geopandas as gpd   shpfile=gpd.read_file('D:/OneDrive/data/countries.shp')

绘制全球地图和统计线：

fig = plt.figure()   proj = ccrs.PlateCarree() #ccrs.Robinson()ccrs.Mollweide()Mollweide()   ax = fig.add_axes([0.1, 0.1, 0.9, 0.9], projection=proj)   levels = np.linspace(-30, 30, num=19)   plot.one_map_flat(ds['data'][0], ax, levels=levels, cmap="RdBu",  mask_ocean=True, add_coastlines=True, add_land=False,  colorbar=True, plotfunc="pcolormesh")   ax.set_ylim([-60, 90])   ax.set_title("Our gap-filling 2017-07", fontsize=15, pad=8)   ax2 = fig.add_axes([1.05, 0.3, 0.15, 0.5])   plot.add_sta(ax2, ds['data'][0].salem.roi(shape=shpfile), [-80,20], 'lat')

利用xarray

xarray写nc相对简单，但不容易向下兼容，因此不推荐。以下是一个代码示例。

temp = 15 + 8 * np.random.randn(2, 2, 3)   precip = 10 * np.random.rand(2, 2, 3)   lon = [[-99.83, -99.32], [-99.79, -99.23]]   lat = [[42.25, 42.21], [42.63, 42.59]]      ds = xr.Dataset(       {           "temperature": (["x", "y", "time"], temp),           "precipitation": (["x", "y", "time"], precip),       },       coords={           "lon": (["x", "y"], lon),           "lat": (["x", "y"], lat),           "time": pd.date_range("2014-09-06", periods=3),           "reference_time": pd.Timestamp("2014-09-05"),       },   )