GDAL读取HDF、NetCDF数据集

ytkz2023-09-182024-09-09

HDF5、NetCDF

默认HDF5、NetCDF这两类文件是相似，读取文件的步骤也是相似的。

HDF5的文件后缀一般是.h5，NetCDF的文件后缀一般是.nc

以下用nc文件代指HDF、NetCDF文件。

对程序员来说，nc文件和zip、jpeg、bmp文件格式类似，都是一种文件格式的标准。netcdf文件开始的目的是用于存储气象科学中的数据，已经成为许多数据采集软件的生成文件的格式。利用NetCDF可以对网格数据进行高效地存储、管理、获取和分发等操作。

由于其灵活性，能够传输海量的面向阵列（array-oriented）数据，广泛应用于大气科学、水文、海洋学、环境模拟、地球物理等诸多领域，例如，NCEP(美国国家环境预报中心)发布的再分析资料，NOAA的CDC(气候数据中心)发布的海洋与大气综合数据集(COADS)均采用NetCDF作为标准。

从数学上来说，netcdf存储的数据就是一个多自变量的单值函数。用公式来说就是f(x,y,z,…)=value,函数的自变量x,y,z等在netcdf中叫做维(dimension)或坐标轴(axis),函数值value在netcdf中叫做变量(Variables).而自变量和函数值在物理学上的一些性质，比如计量单位(量纲)、物理学名称等等在netcdf中就叫属性(Attributes).

gdal读取NC文件的步骤

1.使用gdalinfo获取NC文件的中包含的子数据信息。

2.根据子数据信息，打开对应的数据集。

gdalinfo获取NC子数据信息

在命令行输入 gdalinfo xxx.nc

以下是用sentinel3数据为示例，如下：

打开对应的数据集

因为 HDF5 、nc是 GDAL 支持的格式，所以GDAL是一种打开HDF5的一种方法

from osgeo import gdal

hdf_ds = gdal.Open("/path/to/hdf/hdf_file.h5", gdal.GA_ReadOnly)

# Replace the <subdataset> with the band you want to read, indexing from 0
band_ds = gdal.Open(hdf_ds.GetSubDatasets()[<subdataset>][0], gdal.GA_ReadOnly)

band_array = band_ds.ReadAsArray()

band_array然后将是一个包含子数据集值的 numpy 数组。

可以使用该方法或从命令行.GetSubDatasets()使用来找出要寻址的子数据集。

以下是具体例子，使用GDAL获取Sentinel-3 数据的经度、纬度、高度

Sentinel-3 数据地理空间信息保存在geo_coordinates.nc，该文件共有三个波段，第一个波段到第三个波段分别是高度、纬度、经度。

from osgeo import gdal, osr
geo_file = r'xxx\\geo_coordinates.nc'
ds = gdal.Open(geo_file, gdal.GA_ReadOnly)
altitude_ds = gdal.Open(ds.GetSubDatasets()[0][0], gdal.GA_ReadOnly)
latitude_ds = gdal.Open(ds.GetSubDatasets()[1][0], gdal.GA_ReadOnly)
longitude_ds = gdal.Open(ds.GetSubDatasets()[2][0], gdal.GA_ReadOnly)

altitude_array = altitude_ds.ReadAsArray()
latitude_array = latitude_ds.ReadAsArray() * 0.000001
longitude_array = longitude_ds.ReadAsArray() * 0.000001

在nc文件中保存的数值一般是整型，为了保证数值的精确性一般要进行缩放。例如，假设纬度是21.111578°，如果用浮点型进行保存，存储空间占比大。所以，欧空局对此的操作是，把经纬度放大1000000倍进行保存，变成21111578。相对的，我们要读取数据则要进行缩小1000000倍。这样做的好处在于存储空间大大减少。

缩放系数0.000001，在gdalinfo 读取信息发现。

或者使用gdal中的 ds.GetMetadata()函数获取信息.

from osgeo import gdal, osr
geo_file = r'xxx\\geo_coordinates.nc'
ds = gdal.Open(geo_file, gdal.GA_ReadOnly)
print(ds.GetMetadata())

输出如下：

Out[1]: 
{'absolute_orbit_number': '35981 ',
 'altitude_long_name': 'DEM corrected altitude',
 'altitude_positive': 'up',
 'altitude_standard_name': 'altitude',
 'altitude_units': 'meters',
 'altitude_valid_max': '9000 ',
 'altitude_valid_min': '-1000 ',
 'altitude__FillValue': '-32768 ',
 'columns_CLASS': 'DIMENSION_SCALE',
 'columns_NAME': 'This is a netCDF dimension but not a netCDF variable.      4865',
 'columns_REFERENCE_LIST': '',
 'comment': ' ',
 'contact': 'eosupport@copernicus.esa.int',
 'creation_time': '2022-01-15T21:01:27Z',
 'history': '  2022-02-15T22:00:27Z: PUGCoreProcessor joborder.21059.xml',
 'institution': 'PS1',
 'lat_long_name': 'DEM corrected latitude',
 'lat_scale_factor': '1e-06 ',
 'lat_standard_name': 'latitude',
 'lat_units': 'degrees_north',
 'lat_valid_max': '90000000 ',
 'lat_valid_min': '-90000000 ',
 'lat__FillValue': '-2147483648 ',
 'lon_long_name': 'DEM corrected longitude',
 'lon_scale_factor': '1e-06 ',
 'lon_standard_name': 'longitude',
 'lon_units': 'degrees_east',
 'lon_valid_max': '180000000 ',
 'lon_valid_min': '-180000000 ',
 'lon__FillValue': '-2147483648 ',
 'netCDF_version': '4.2 of Mar 13 2018 10:14:33 $',
 'processing_baseline': 'SYN_L2_.002.16.00',
 'product_name': 'S3A_SY_2_SYN___XXXXX.SEN3',
 'references': 'S3IPF PDS 006 - i1r14 - Product Data Format Specification - SYNERGY, S3IPF PDS 002 - i1r8 - Product Data Format Specification - Product Structures, S3IPF DPM 004 - i1r11 - Detailed Processing Model - SYN',
 'resolution': '[ 300 300 ]',
 'rows_CLASS': 'DIMENSION_SCALE',
 'rows_NAME': 'This is a netCDF dimension but not a netCDF variable.      4091',
 'rows_REFERENCE_LIST': '',
 'source': 'IPF-SY-2 06.23',
 'start_time': '2022-02-14T02:42:05.084552Z',
 'stop_time': '2022-02-14T02:45:05.049627Z',
 'title': 'SYN L2, high resolution georeferencing data'}