H8/AHI数据预处理-生成带投影信息的tif格式数据

H8数据预处理的步骤为：

读取所有反射率/亮度温度数据：
反射率数据 = 反射率数据0.0001
亮度温度数据 = 亮度温度数据0.01+273.15
读取经度和纬度数据，投影方式为GLL, 即等经纬度投影，
等经纬度投影的经纬度范围为（0：360，90：-90）
其中左上角经纬度：80，60，右下角经纬度：200，-60
像素分辨率的0.02°，空间分辨率为2km
读取其他数据，如果需要（SOA,SZA,SAA,SAZ）
左上角x坐标，水平分辨率，旋转参数，左上角y坐标，旋转参数，竖直分辨率
Geotransform = (80,0.02,0,60,0,0.02)

1.NC文件读取

废话不多说，直接上代码！ filename为NC文件全路径，sds为读取的数据集，需要加引号

def ReadNC(filename, sds):
    if os.path.exists(filename):
        nc_obj = Dataset(filename)
        data = nc_obj.variables[sds][:]

    else:
        print("%s is not exists, Read HDF Data fileld!" % (filename))

    return data

2.2.反射率和亮温定标

我这里比较暴力，直接挨个读取！，还有其他的数据，如4个角度（方位角、天顶角需要可以自行添加）

 # 读取 反射率数据并预处理
    Ref_01 = ReadNC(filename, 'albedo_01')*0.0001
    Ref_02 = ReadNC(filename, 'albedo_02')*0.0001
    Ref_03 = ReadNC(filename, 'albedo_03')*0.0001
    Ref_04 = ReadNC(filename, 'albedo_04')*0.0001
    Ref_05 = ReadNC(filename, 'albedo_05')*0.0001
    Ref_06 = ReadNC(filename, 'albedo_06')*0.0001
    # 读取 亮温数据并预处理
    # 亮温：Tbb*0.01+273.15
    Tbb_07 = ReadNC(filename, 'tbb_07')*0.01+273.15
    Tbb_08 = ReadNC(filename, 'tbb_08')*0.01+273.15
    Tbb_09 = ReadNC(filename, 'tbb_09')*0.01+273.15
    Tbb_10 = ReadNC(filename, 'tbb_10')*0.01+273.15
    Tbb_11 = ReadNC(filename, 'tbb_11')*0.01+273.15
    Tbb_12 = ReadNC(filename, 'tbb_12')*0.01+273.15
    Tbb_13 = ReadNC(filename, 'tbb_13')*0.01+273.15
    Tbb_14 = ReadNC(filename, 'tbb_14')*0.01+273.15
    Tbb_15 = ReadNC(filename, 'tbb_15')*0.01+273.15
    Tbb_16 = ReadNC(filename, 'tbb_16')*0.01+273.15

3.生成tif格式数据

生成tif前，需要先获取投影信息（6个参数：左上角经纬度、经纬度分辨率，起始位置），定义坐标系统（WGS84）

def GetTiff(TifName, array):
    cols = array.shape[1]  # 矩阵列数
    rows = array.shape[0]  # 矩阵行数
    bands = array.shape[2]
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(TifName, cols, rows, bands, gdal.GDT_Float32)
    # 括号中两个0表示起始像元的行列号从(0,0)开始，(80,60)为起始经纬度，0.02为像素经纬度分辨率（注意：H8数据为2km分辨率，对应像素为0.02°；简单说下来源：地球赤道一圈大约40000Km，对应360°，因此，1°为111Km左右，为方便计算，一般习惯上100Km对应1°）
    outRaster.SetGeoTransform([80, 0.02, 0, 60, 0, 0.02])
    # 获取数据集第一个波段，从1开始
    for i in range(bands):
        outband = outRaster.GetRasterBand(i + 1)
        outband.WriteArray(array[:, :, i])
    outRasterSRS = osr.SpatialReference()
    # 4326表示WGS84坐标
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()