# -*- coding: UTF-8 -*-
"""
@Project:microproduct 
@File:SalinityALg.py
@Function:土壤粗糙度计算方法
@Contact:
@Author:SHJ
@Date:2021/8/23
@Version:1.0.0
"""
import logging
from tool.algorithm.image.ImageHandle import ImageHandler
import numpy as np
import gc
import math

from tool.algorithm.algtools.oh2004 import oh2004

logger = logging.getLogger("mylog")


class RoughnessAlg:
    def __init__(self,):
        pass

    @staticmethod
    def cal_dielectric_const_factor(dielectric_const_factor_path, dielectric_const_path, arrival_angle_path,
                                    polarization_mode):
        """
       :param dielectric_const_factor_path:土壤介电常数因子影像路径
       :param dielectric_const_path:土壤介电常数
       :param arrival_angle_path:入射角影像文件
       :param polarization_mode:输入“HH”或者"VV"
       :return: True or False
       """
        image_handler = ImageHandler()
        proj = image_handler.get_projection(dielectric_const_path)
        geotrans = image_handler.get_geotransform(dielectric_const_path)
        dielectric_const = image_handler.get_band_array(dielectric_const_path, 1)
        angle_array = image_handler.get_band_array(arrival_angle_path, 1)

        where_dcf_1 = np.where(dielectric_const <= 1)
        dielectric_const[where_dcf_1] = 1
        try:
            if polarization_mode == "HH":
                tmp = np.power(np.sin(angle_array), 2)
                tmp = np.power(dielectric_const - tmp, 0.5)
                tmp = np.power(np.cos(angle_array) + tmp, 2)
                tmp = (dielectric_const - 1)/tmp
                tmp = np.power(tmp, 2)
                tmp[where_dcf_1] = 0 #异常值设为0
                image_handler.write_img(dielectric_const_factor_path, proj, geotrans, tmp)
            elif polarization_mode == "VV":
                tmp = -1 * (dielectric_const - 1) * (dielectric_const * (1 + np.power(np.sin(angle_array), 2)
                                                                         - np.power(np.sin(angle_array), 2)))
                tmp1 = np.power(np.sin(angle_array), 2)
                tmp1 = np.power(dielectric_const - tmp1, 0.5)
                tmp1 = np.power(dielectric_const * np.cos(angle_array) + tmp1, 2)
                tmp = np.power(tmp/tmp1, 2)
                tmp[where_dcf_1] = 0 #异常值设为0
                image_handler.write_img(dielectric_const_factor_path, proj, geotrans, tmp)
        except Exception:
            raise Exception("cal_dielectric_const_factor error!")

    @staticmethod
    def cal_soil_roughness(soil_roughness_path, tif_path, dielectric_const_factor_path, arrival_angle_path, radar_center_frequency):
        """
       :param soil_roughness_path:土壤粗糙度产品
       :param tif_path:极化影像路径
       :param dielectric_const_factor_path:土壤介电常数因子影像路径
       :param arrival_angle_path:入射角影像文件
       :param radar_center_frequency:微波波长
       :return: True or False
       """
        image_handler = ImageHandler()
        proj = image_handler.get_projection(tif_path)
        geotrans = image_handler.get_geotransform(tif_path)
        if radar_center_frequency == 0.0:
            raise ValueError (" RadarCenterFrequency == 0.0!")
        try:
            angle_array = image_handler.get_band_array(arrival_angle_path, 1)
            imaging_para = 8 * np.power(2 * math.pi / radar_center_frequency, 4) * np.power(np.cos(angle_array), 4)
            where_imaging_0 = np.where(imaging_para == 0.0)
            imaging_para[where_imaging_0] = 1

            tif_array = image_handler.get_band_array(tif_path, 1)
            dcf_array = image_handler.get_band_array(dielectric_const_factor_path, 1)
            where_dcf_0 = np.where(dcf_array == 0.0)
            dcf_array[where_dcf_0] = 1

            roughness = tif_array/(imaging_para * dcf_array)
            roughness[where_dcf_0] = 0
            roughness[where_imaging_0] = 0
        except Exception:
            raise Exception("cal_soil_roughness error!")
        image_handler.write_img(soil_roughness_path, proj, geotrans, roughness)


    @staticmethod
    def cal_soil_roug(surface_roug_path, hh_bsc_path, vv_bsc_path, hv_bsc_path, vh_bsc_path, angle_path, mask_path, wl):
        """

        :param surface_roug_path: 地表粗糙度
        :param hh_bsc_path:   hh极化路径
        :param vv_bsc_path:   vv极化路径
        :param hv_bsc_path:   hv极化路径
        :param vh_bsc_path:   vh极化路径
        :param angle_path:    入射角
        :param wl:  波长
        :return:
        """

        image_handler = ImageHandler()
        proj = image_handler.get_projection(hh_bsc_path)
        geotrans = image_handler.get_geotransform(hh_bsc_path)


        hh_arr = ImageHandler.get_band_array(hh_bsc_path, 1)
        hv_arr = ImageHandler.get_band_array(hv_bsc_path, 1)
        vh_arr = ImageHandler.get_band_array(vh_bsc_path, 1)
        vv_arr = ImageHandler.get_band_array(vv_bsc_path, 1)
        angle_arr = ImageHandler.get_band_array(angle_path, 1)
        mask_arr = ImageHandler.get_band_array(mask_path, 1)
        f = oh2004.lamda2freq(wl / 100) / 1e9       # GHZ
        n = np.array(hh_arr).shape[0] * np.array(hh_arr).shape[1]

        hh = np.array(hh_arr).flatten().astype(np.double)
        hh[np.where(hh == -9999)] = np.nan
        hv = np.array(hv_arr).flatten().astype(np.double)
        hv[np.where(hv == -9999)] = np.nan
        vv = np.array(vv_arr).flatten().astype(np.double)
        vv[np.where(vv == -9999)] = np.nan
        angle = np.array(angle_arr).flatten().astype(np.double)
        angle[np.where(angle == -9999)] = np.nan
        mask = np.array(mask_arr).flatten().astype(np.int32)

        mv = np.zeros(np.array(hh_arr).shape, dtype=np.double).flatten()
        h = np.zeros(np.array(hh_arr).shape, dtype=np.double).flatten()
        angle = angle * 180 / 3.1415926

        oh2004.retrieve_oh2004_main(n, mv, h, mask, vv, hh, hv, hv, angle, f)

        roughness = h.reshape(np.array(hh_arr).shape)

        image_handler.write_img(surface_roug_path, proj, geotrans, roughness)
        return True