corning410_radiance_calibration/1.1radiance_calibration_300tc.py

'''
1、此版本所需文件：dn影像、dark影像和asd测定的积分球辐亮度曲线
2、处理过程：1）扣除暗电流；2）通过辐亮度曲线生成定标文件
'''
import numpy as np
import pandas as pd
import sys, os, traceback, re
import shutil
from radiance_calibration_ui import Ui_MainWindow
from library.multithread import Worker
import sys, traceback
from osgeo import gdal

from PyQt5.QtWidgets import QMainWindow, QFileDialog, QApplication
from PyQt5.QtCore import Qt


class EnterWindow(QMainWindow, Ui_MainWindow):
    def __init__(self, parent=None):
        super(EnterWindow, self).__init__(parent)
        self.setupUi(self)

        self.radiance_calibration_object = RadianceCalibration()

        # self.setWindowState(Qt.WindowMaximized)  # 初始化时就最大化窗口

        self.rad_bt.clicked.connect(self.select_rad)
        self.dn_bt.clicked.connect(self.select_dn)
        self.dark_bt.clicked.connect(self.select_dark)
        self.out_file_bt.clicked.connect(self.select_out_file)

        self.operate_bt.clicked.connect(self.operate)

        self.tmp_rad_file_path = None  # 用于保存当前路径，下次打开就默认在此路径下
        self.tmp_dn_file_path = None  # 用于保存当前路径，下次打开就默认在此路径下
        self.tmp_dark_file_path = None  # 用于保存当前路径，下次打开就默认在此路径下
        self.tmp_out_file_path = None  # 用于保存当前路径，下次打开就默认在此路径下

    def select_rad(self):
        if self.tmp_rad_file_path == None:
            rad_file_path = QFileDialog.getOpenFileName(self, '选择asd辐亮度文件', os.path.dirname(__file__))[0]

            self.tmp_rad_file_path = os.path.dirname(rad_file_path)
        elif self.tmp_rad_file_path is not None:
            rad_file_path = QFileDialog.getOpenFileName(self, '选择asd辐亮度文件', self.tmp_rad_file_path)[0]

            self.tmp_rad_file_path = os.path.dirname(rad_file_path)

        self.radiance_calibration_object.asd_radiance_file_path = rad_file_path
        self.rad_le.setText(rad_file_path)

    def select_dn(self):
        if self.tmp_dn_file_path == None:
            dn_file_path = QFileDialog.getOpenFileName(self, '选择dn影像', os.path.dirname(__file__))[0]

            self.tmp_dn_file_path = os.path.dirname(dn_file_path)
        elif self.tmp_dn_file_path is not None:
            dn_file_path = QFileDialog.getOpenFileName(self, '选择dn影像', self.tmp_dn_file_path)[0]

            self.tmp_dn_file_path = os.path.dirname(dn_file_path)

        self.radiance_calibration_object.dn_file_path = dn_file_path
        self.dn_le.setText(dn_file_path)

    def select_dark(self):
        if self.tmp_dark_file_path == None:
            dark_file_path = QFileDialog.getOpenFileName(self, '选择dark影像', os.path.dirname(__file__))[0]

            self.tmp_dark_file_path = os.path.dirname(dark_file_path)
        elif self.tmp_dark_file_path is not None:
            dark_file_path = QFileDialog.getOpenFileName(self, '选择dark影像', self.tmp_dark_file_path)[0]

            self.tmp_dark_file_path = os.path.dirname(dark_file_path)

        self.radiance_calibration_object.dark_file_path = dark_file_path
        self.dark_le.setText(dark_file_path)

    def select_out_file(self):
        if self.tmp_out_file_path == None:
            out_file_path = QFileDialog.getExistingDirectory(self, '选择输出路径', os.path.dirname(__file__))

            self.tmp_out_file_path = out_file_path
        elif self.tmp_out_file_path is not None:
            out_file_path = QFileDialog.getExistingDirectory(self, '选择输出路径', self.tmp_out_file_path)

            self.tmp_out_file_path = out_file_path

        self.radiance_calibration_object.out_file_path = out_file_path
        self.out_file_le.setText(out_file_path)

    def operate(self):
        self.radiance_calibration_object.operate()


