148 lines
6.1 KiB
Python
148 lines
6.1 KiB
Python
# -*- 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
|
|
|