PostProcessing

  1# coding: utf-8
  2# Author: Jianbo Qi
  3# Date: 2020/7/6
  4import os
  5import shutil
  6from osgeo import gdal
  7import math
  8import sys
  9import numpy as np
 10
 11
 12class PostProcessing:
 13    @staticmethod
 14    def readIrr(irrFilePath):
 15        # read irrfile
 16        f = open(irrFilePath)
 17        index = 0
 18        sunirr = []
 19        skyirr = None
 20        for line in f:
 21            if index == 1:
 22                arr = line.replace("\n", "").strip().split(" ")
 23                sunirr = list(map(lambda x: float(x), arr[1:]))
 24            if index == 2:
 25                arr = line.replace("\n", "").strip().split(" ")
 26                skyirr = list(map(lambda x: float(x), arr[1:]))
 27            index = index + 1
 28
 29        if skyirr is not None:
 30            sunirr = list(map(lambda x, y: x + y, sunirr, skyirr))
 31        f.close()
 32        return sunirr
 33
 34    @staticmethod
 35    def brf_single_img_processing(radianceFilePath, sunirr, outFilePath):
 36        dataset = gdal.Open(radianceFilePath)
 37        band_num = dataset.RasterCount
 38        XSize, YSize = dataset.RasterXSize, dataset.RasterYSize
 39
 40        format = "ENVI"
 41        driver = gdal.GetDriverByName(format)
 42        dst_ds = driver.Create(outFilePath, XSize, YSize, band_num, gdal.GDT_Float32)
 43        meanBRFs = []
 44        for i in range(0, band_num):
 45            band = dataset.GetRasterBand(i + 1)
 46            dataarr = band.ReadAsArray(0, 0, band.XSize, band.YSize)
 47            dataarr = dataarr / sunirr[i] * math.pi
 48            dst_ds.GetRasterBand(i + 1).WriteArray(dataarr)
 49            tmp = dataarr[dataarr >= 0]
 50            meanBRF = tmp.mean()
 51            meanBRFs.append(meanBRF)
 52        dst_ds = None
 53        # processing header
 54        shutil.copy(radianceFilePath + ".hdr", outFilePath + ".hdr")
 55
 56        return meanBRFs
 57
 58    @staticmethod
 59    def fourcomponents2fvc(fourcomponents_file, out_file):
 60        dataset = gdal.Open(fourcomponents_file)
 61        band_num = dataset.RasterCount
 62        XSize = dataset.RasterXSize
 63        YSize = dataset.RasterYSize
 64        if band_num != 5:
 65            print("Fourcomponent file is not valid, please check the band number.")
 66            return
 67        ill_obj_data = dataset.GetRasterBand(3).ReadAsArray(0, 0, XSize, YSize)
 68        shade_obj_data = dataset.GetRasterBand(5).ReadAsArray(0, 0, XSize, YSize)
 69        fvc = ill_obj_data + shade_obj_data
 70        format = "ENVI"
 71        driver = gdal.GetDriverByName(format)
 72        dst_ds = driver.Create(out_file, XSize, YSize, 1, gdal.GDT_Float32)
 73        dst_ds.GetRasterBand(1).WriteArray(fvc)
 74        dst_ds = None
 75
 76
 77    @staticmethod
 78    def radiance2brf(sim_dir, input_radiace_file="", output_brf_file=""):
 79        if input_radiace_file != "" and output_brf_file != "":
 80            sunirr = PostProcessing.readIrr(os.path.join(sim_dir, "Results", "Irradiance.txt"))
 81            if input_radiace_file == "" or output_brf_file == "":
 82                return
 83            else:
 84                if os.path.exists(input_radiace_file):
 85                    meanBRFs = PostProcessing.brf_single_img_processing(input_radiace_file, sunirr, output_brf_file)
 86        elif input_radiace_file == "" and output_brf_file == "":
 87            spectral_info_path = os.path.join(sim_dir, "Results", "spectral.txt")
 88            if os.path.exists(spectral_info_path):
 89                fs = open(spectral_info_path, "r")
 90                f = open(os.path.join(sim_dir, "Results", "spectral_BRF.txt"), 'w')
 91                sunirr = PostProcessing.readIrr(os.path.join(sim_dir, "Results", "Irradiance.txt"))
 92                for line in fs:
 93                    radiance_path = os.path.join(sim_dir, "Results", line)
 94                    if os.path.exists(radiance_path):
 95                        print("INFO: Processing Image: " + line)
 96                        output_file = radiance_path + "_BRF"
 97                        meanBRFs = PostProcessing.brf_single_img_processing(radiance_path, sunirr,
 98                                                                            output_file)  # mean BRF for each band
 99                        f.write(line + " ")
