'''
https://blog.csdn.net/vonuo/article/details/74783291
本方法是使用asd的方法来进行corning的辐亮度定标（resonon的方法感觉和此方法差不多）
'''
from osgeo import gdal
import os, math
import sys
import easygui
import numpy as np
import pandas as pd
import xlwt

# 读写影像类
class Grid(object):

    #读图像文件
    @classmethod
    def read_img(cls, filename, xoff=0, yoff=0, im_width=None, im_height=None):
        try:
            dataset = gdal.Open(filename)  # 打开文件
            if im_width == None:
                im_width = dataset.RasterXSize  # 栅格矩阵的列数
            if im_height == None:
                im_height = dataset.RasterYSize  # 栅格矩阵的行数
            num_bands = dataset.RasterCount  # 栅格矩阵的波段数
            im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵
            im_proj = dataset.GetProjection()  # 地图投影信息
            im_data = dataset.ReadAsArray(xoff, yoff, im_width, im_height)  # 将数据写成数组，对应栅格矩阵
            del dataset
            return im_proj, im_geotrans, im_data
        except:
            sys.exit()

    #写文件，以写成tif为例
    @classmethod
    def write_img(cls, dst_filename, data):
        format = "ENVI"
        driver = gdal.GetDriverByName(format)
        RasterXSize = data.shape[2]  # 遥感影像的sample（列数）
        RasterYSize = data.shape[1]  # 遥感影像的line（行数）
        band = data.shape[0]
        # driver.Create()函数中RasterXSize代表影像的sample（列数），RasterYSize代表影像的line（行数）
        dst_ds = driver.Create(dst_filename, RasterXSize, RasterYSize, band, gdal.GDT_Float32)
        for i in range(band):
            dst_ds.GetRasterBand(i + 1).WriteArray(data[i, :, :])  # gdal的band从1开始，所以dst_ds.GetRasterBand(i+1)
        dst_ds = None

# 是否转辐亮度：0→不转，1→转
rad_switch = 0

img = r'D:\py_program\corning410\record_system_v24\baiban_record'
img_baiban = r'D:\py_program\corning410\record_system_v24\baiban'
img_gain = r'D:\py_program\corning410\corning410_radiance_calibration\jfq_dn_gain'
img_offset = r'D:\py_program\corning410\record_system_v24\dark'


dirpath = os.path.splitext(img)

# 读取影像
proj, geotrans, data = Grid.read_img(img)
proj_baiban, geotrans_baiban, data_baiban = Grid.read_img(img_baiban)
proj_gain, geotrans_gain, data_gain = Grid.read_img(img_gain)
proj_offset, geotrans_offset, data_offset = Grid.read_img(img_offset)

data_baiban = np.mean(data_baiban, axis=1)
data_offset = np.mean(data_offset, axis=1)

if rad_switch == 1:
    # 计算辐射定标参数
    cal_it = 6059
    target_it = 200004
    gain_scale = cal_it / target_it
    data_gain_adjust = data_gain * gain_scale

    # 影像和白板：1、扣除暗电流；2、转换成辐亮度；
    data_baiban = data_baiban - data_offset  # 白板扣除暗电流
    data_baiban_rad = data_baiban * data_gain_adjust[:, 0, :]  # 白板转辐亮度

    data_rad = np.empty(data.shape)
    for i in range(data.shape[1]):
        data_rad[:, i, :] = (data[:, i, :] - data_offset) * data_gain_adjust[:, 0, :]  # 转辐亮度
    Grid.write_img(dirpath[0] + '_rad', data_rad)

    # 转换成反射率
    data_rad_ref = np.empty(data.shape)
    for i in range(data.shape[1]):
        data_rad_ref[:, i, :] = data_rad[:, i, :] / data_baiban_rad  # 转反射率

    Grid.write_img(dirpath[0] + '_rad_ref', data_rad_ref)
elif rad_switch == 0:
    # 影像和白板扣除暗电流
    data_baiban = data_baiban - data_offset
    data_rmdark = np.empty(data.shape)
    for i in range(data.shape[1]):
        data_rmdark[:, i, :] = (data[:, i, :] - data_offset)

    # 转反射率
    data_ref = np.empty(data.shape)
    for i in range(data.shape[1]):
        data_ref[:, i, :] = data_rmdark[:, i, :] / data_baiban  # 转反射率
    Grid.write_img(dirpath[0] + '_ref', data_ref)