python读取Himawari8数据

影像处理过程中将500m和1km分辨率的都采样到了2km

def read_hsd(inputfile):
    resolution = int(inputfile[-12:-10])
    if resolution == 10:
        rows = 1100
        cols = 11000
    elif resolution == 20:
        rows = 550
        cols = 5500
    else:
        rows = 2200
        cols = 22000
    # 打开文件
    f = open(inputfile, 'rb')
    # 获取Observation start time
    f.seek(46)
    imgtime = struct.unpack('d', f.read(8))[0]
    # 获取Total header length
    f.seek(70)
    header_length = struct.unpack('i', f.read(4))[0]
    # 获取影像
    formation = [('header', 'S1', header_length), ('pixel', 'i2', rows*cols)]
    imgdata = np.fromfile(inputfile, dtype=formation)['pixel'].reshape(rows, cols)
    if resolution != 20:
        # 重采样至550行，5500列
        imgdata = imgdata.reshape(550, int(20/resolution), 5500, int(20/resolution)).mean(axis=(1, 3))
    # 获取Sun's position
    f.seek(510)
    sun_pos1 = struct.unpack('d', f.read(8))[0]
    # print(sun_pos1)
    f.seek(510+8)
    sun_pos2 = struct.unpack('d', f.read(8))[0]
    # print(sun_pos2)
    f.seek(510+8+8)
    sun_pos3 = struct.unpack('d', f.read(8))[0]
    # print(sun_pos3)
    # 获取Band number
    f.seek(601)
    band_num = int.from_bytes(f.read(2), byteorder='little', signed=False)
    # 获取Gain
    f.seek(617)
    Gain = struct.unpack('d', f.read(8))[0]
    # print(Gain)
    # 获取Offset
    f.seek(625)
    Offset = struct.unpack('d', f.read(8))[0]
    # print(Offset)
    # 计算radiance
    radiance = imgdata*Gain+Offset
    del imgdata
    # print(radiance[0][1580])
    # 对前6波段定标成反射率
    if band_num <= 6:
        # 获取radiance to albedo
        f.seek(633)
        cc = struct.unpack('d', f.read(8))[0]
        f.close
        # 计算反射率
        albedo = radiance * cc
        outdata = albedo
        del albedo, radiance
    # 后面的波段定标成计算亮温
    else:
        # 获取Central wave length
        f.seek(603)
        wv = struct.unpack('d', f.read(8))[0]
        # 获取radiance to brightness temperature(c0)
        f.seek(633)
        c0 = struct.unpack('d', f.read(8))[0]
        # 获取radiance to brightness temperature(c1)
        f.seek(641)
        c1 = struct.unpack('d', f.read(8))[0]
        # 获取radiance to brightness temperature(c2)
        f.seek(649)
        c2 = struct.unpack('d', f.read(8))[0]
        # 获取Speed of light(c)
        f.seek(681)
        c = struct.unpack('d', f.read(8))[0]
        # 获取Planck constant(h)
        f.seek(689)
        h = struct.unpack('d', f.read(8))[0]
        # 获取Boltzmann constant(k)
        f.seek(697)
        k = struct.unpack('d', f.read(8))[0]
        f.close

        # 计算亮温
        wv = wv * 1e-6
        rad = radiance * 1e6
        del radiance
        Te = h*c/k/wv/(np.log(2*h*c*c/((wv**5)*rad)+1))
        del rad
        BT = c0 + c1 * Te + c2 * Te * Te
        del Te
        # 返回值
        outdata = BT
        del BT
    # 返回：albedo / BT, 时间, 太阳坐标
    sunpos = [sun_pos1, sun_pos2, sun_pos3]

    out = {'pixels': list(outdata.flatten()), 'time': imgtime, 'sun_pos': sunpos}
    del outdata, sunpos
    return out

# 获取影像数据(影像、时间、太阳坐标)
def H8_Getdata(inputfolder):
    count = 0
    for band in range(1, 17):
        for seg in range(1, 11):
            seg_file = glob.glob(inputfolder + '/HS_H08*_B' + ('0'+str(band))[-2:] + \
                                   '_FLDK_*_S' + ('0'+str(seg))[-2:] + '10.DAT')
            if seg_file != '':
                count += 1
    if count == 160:
        R20_data = np.empty([5500, 5500, 16], dtype=np.float)
        time_data = np.empty([5500, 5500, 4], dtype=np.float)
        for band in range(16):
            # print(band + 1)
            band_data = []
            obstime = []
            sunpos1 = []
            sunpos2 = []
            sunpos3 = []
            for seg in range(1, 11):
                inputfile = glob.glob(inputfolder + '/HS_H08*_B' + ('0'+str(band+1))[-2:] + \
                                   '_FLDK_R*_S' + ('0'+str(seg))[-2:] + '10.DAT')
                # print(inputfile)
                segment = read_hsd(inputfile[0])
                imgdata = segment['pixels']
                band_data = band_data + imgdata
                if band == 15:
                    ones = np.ones([1, 550*5500], dtype=np.float)
                    imgtime = segment['time']
                    sun_pos = segment['sun_pos']
                    obstime = obstime + list(imgtime * ones.flatten())
                    sunpos1 = sunpos1 + list(sun_pos[0] * ones.flatten())
                    sunpos2 = sunpos2 + list(sun_pos[1] * ones.flatten())
                    sunpos3 = sunpos3 + list(sun_pos[2] * ones.flatten())
                    del ones
                del segment

            # 写入影像数据
            R20_data[:, :, band] = np.array(band_data).reshape(5500, 5500)
            # 写入时间和太阳坐标
            if obstime != []:
                time_data[:, :, 0] = np.array(obstime).reshape(5500, 5500)
                time_data[:, :, 1] = np.array(sunpos1).reshape(5500, 5500)
                time_data[:, :, 2] = np.array(sunpos2).reshape(5500, 5500)
                time_data[:, :, 3] = np.array(sunpos3).reshape(5500, 5500)
            # 删除变量
            del band_data, obstime, sunpos1, sunpos2, sunpos3

        # 返回值
        return {'imgdata': R20_data, 'time_data': time_data}

    else:
        return 0