100                        for i in range(0, len(meanBRFs)):
101                            f.write("%.5f " % meanBRFs[i])
102                        f.write("\n")
103                f.close()
104        else:
105            print("Error: input or output file is not specified.")
106            sys.exit()
107
108
109def planck_invert(radiance, wavelength):
110    if radiance <= 0:
111        return 0
112    kb = 1.38064852e-23  # Boltzmann constant
113    hp = 6.626070040e-34  # Planck constant
114    c = 299792458
115    wavelength *= math.pow(10, -9)
116    radiance *= math.pow(10, 9)
117    inside = 1+2*hp*c*c/(math.pow(wavelength,5)*radiance)
118    down = wavelength*kb*math.log(inside)
119    up = hp*c
120    return up/down
121
122
123class RasterHelper:
124    @staticmethod
125    def saveToHdr_no_transform(npArray, dstFilePath, wlist, output_format):
126        dshape = npArray.shape
127        if len(dshape) == 3:
128            bandnum = dshape[2]
129        else:
130            bandnum = 1
131        # 从hdrHeaderPath中提取投影信息
132        if output_format == "ENVI":
133            format = "ENVI"
134        else:
135            format = "GTiff"
136            dstFilePath += ".tif"
137        driver = gdal.GetDriverByName(format)
138        dst_ds = driver.Create(dstFilePath, dshape[1], dshape[0], bandnum, gdal.GDT_Float32)
139        #     npArray = linear_stretch_3d(npArray)
140        if bandnum > 1:
141            for i in range(1, bandnum + 1):
142                dst_ds.GetRasterBand(i).WriteArray(npArray[:, :, i - 1])
143        else:
144            dst_ds.GetRasterBand(1).WriteArray(npArray)
145        dst_ds = None
146
147        if output_format == "ENVI" and len(wlist) > 0:
148            # wirte wavelength
149            f = open(dstFilePath + ".hdr", 'r')
150            text = f.read()
151            f.close()
152            wstr = "\nwavelength = {"
153            for i in range(0, len(wlist)):
154                wstr += wlist[i].split(":")[0] + ","
155            wstr = wstr[0:len(wstr) - 1] + "}"
156            f = open(dstFilePath + ".hdr", 'w')
157            text = text + wstr
158            f.write(text)
159            f.close()
160
161    @staticmethod
162    def convertRadiance2BTimage(data, wlist):
163        wavelengths = []
164        for i in range(0, len(wlist)):
165            wavelengths.append(float(wlist[i].split(":")[0]))
166
167        dshape = data.shape
168
169        re = np.zeros(dshape)
170        if len(dshape) == 3:  # multiple bands
171            rows, cols, depth = dshape
172            for i in range(0, dshape[2]):
173                banddata = data[:, :, i]
174                for r in range(rows):
175                    for c in range(cols):
176                        re[r][c][i] = planck_invert(banddata[r][c], wavelengths[i])
177        else:
178            rows, cols = dshape
179            for r in range(rows):
180                for c in range(cols):
181                    re[r][c] = planck_invert(data[r][c], wavelengths[0])
182        return re