Python处理NetCDF(nc)数据：从GDAL查看到rasterio写入GeoTIFF实战
本文详解如何使用Python处理NetCDF(nc)数据。首先通过GDAL命令行查看nc文件中的SUBDATASET信息，随后介绍利用netCDF4库读取数据为数组进行统计计算的方法。重点讲解使用rasterio库直接读取nc子数据集，并通过添加空间参考和变换信息将其转换为GeoTI格式，实现科学观测数据向地图可视化的完整流程。
前言
近期协助朋友处理了一批nc格式数据，原本以为操作简便，实际过程中却遇到了诸多细节与陷阱。无论是‘知难行易’还是‘知易行难’，似乎都无法完全概括其中的复杂性，‘知行合一’或许更为贴切——既要理解理论原理，也要掌握具体实现。经过初步研究，我将处理nc数据的完整步骤总结如下，供参考。
一、nc数据基础介绍
nc全称为NetCDF（The Network Common Data Form），其核心功能是存储一系列数组，结构相对简单（参考https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_introduction.html）。
由于nc格式擅长存储数组序列，它常被用于保存科学观测数据，尤其是长时间序列的数据。
假设科学家每隔一分钟采集一次实验数据并存储，若不使用nc格式，长期积累可能需要创建大量的csv或txt文件。而采用nc格式，单个文件即可容纳所有数据，显著提升了存储效率。
此外，若实验涉及气象、水文或温度等地理信息相关领域，nc格式同样适用。虽然GeoTiff能存储多个波段（可理解为气象、水文等不同信号），但nc格式能更好地存储不同波段的长时间观测结果，优势明显。
可以使用gdal工具查看数据信息，执行以下命令：
gdalinfo name.nc
执行后将显示如下信息：
Driver: netCDF/Network Common Data Format Files: test.nc Size is 512, 512 Coordinate System is `' Subdatasets: SUBDATASET_1_NAME=NETCDF:"test.nc":T2 SUBDATASET_1_DESC=[696x130x120] T2 (32-bit floating-point) SUBDATASET_2_NAME=NETCDF:"test.nc":PSFC SUBDATASET_2_DESC=[696x130x120] PSFC (32-bit floating-point) SUBDATASET_3_NAME=NETCDF:"test.nc":Q2 SUBDATASET_3_DESC=[696x130x120] Q2 (32-bit floating-point) SUBDATASET_4_NAME=NETCDF:"test.nc":U10 SUBDATASET_4_DESC=[696x130x120] U10 (32-bit floating-point) SUBDATASET_5_NAME=NETCDF:"test.nc":V10 SUBDATASET_5_DESC=[696x130x120] V10 (32-bit floating-point) SUBDATASET_6_NAME=NETCDF:"test.nc":RAINC SUBDATASET_6_DESC=[696x130x120] RAINC (32-bit floating-point) SUBDATASET_7_NAME=NETCDF:"test.nc":SWDOWN SUBDATASET_7_DESC=[696x130x120] SWDOWN (32-bit floating-point) SUBDATASET_8_NAME=NETCDF:"test.nc":GLW SUBDATASET_8_DESC=[696x130x120] GLW (32-bit floating-point) SUBDATASET_9_NAME=NETCDF:"test.nc":LAT SUBDATASET_9_DESC=[130x120] LAT (32-bit floating-point) SUBDATASET_10_NAME=NETCDF:"test.nc":LONG SUBDATASET_10_DESC=[130x120] LONG (32-bit floating-point) Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 512.0) Upper Right ( 512.0, 0.0) Lower Right ( 512.0, 512.0) Center ( 256.0, 256.0)
每个SUBDATASET代表一种格式的数据（如气象、水文等）。若要查看特定SUBDATASET的详细信息，可执行：
gdalinfo NETCDF:name.nc:SUBDATASET_NAME
此处的SUBDATASET_NAME对应上述输出中的T2、PSFC等字段，执行后可获得如下详细信息：
Driver: netCDF/Network Common Data Format Files: test.nc Size is 120, 130 Coordinate System is `' Metadata: LAT#description=LATITUDE, SOUTH IS NEGATIVE LAT#FieldType=104 LAT#MemoryOrder=XY LAT#stagger= LAT#units=degree_north Corner Coordinates: Upper Left ( 0.0, 0.0) Lower Left ( 0.0, 130.0) Upper Right ( 120.0, 0.0) Lower Right ( 120.0, 130.0) Center ( 60.0, 65.0) Band 1 Block=120x1 Type=Float32, ColorInterp=Undefined NoData Value=9.96920996838686905e+36 Unit Type: degree_north Metadata: description=LATITUDE, SOUTH IS NEGATIVE FieldType=104 MemoryOrder=XY NETCDF_VARNAME=LAT stagger= units=degree_north
此处仅包含一个Band，每个Band记录了一个时间点（或类似区分方式）的一条记录，该记录本质上是一个数组。
至此，读者应已明白，可以直接使用GDAL处理nc数据，例如直接使用gdalwarp将某个SUBDATASET转换为GeoTiff等格式。此处暂不展开，读者可查阅gdalwarp手册了解具体处理方法。
理解上述信息后，处理此类数据的基本思路便清晰了。
二、数据处理方法
Python应用广泛，其下的类库生态非常丰富。
2.1 netCDF4库
该框架可直接将nc文件读取为数组（详细信息参考https://github.com/Unidata/netcdf4-python（本地下载））。读取方式如下：
dataset = netCDF4.Dataset('name.nc') # open the dataset
如此即可读出整个nc文件中的数据信息。若需获取特定SUBDATASET，只需使用dataset[SUBDATASET_NAME]，返回的是一个三维数组，表示不同时间段（或其他区分方式）下的数据信息。
我们可以对该数组执行各种操作，如求平均值、方差等，这让人想起大学里那些枯燥却有趣的实验课程。当然，若使用numpy框架处理，效果更佳，例如求长时间序列的平均值：
np_arr = np.asarray(dataset[SUBDATASET_NAME]) average_arr = np.average(np_arr, axis=0)
至此，从事地理信息系统（GIS）工作的读者可能会发现一个问题：该框架仅能处理数据，无法进行与位置相关的操作，导致数据难以直接转化为地图可视化效果。其实数据是相通的，我们可以先以此方式处理，再转换为GeoTiff等格式；或者直接使用GeoTiff的处理流程进行操作。
2.2 rasterio库
rasterio是Mapbox开源的空间数据处理框架，功能强大，此处不赘述其所有功能，仅介绍如何处理我们的nc数据。
第一种方式是使用netCDF4处理完后，再用此框架写入GeoTiff，但这不够优雅，且涉及两个框架，略显繁琐。我们直接使用此框架从读取数据开始处理。
读取时需注意技巧，应像使用gdalinfo读取SUBDATASET一样，直接读取该SUBDATASET数据，如下：
with rio.open('NETCDF:name.nc:SUBDATASET_NAME') as src:
print(src.meta)
dim = int(src.meta['count'])
src.read(range(1, dim + 1))
即给open函数传入NETCDF:name.nc:SUBDATASET_NAME，采用src.read(range(1, dim + 1))可直接读出此范围内所有Band（时间点）的信息，范围可自行设定，注意从0开始，当然也可以仅读取某个Band的信息。
src.meta记录了此SUBDATASET的元数据信息，与gdalinfo看到的基本相同。
这样我们就可以继续将此数据使用numpy等框架进行处理，处理完之后更重要的是要写入GeoTiff中（直白的说就是添加空间信息）。
也很简单，如下即可：
with rio.open(newfile, 'w', **out_meta) as dst: dst.write_band(1, res_arr)
newfile为存储路径，res_arr为计算结果数组，注意尺寸不要发生变化（width*height），out_meta为目标文件的元数据描述信息，可以直接将上面src.meta进行简单处理即可。
out_meta =
meta.update({"driver": "GTiff",
"dtype": "float32",
'count': 1,
'crs': 'Proj4: +proj=longlat +datum=WGS84 +no_defs',
'transform': rasterio.transform.from_bounds(west, south, east, north, width, height)
})
crs表示目标数据空间投影信息，transform表示目标文件空间范围信息，可以通过经纬度信息和图像尺寸等计算得到。
dst.write_band将数据写入对应波段，当然此处也可以写入多个波段，根据计算结果而定，同样从1开始。
三、总结
本文简单介绍了nc数据的特点及如何使用python处理nc数据。每个目标都有多条路可以达到，重要的是找到那条自己喜欢的和适合自己的路，然而话又说回来，即使走的不是想要的那条路，不是一样可以达到目标嘛！所以关键是要找到自己的目标。
好了，以上就是这篇文章的全部内容了，希望本文的内容对大家的学习或者工作具有一定的参考学习价值，如果有疑问大家可以留言交流，谢谢大家对本站的支持。
免责声明：本文内容仅供参考，不代表本站观点。转载或引用请注明出处。
