python FY-3D MWRI亮温每轨影像hdf转tif-526互联

FY-3D MWRI亮温数据包含单独的经纬度数据集和亮温数据集（维度大小一样），需要把HDF数据转成0.1度的WGS84 tif数据。

方法：gdal.Translate + gdal.Warp

1. 先使用gdal.Translate 函数构建VRT文件

gdal.Translate(vrtDir, subDataset, format='vrt')

vrtDir为VRT文件的输出路径，subDataset为需要转换的数据集。
接着需要在VRT文件中写入GEOLOCATION元数据，X_DATASET和Y_DATASET分别指定了经度和纬度数据的存储位置。

        lines = []
        with open(vrtDir, 'r') as f:
            for line in f:
                lines.append(line)
        lines.insert(1, '<Metadata domain="GEOLOCATION">\n')
        lines.insert(2, ' <MDI key="LINE_OFFSET">1</MDI>\n')
        lines.insert(3, ' <MDI key="LINE_STEP">1</MDI>\n')
        lines.insert(4, ' <MDI key="PIXEL_OFFSET">1</MDI>\n')
        lines.insert(5, ' <MDI key="PIXEL_STEP">1</MDI>\n')
        lines.insert(6,
                     ' <MDI key="SRS">GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],'
                     'TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],'
                     'UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]</MDI>\n')
        lines.insert(7, ' <MDI key="X_BAND">1</MDI>')
        lines.insert(8, ' <MDI key="X_DATASET">HDF5:"{}"://GEOLOCATION/Longitude</MDI>\n'.format(geo_file))
        lines.insert(9, ' <MDI key="Y_BAND">1</MDI>\n')
        lines.insert(10, ' <MDI key="Y_DATASET">HDF5:"{}"://GEOLOCATION/Latitude</MDI>\n'.format(geo_file))
        lines.insert(11, '</Metadata>\n')
        if 'GEO' in vrtDir and 'vrt' in vrtDir:
            with open(vrtDir, 'w') as f:
                for line in lines:
                    f.writelines(line)

但需要注意的是，FY-3D MWRI文件中经纬度数据存在填充值65535，这样会导致后面gdal.Warp转换出错，所以需要去掉经纬度数据集中的填充值后，再将其写入VRT文件。

# 删除65535所在的行
        data = h5py.File(inputfile)
        Lon_set = data['Geolocation/Longitude']
        lon = Lon_set[:, :]
        lon = np.where(lon == 65535, np.nan, lon)
        lon_xin = lon[~np.isnan(lon).any(axis=1)]

        Lat_set = data['Geolocation/Latitude']
        lat = Lat_set[:, :]
        lat = np.where(lat == 65535, np.nan, lat)
        lat_xin = lat[~np.isnan(lat).any(axis=1)]

# 将修改后的经纬度矩阵写入新HDF

with h5py.File(geofile, 'w') as hf:
hf.create_dataset("lon-dataset", data=lon_xin)
hf.create_dataset("lat-dataset", data=lat_xin)

然后将包含修改后经纬度的写入VRT文件，即修改X_DATASET和Y_DATASET

lines.insert(8, ' <MDI key="X_DATASET">HDF5:"{}"://lon-dataset</MDI>\n'.format(geo_file))
lines.insert(10, ' <MDI key="Y_DATASET">HDF5:"{}"://lat-dataset</MDI>\n'.format(geo_file))

2. 使用gdal.Warp 函数进行几何校正

VrtData = gdal.Warp(tiff_file, vrtDir,
                            dstNodata=-32768,
                            format='GTiff', geoloc=True, outputType=dst_dtype,
                            resampleAlg=resample_alg, xRes=Res_x, yRes=Res_y)

参考：https://blog.csdn.net/tk20190411/article/details/108547883