c. HDF文件读取,MODIS波段简介
卫星云图数据存储的格式相当多样,其中不少格式需要较复杂的脚本才能读取(比如Himawari-8/9最广泛使用的DAT格式)。相对容易读取的数据存储方式包括HDF4,HDF5,NetCDF等。MODIS的L1B数据正是用HDF4格式存储的。HDF4有多种读取方式,这里我们选择使用第三方库pyhdf。HDF5和NetCDF我使用的分别是h5py和netCDF4,在此不作详细介绍。
pyhdf的介绍:
https://fhs.github.io/pyhdf/modules/SD.html,下图列举了部分会用到的功能
from pyhdf.SD import SD, SDC
import numpy as np
首先导入pyhdf与numpy,此处pyhdf只导入了其中的一部分模块,numpy则是导入时命名为了np,后续反复调用的时候少打几个字母。
file = SD(r'E:\Data\MOD021KM.A2015296.0515.061.2017322234443.hdf', SDC.READ)
在只读模式(SDC.READ)下打开文件,文件路径请自行替换,这里是用了Patricia10.23 0515z的Terra,此处路径前的r是为了直接拷贝过来的绝对路径能够正常读取。
我们先使用datasets()看看文件中放了哪些数据,这一步相当于列出数据集的“目录”,还没有真正将数据取出来:
data = file.datasets()
print(data)
print(data.keys())
直接print(data)会返回每个数据集的名称与一些附带信息,加上keys()能仅返回数据集名称,以下是print(data.keys())的结果:
dict_keys(['Latitude', 'Longitude', 'EV_1KM_RefSB', 'EV_1KM_RefSB_Uncert_Indexes', 'EV_1KM_Emissive', 'EV_1KM_Emissive_Uncert_Indexes', 'EV_250_Aggr1km_RefSB', 'EV_250_Aggr1km_RefSB_Uncert_Indexes', 'EV_250_Aggr1km_RefSB_Samples_Used', 'EV_500_Aggr1km_RefSB', 'EV_500_Aggr1km_RefSB_Uncert_Indexes', 'EV_500_Aggr1km_RefSB_Samples_Used', 'Height', 'SensorZenith', 'SensorAzimuth', 'Range', 'SolarZenith', 'SolarAzimuth', 'gflags', 'EV_Band26', 'EV_Band26_Uncert_Indexes', 'Band_250M', 'Band_500M', 'Band_1KM_RefSB', 'Band_1KM_Emissive', 'Noise in Thermal Detectors', 'Change in relative responses of thermal detectors', 'DC Restore Change for Thermal Bands', 'DC Restore Change for Reflective 250m Bands', 'DC Restore Change for Reflective 500m Bands', 'DC Restore Change for Reflective 1km Bands'])
其中'Latitude','Longitude'是经纬度数据,'EV_1KM_RefSB'是反射波段数据,这里暂时用不到,'EV_1KM_Emissive'是热辐射波段数据,包含了我们所需要的各个红外波段。
使用select将'Latitude','Longitude','EV_1KM_Emissive'读取出来,并用get()读取数组,用np.shape()获取每个数组的形状(不过其实print(data)也会将数组形状作附带信息返回):
lon = file.select('Longitude').get()
lat = file.select('Latitude').get()
Emissive = file.select('EV_1KM_Emissive').get()
print(lon, lat, Emissive)
print(np.shape(lon), np.shape(lat), np.shape(Emissive))
可以看到热辐射波段数组的形状是(16, 2030, 1354),代表有16个不同波段堆叠在一起,每个波段的图像有2030×1354个像素,而数组的值都是0~65535的整数,显然不是我们最终需要的亮温,需经过辐射定标与亮温反演(见2.d);经纬度数组的值是正确的,但数组形状是(406, 271),像素比波段数据少,这是因为文件中存储的经纬度降过分辨率,需插值后才能使用(见2.e)。
'EV_1KM_Emissive'这个数据集并非只包含上面读取的数组,HDF文件可以为数据集附上attributes,其中往往包含辐射定标的关键信息,查看方式如下:
print(file.select('EV_1KM_Emissive').attributes())
print(file.select('EV_1KM_Emissive').attributes().keys())
以下是加上keys(),只获取attributes名称的结果:
dict_keys(['long_name', 'units', 'valid_range', '_FillValue', 'band_names', 'radiance_scales', 'radiance_offsets', 'radiance_units'])
其中'long_name'是数据集全称, 'band_names'是波段号,这两个attribute可以在看官方文献看晕掉的时候让你知道读的是什么数据;'radiance_scales','radiance_offsets'是辐射定标的两个系数(见2.d),'radiance_units'是辐射定标后得到的辐亮度的单位,有时候不同卫星数据会采用不同的辐亮度单位,反演的亮温如果偏差很大则需要检查该单位并调整黑体辐射公式的形式。
获取attributes中的数组同样是使用get()即可:
radiance_scales = file.select('EV_1KM_Emissive').attributes().get('radiance_scales')
radiance_offsets = file.select('EV_1KM_Emissive').attributes().get('radiance_offsets')
MODIS波段介绍:
https://modis.gsfc.nasa.gov/about/specifications.php先前查看热辐射波段attributes的时候16个波段号分别为20,21,22,23,24,25,27,28,29,30,31,32,33,34,35,36。细心的读者可能已经发现了,Band26不在其中。(不细心的作者因为这个在伪VIS教程篇里面把中波红外要用的Band20错写成Band21了,前两天刚勘误)Band26是短波红外反射波段,在2.d中读取指定波段数据的时候需要注意到这一点,不然可能会读错波段。
16个热辐射波段中,Band31用于渲染我们最熟知的IR云图,Band27则是标准水汽波段,Band28是响应层面偏低一些的水汽波段,反映中低层大气的暖心和水汽吸收情况。伪VIS所用到的波段则是Band31,Band32与Band20,波段选择的具体原因可以参考《伪VIS算法原理与教程 - 理论篇》。值得一提的是,Band20与Band21,22的波长非常相近,选择Band20的主要原因是Band20接收的波长宽一些,更容易接收到较冷云顶上微弱的辐射,云顶噪点相对而言少一些。