class RadianceCalibration():
    def __init__(self):
        # asd_radiance_file = os.path.dirname(__file__) + '\\corning_radance.txt'
        # img_file_path = os.path.dirname(__file__) + '\\jfq_dn'
        # img_dark_file_path = os.path.dirname(__file__) + '\\jfq_dn_dark'
        # out_file = os.path.dirname(__file__) + '\\jfq_dn_gain'

        self.asd_radiance_file_path = None
        self.dn_file_path = None
        self.dark_file_path = None

        self.img_data = None
        self.img_dark_data = None

        self.out_file_path = None

    def get_envi_header_dict(self, hdr_file):
        with open(hdr_file, encoding='utf-8') as file_obj:
            contents = file_obj.read()

            # Get all "key = {val}" type matches
            regex = re.compile(r'^(.+?)\s*=\s*({\s*.*?\n*.*?})$', re.M | re.I)
            matches = regex.findall(contents)

            # Remove them from the header
            subhdr = regex.sub('', contents)

            # Get all "key = val" type matches
            regex = re.compile(r'^(.+?)\s*=\s*(.*?)$', re.M | re.I)
            matches.extend(regex.findall(subhdr))

            hdr_fields = dict(matches)
            keys_delete = []
            for key in hdr_fields.keys():
                if (
                        key == "interleave" or key == "byte order" or key == "data type" or key == "samples" or key == "bands" or key == "lines" or key == "framerate"):
                    keys_delete.append(key)

            for x in range(len(keys_delete)):
                hdr_fields.pop(keys_delete[x])

            return hdr_fields

    def get_hdr_filename(self, file_path):
        return os.path.splitext(file_path)[0] + ".hdr"

    def get_hdr_filename_with_bip(self, file_path):
        return os.path.splitext(file_path)[0] + ".bip.hdr"

    def write_fields_to_hdrfile(self, fields, hdr_file):
        with open(hdr_file, "a", encoding='utf-8') as f:
            for key in fields.keys():
                f.write(key + " = " + fields[key] + "\n")

    def process_hdr(self, hdr_file_path, envi_header_dict):
        self.write_fields_to_hdrfile(envi_header_dict, self.get_hdr_filename(hdr_file_path))
        shutil.copyfile(self.get_hdr_filename(hdr_file_path), self.get_hdr_filename_with_bip(hdr_file_path))
        os.remove(self.get_hdr_filename(hdr_file_path))

    # 读图像文件
    def read_img(self, 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 AttributeError:
            print('gdal打开影像：没有文件！')
        except Exception:
            traceback.print_exc()
            pass

    # 写文件，以写成tif为例
    def write_img(self, 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_Float64,
                               options=["INTERLEAVE=BIP"])
        for i in range(band):
            dst_ds.GetRasterBand(i + 1).WriteArray(data[i, :, :])  # gdal的band从1开始，所以dst_ds.GetRasterBand(i+1)
        dst_ds = None

    def operate(self):
        try:
            # 读取asd辐亮度数据
            data = pd.read_csv(self.asd_radiance_file_path, sep='\t', dtype=np.float64, header=None)
            self.asd_radiance = np.array(data)

            # 读取影像
            img_proj, img_geotrans, self.img_data = self.read_img(self.dn_file_path)
            dn_hdr = self.get_envi_header_dict(self.get_hdr_filename(self.dn_file_path))
            img_dark_proj, img_dark_geotrans, self.img_dark_data = self.read_img(self.dark_file_path)
            dark_hdr = self.get_envi_header_dict(self.get_hdr_filename(self.dark_file_path))

            # 将影像所有行平均，得到一行（帧）影像
            img_data_ave = np.mean(self.img_data, axis=1)
            img_dark_data_ave = np.mean(self.img_dark_data, axis=1)

            # 去除暗电流
            img_data_ave_rmdark = img_data_ave - img_dark_data_ave
            img_dark_data_ave = img_dark_data_ave + 50

            gain = np.empty((img_data_ave_rmdark.shape[0], 1, img_data_ave_rmdark.shape[1]))
            offset = np.empty((img_data_ave_rmdark.shape[0], 1, img_data_ave_rmdark.shape[1]))
            for i in range(gain.shape[2]):
                gain[:, 0, i] = self.asd_radiance[:, 1] / img_data_ave_rmdark[:, i]
                offset[:, 0, i] = img_dark_data_ave[:, i]

            gain_name = os.path.join(self.out_file_path, "gain.bip")
            offset_name = os.path.join(self.out_file_path, "offset.bip")
            self.write_img(gain_name, gain)
            self.process_hdr(gain_name, dn_hdr)
            self.write_img(offset_name, offset)
            self.process_hdr(offset_name, dark_hdr)
        except:
            traceback.print_exc()

    def operate_without_dark(self):
        try:
            # 读取asd辐亮度数据
            data = pd.read_csv(self.asd_radiance_file_path, sep='\t', dtype=np.float64, header=None)
            self.asd_radiance = np.array(data)

            # 读取影像
            img_proj, img_geotrans, self.img_data = self.read_img(self.dn_file_path)

            # 将影像所有行平均，得到一行（帧）影像
            img_data_ave = np.mean(self.img_data, axis=1)

            gain = np.empty((img_data_ave.shape[0], 1, img_data_ave.shape[1]))
            for i in range(gain.shape[2]):
                gain[:, 0, i] = self.asd_radiance[:, 1] / img_data_ave[:, i]

            self.write_img(self.out_file_path, gain)
        except:
            traceback.print_exc()


if __name__ == '__main__':
    app = QApplication(sys.argv)

    # 实例化主窗口
    enter_window_instance = EnterWindow()
    enter_window_instance.show()

    sys.exit(app.exec_())