# -*- coding: UTF-8 -*- """ @Project:microproduct @File:SoilMoistureALg.py @Function:实现土壤水分计算的算法 @Contact: @Author:SHJ @Date:2021/8/10 10:01 @Version:1.0.0 """ import logging from tool.algorithm.image.ImageHandle import ImageHandler import numpy as np import time import scipy import scipy.spatial.transform._rotation_groups # 解决打包的问题 from scipy.optimize import fmin import math from tool.algorithm.algtools.oh2004 import oh2004 logger = logging.getLogger("mylog") class MoistureAlg: def __init__(self,): pass @staticmethod def cal_vwc(vwc_path, ndwi_path, e1, e2): """ :param vwc_path: :param ndwi_path: :param e1: :param e2: :return: True or False """ image_handler = ImageHandler() proj = image_handler.get_projection(ndwi_path) geotrans = image_handler.get_geotransform(ndwi_path) array = image_handler.get_band_array(ndwi_path, 1) # 原方案的计算方式 vwc = e1 * np.square(array) + e2 * array # 论文《基于Radarsat-2全极化数据的高原牧草覆盖地表土壤水分反演》的计算方式 #vwc = e1 * array + e2 image_handler.write_img(vwc_path, proj, geotrans, vwc) return True @staticmethod def cal_bare_soil_bsc(bsc_path, tif_path, vwc_path, arrival_angle_path, c, d): """ :param bsc_path:裸土后向散射系数(线性值） :param tif_path:输入影像 会自动完成 dB 到线性值变换 :param vwc_path:vwc影像路径 :param arrival_angle_path:入射角影像文件（单位：弧度） :param c:经验系数 a :param d:经验系数 b """ image_handler = ImageHandler() proj = image_handler.get_projection(tif_path) geotrans = image_handler.get_geotransform(tif_path) tif_array = image_handler.get_band_array(tif_path, 1) vwc_array = image_handler.get_band_array(vwc_path, 1) angle_array = image_handler.get_band_array(arrival_angle_path, 1) tif_array=np.power(10.0,tif_array/10.0) # dB --> 线性值 后向散射系数 try: cos_angle = np.cos(angle_array) tmp1 = tif_array - c * vwc_array * cos_angle tmp2 = np.exp(-2 * d * vwc_array * (1/cos_angle)) bsc_array =c * vwc_array * cos_angle + tmp1/tmp2 # --- 公式计算错误 except BaseException as msg: logger.error(msg) return False bsc_array = np.where(bsc_array > 0, bsc_array, 0) bsc_array[tif_array==1] = 0 image_handler.write_img(bsc_path, proj, geotrans, bsc_array, '0') return True @staticmethod def cal_soil_moisture(soil_moisture_path, hh_bsc_path, vv_bsc_path, arrival_angle_path, mask_path, λ, f=5.3, T=24.5,bd=1, vsand=0.54, vclay=0.19,block_num=0): """ :param soil_moisture_path:土壤水分产品路径 :param hh_bsc_path:hh裸土后向散射系数 :param vv_bsc_path:vv裸土后向散射系数 :param arrival_angle_path:入射角影像文件 :param mask_path:掩膜影像文件 :param λ:经验系数 :param f:微波频率，单位GHz :param T:温度，摄氏度 :param bd:土壤容重 :param vsand:沙土含量，范围0-1 :param vclay:黏土含量，范围0-1 :return: True or False """ image_handler = ImageHandler() proj = image_handler.get_projection(hh_bsc_path) geotrans = image_handler.get_geotransform(hh_bsc_path) angle_array = image_handler.get_band_array(arrival_angle_path, 1) try: # 计算土壤介电常数 hh_bsc_array = image_handler.get_band_array(hh_bsc_path, 1) if np.any(hh_bsc_array < 0): logger.error("hh_bsc_array include negative value!") return False tmp = np.power(10, 0.19) * np.power(float(λ), 0.15) * hh_bsc_array # 处理异常值 where_tmp1_0 = np.where(tmp == 0.0) tmp[where_tmp1_0] = 1 except BaseException as msg: logger.error(msg) return False try: vv_bsc_array = image_handler.get_band_array(vv_bsc_path, 1) if np.any(vv_bsc_array < 0): logger.error("vv_bsc_array include negative value!") return False tmp2 = np.power(np.cos(angle_array), 1.82) * np.power(np.sin(angle_array), 0.93) *\ np.power(vv_bsc_array, 0.786) # 处理异常值 where_tmp2_0 = np.where(tmp2 == 0.0) tmp2[where_tmp2_0] = 1 tmp = tmp/tmp2 # 土壤介电常数 soil_dielectric = (1/(0.024 * np.tan(angle_array))) * np.log10(tmp) soil_dielectric_path = soil_moisture_path.replace("soil_moisture","soil_dielectric") image_handler.write_img(soil_dielectric_path, proj, geotrans, soil_dielectric) except BaseException as msg: logger.error(msg) return False # mask_array = ImageHandler.get_band_array(mask_path, 1) # soil_dielectric[np.where(mask_array == 0)] = np.nan try: # Dobsen模型计算土壤水分 # soil_moisture = dobsen_inverse(f, T, bd, vsand, vclay, soil_dielectric, block_num) # topp模型计算土壤水分 soil_moisture = -5.3 * np.power(10.0, -2) + 2.92 * np.power(10.0, -2) * soil_dielectric - \ 5.5 * np.power(10.0, -4) * np.power(soil_dielectric, 2) + \ 4.3 * np.power(10.0, -6) * np.power(soil_dielectric, 3) # 处理异常值 # soil_moisture[where_tmp1_0] = 0 # soil_moisture[where_tmp2_0] = 0 except BaseException as msg: logger.error(msg) return False image_handler.write_img(soil_moisture_path, proj, geotrans, soil_moisture) return True @staticmethod def cal_soilM(soil_moisture_path, hh_bsc_path, vv_bsc_path, hv_bsc_path, vh_bsc_path, angle_path,mask_path, wl): """ :param soil_moisture_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 # wl 原来是cm ，得转成m 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) # foo_retrieve = inverse_radar(hh_arr, hv_arr, vh_arr, vv_arr, angle_arr, wl) # mv, h = foo_retrieve.retrieve_oh2004_Cython() 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) soil = np.array(mv).reshape(np.array(hh_arr).shape) image_handler.write_img(soil_moisture_path, proj, geotrans, soil) return True @staticmethod def cal_roughness_oh2004(surface_roughness_path, hh_bsc_path, vv_bsc_path, hv_bsc_path, angle_path, Freq, thres1, thres2): Nlig = ImageHandler.get_img_height(hh_bsc_path) Ncol = ImageHandler.get_img_width(hh_bsc_path) theta = ImageHandler.get_data(angle_path) Shvhv = ImageHandler.get_data(hv_bsc_path) Shhhh = ImageHandler.get_data(hh_bsc_path) Svvvv = ImageHandler.get_data(vv_bsc_path) im_proj, im_geotrans, _ = ImageHandler.read_img(hh_bsc_path) msk_out = np.zeros((Nlig, Ncol), dtype=np.float32) mv_oh = np.zeros((Nlig, Ncol), dtype=np.float32) s_oh = np.zeros((Nlig, Ncol), dtype=np.float32) msk_valid = np.zeros((Nlig, Ncol), dtype=np.float32) mv_inv1 = np.zeros((Nlig, Ncol), dtype=np.float32) mv_inv2 = np.zeros((Nlig, Ncol), dtype=np.float32) mv_inv3 = np.zeros((Nlig, Ncol), dtype=np.float32) ks_inv1 = np.zeros((Nlig, Ncol), dtype=np.float32) ks_inv2 = np.zeros((Nlig, Ncol), dtype=np.float32) for ii in range(Nlig): for jj in range(Ncol): f = 1e9 * Freq # p_max = 1 - math.exp(0.35 * math.log(theta[ii, jj] / (math.pi / 2)) * (thres1 ** -0.65)) * math.exp( # -0.4 * ((2 * math.pi / (3e8 / f) * thres2) ** 1.4)) msk_valid[ii, jj] = 1 mv_inv = 0.5 for tt in range(20): xx = 1. - Shvhv[ii, jj] / (0.11 * (mv_inv ** 0.7) * (math.cos(theta[ii, jj]) ** 2.2)) if xx <= 0: mv_inv = 0 break else: xx = -3.125 * math.log(xx) if mv_inv <= 0 or xx <= 0: mv_inv = 0 break else: G = (1 - math.exp( 0.35 * math.log(theta[ii, jj] / (math.pi / 2)) * (mv_inv ** -0.65)) * math.exp( -0.4 * (xx ** (0.556 * 1.4))) - Shhhh[ii, jj] / Svvvv[ii, jj]) G1 = -math.exp( 0.35 * math.log(theta[ii, jj] / (math.pi / 2)) * (mv_inv ** -0.65)) * math.log( theta[ii, jj] / (math.pi / 2)) * 0.35 * (-0.65) * (mv_inv ** -1.65) * math.exp( -0.4 * (xx ** (0.556 * 1.4))) G1 -= math.exp( 0.35 * math.log(theta[ii, jj] / (math.pi / 2)) * (mv_inv ** -0.65)) * math.exp( -0.4 * (xx ** (0.556 * 1.4))) * (-0.4) * 0.556 * 1.4 * (xx ** (0.556 * 1.4 - 1)) * ( -3.125) / (1 - Shvhv[ii, jj] / ( 0.11 * (mv_inv ** 0.7) * (math.cos(theta[ii, jj]) ** 2.2))) * 0.7 * Shvhv[ ii, jj] / 0.11 / (math.cos(theta[ii, jj]) ** 2.2) * (mv_inv ** -1.7) if abs(G / G1) < 1e-10: break else: if G1 != 0: mv_inv -= G / G1 else: continue mv_inv1[ii, jj] = mv_inv if mv_inv == 0: ks_inv1[ii, jj] = 0 else: xx = 1. - Shvhv[ii, jj] / (0.11 * (mv_inv1[ii, jj] ** 0.7) * (math.cos(theta[ii, jj]) ** 2.2)) if xx <= 0: ks_inv1[ii, jj] = 0 else: xx = -3.125 * math.log(xx) if xx <= 0: ks_inv1[ii, jj] = 0 else: ks_inv1[ii, jj] = xx ** 0.556 xx = 1. - (Shvhv[ii, jj] / Svvvv[ii, jj]) / (0.095 * np.power(0.13 + np.sin(1.5 * theta[ii, jj]), 1.4)) if xx <= 0: ks_inv2[ii, jj] = 0 else: xx = np.log(xx) / (-1.3) ks_inv2[ii, jj] = np.power(xx, 10. / 9.) if ks_inv2[ii, jj] == 0: mv_inv2[ii, jj] = 0 else: xx = (1 - Shhhh[ii, jj] / Svvvv[ii, jj]) / np.exp(-0.4 * np.power(ks_inv2[ii, jj], 1.4)) if xx <= 0: mv_inv2[ii, jj] = 0 else: xx = np.log(xx) / (0.35 * np.log(theta[ii, jj] / (np.pi / 2))) mv_inv2[ii, jj] = np.power(xx, -100. / 65.) if ks_inv2[ii, jj] == 0: mv_inv3[ii, jj] = 0 else: xx = Shvhv[ii, jj] / (0.11 * np.power(np.cos(theta[ii, jj]), 2.2) * ( 1 - np.exp(-0.32 * np.power(ks_inv2[ii, jj], 1.8)))) if xx <= 0: mv_inv3[ii, jj] = 0 else: mv_inv3[ii, jj] = np.power(xx, 10. / 7.) a = 1 b = 1 c = 1 if a == 0 and b == 0 and c == 0: mv_inv = 0 else: mv_inv = (a * mv_inv1[ii, jj] + b * mv_inv2[ii, jj] + c * mv_inv3[ii, jj]) / (a + b + c) m = 1 n = 1 if m == 0 and n == 0: ks_inv = 0 else: ks_inv = (m * ks_inv1[ii, jj] + 0.25 * n * ks_inv2[ii, jj]) / (m + 0.25 * n) msk_out[ii, jj] = 1 mv_oh[ii, jj] = mv_inv * msk_out[ii, jj] * 100. s_oh[ii, jj] = ks_inv * msk_out[ii, jj] / (2 * np.pi / (3e8 / f)) * 100. ImageHandler.write_img(surface_roughness_path, im_proj, im_geotrans, s_oh, '0') def dobsen_inverse(f, T, bd, vsand, vclay, soil_dielectric_array,block_num, x0=0): """ Dobsen模型，土壤水分反演 :param f:微波频率，单位GHz :param T:温度，摄氏度 :param bd:土壤容重 :param vsand:沙土含量，范围0-1 :param vclay:黏土含量，范围0-1 :param soil_dielectric_array:土壤介电常数 :param x0 :土壤水分寻优，初始值 :return: True or False """ alpha = 0.65 sd = 2.65 dcs = (1.01 + 0.44 * sd) ** 2 - 0.062 dc0 = 0.008854 dcw0 = 88.045 - 0.4147 * T + 6.295e-4 * (T ** 2) + 1.075e-5 * (T ** 3) tpt = 0.11109 - 3.824e-3 * T + 6.938e-5 * (T ** 2) - 5.096e-7 * (T ** 3) dcwinf = 4.9 if f >= 1.4: sigma = -1.645 + 1.939 * bd - 2.013e-2 * vsand + 1.594e-2 * vclay else: sigma = 0.0467 + 0.2204 * bd - 4.111e-3 * vsand + 6.614e-3 * vclay dcwr = dcwinf + ((dcw0 - dcwinf) / (1 + (tpt * f) ** 2)) betar = 1.2748 - 0.00519 * vsand - 0.00152 * vclay betai = 1.33797 - 0.00603 * vsand - 0.00166 * vclay # fun = lambda vwc: np.abs(np.sqrt( # ((1.0 + (bd / sd) * ((dcs ** alpha) - 1.0) + (vwc ** betar) * (dcwr ** alpha) - vwc) ** (1 / alpha)) ** 2 + # ((vwc ** (betai / alpha)) * ((tpt * f * (dcw0 - dcwinf)) / (1 + (tpt * f) ** 2) + sigma * (1.0 - (bd / sd)) / (2 * np.pi * dc0 * f * vwc))) ** 2 # ) - soil_dielectric) fun = lambda vwc: np.abs(np.sqrt( ((1.0 + (bd / sd) * ((dcs ** alpha) - 1.0) + (vwc ** betar) * (dcwr ** alpha) - vwc) ** (1 / alpha)) ** 2 + ((vwc ** (betai / alpha)) * ((tpt * f * (dcw0 - dcwinf)) / (1 + (tpt * f) ** 2) + sigma * (1.0 - (bd / sd)) / (8 * math.atan(1.0) * dc0 * f * vwc))) ** 2 ) - soil_dielectric) soil_dielectric_array[np.isnan(soil_dielectric_array)] = 0 w = np.where((soil_dielectric_array == 0) == False) num = len(w[0]) soil_mois = np.zeros(num) bar_list = [0, int(0.1*num), int(0.2*num), int(0.3*num), int(0.4*num), int(0.5*num), int(0.6*num), int(0.7*num), int(0.8*num), int(0.9*num), num-1] start = time.perf_counter() for soil_dielectric, n in zip(soil_dielectric_array[w], range(num)): soil_mois[n] = fmin(fun, x0, disp=0) if n in bar_list: logger.info('block:{},cal soil_moisture proggress bar:{}%,use time :{}s'.format(block_num, round(n/num * 100), int(time.perf_counter()-start))) soil_dielectric_array[w] = soil_mois # 土壤水分 return soil_dielectric_array