diff --git a/surfaceRoughness_oh2004/SurfaceRoughnessAlg.py b/surfaceRoughness_oh2004/SurfaceRoughnessAlg.py new file mode 100644 index 00000000..58471a83 --- /dev/null +++ b/surfaceRoughness_oh2004/SurfaceRoughnessAlg.py @@ -0,0 +1,261 @@ +# -*- 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:输入影像 + :param vwc_path:vwc影像路径 + :param arrival_angle_path:入射角影像文件（单位：弧度） + :param c:经验系数 + :param d:经验系数 + """ + 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) + + 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 = 1 + tmp1/tmp2 + + except BaseException as msg: + logger.error(msg) + return False + bsc_array = np.where(bsc_array > 0, bsc_array, 0) + + image_handler.write_img(bsc_path, proj, geotrans, bsc_array) + 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 + +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 + + + + + diff --git a/surfaceRoughness_oh2004/SoilMoistureMain.py b/surfaceRoughness_oh2004/SurfaceRoughnessMain.py similarity index 91% rename from surfaceRoughness_oh2004/SoilMoistureMain.py rename to surfaceRoughness_oh2004/SurfaceRoughnessMain.py index 9b285920..723cdc2f 100644 --- a/surfaceRoughness_oh2004/SoilMoistureMain.py +++ b/surfaceRoughness_oh2004/SurfaceRoughnessMain.py @@ -19,16 +19,16 @@ from tool.algorithm.transforml1a.transHandle import TransImgL1A from tool.algorithm.xml.AlgXmlHandle import ManageAlgXML, CheckSource,InitPara # 导入xml文件读取与检查文件 from tool.algorithm.image.ImageHandle import ImageHandler from tool.algorithm.algtools.logHandler import LogHandler -from SoilMoistureALg import MoistureAlg as alg +from SurfaceRoughnessAlg import MoistureAlg as alg from tool.algorithm.block.blockprocess import BlockProcess from tool.algorithm.algtools.MetaDataHandler import MetaDataHandler from tool.config.ConfigeHandle import Config as cf from tool.algorithm.xml.CreatMetafile import CreateMetafile from tool.algorithm.algtools.ROIAlg import ROIAlg as roi -from SoilMoistureXmlInfo import CreateDict, CreateStadardXmlFile +from SurfaceRoughnessXmlInfo import CreateDict, CreateStadardXmlFile from tool.algorithm.algtools.filter.lee_Filter import Filter from tool.file.fileHandle import fileHandle -from SoilMoistureTool import SoilMoistureTool +from SurfaceRoughnessTool import SoilMoistureTool import logging import os import shutil @@ -114,6 +114,8 @@ class MoistureMain: # para_dic.update(InitPara.get_meta_dic(InitPara.get_meta_paths(file_dir, name), name)) para_dic.update(InitPara.get_meta_dic_new(InitPara.get_meta_paths(file_dir, name), name)) # tif路径字典 + parameter_path = os.path.join(file_dir, "orth_para.txt") + para_dic.update({"paraMeter": parameter_path}) pol_dic = InitPara.get_polarization_mode(InitPara.get_tif_paths(file_dir, name)) flag_list = [0, 0, 0, 0] @@ -136,6 +138,8 @@ class MoistureMain: para_dic.update({'inc_angle': in_tif_path}) elif 'ori_sim' == key: para_dic.update({'ori_sim': in_tif_path}) + elif 'sim_ori' == key: + para_dic.update({'sim_ori': in_tif_path}) elif 'LocalIncidenceAngle' == key: para_dic.update({'LocalIncidenceAngle': in_tif_path}) elif 'inci_Angle-ortho' == key: @@ -173,13 +177,13 @@ class MoistureMain: # self.__preprocessed_paras, scopes_roi = p.preprocessing_oh2004(para_names, self.__processing_paras, # self.__workspace_preprocessing_path, self.__workspace_preprocessed_path) - para_names_geo = ['Covering', 'NDVI'] + para_names_geo = ['Covering', 'NDVI', 'inc_angle', 'sim_ori'] p = pp() cutted_img_paths, scopes_roi = p.cut_geoimg(self.__workspace_preprocessing_path, para_names_geo, self.__processing_paras) self.__preprocessed_paras.update(cutted_img_paths) - para_names_l1a = ["HH", "VV", "HV", "VH", 'inci_Angle-ortho', 'ori_sim'] + para_names_l1a = ["HH", "VV", "HV", "VH", 'ori_sim'] self._tr = TransImgL1A(self.__processing_paras['ori_sim'], scopes_roi) for name in para_names_l1a: @@ -254,15 +258,14 @@ class MoistureMain: pp.resampling_by_scale(self.__preprocessed_paras["Covering"], cover_rampling_path, refer_img_path) self.__preprocessed_paras["Covering"] = cover_rampling_path - def inter_Range2Geo(self, lon_lat_path, data_tiff, grid_path, space): + def calInterpolation_bil_Wgs84_rc_sar_sigma(self, parameter_path, dem_rc, in_sar, out_sar): ''' - # std::cout << "mode 10"; - # std::cout << "SIMOrthoProgram.exe 10 lon_lat_path data_tiff grid_path space"; + # std::cout << "mode 11"; + # std::cout << "SIMOrthoProgram.exe 11 in_parameter_path in_rc_wgs84_path in_ori_sar_path out_orth_sar_path"; ''' exe_path = r".\baseTool\x64\Release\SIMOrthoProgram.exe" - exe_cmd = r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5}".format(exe_path, 10, - lon_lat_path, data_tiff, - grid_path, space) + exe_cmd = r"set PROJ_LIB=.\baseTool\x64\Release; & {0} {1} {2} {3} {4} {5}".format(exe_path, 11, parameter_path, + dem_rc, in_sar, out_sar) print(exe_cmd) print(os.system(exe_cmd)) print("==========================================================================") @@ -275,7 +278,7 @@ class MoistureMain: tem_folder = self.__workspace_path + EXE_NAME + r"\Temporary""\\" soilOh2004 = SoilMoistureTool(self.__workspace_preprocessed_path, self.__workspace_processing_path, self.__cols, - self.__rows, self.__preprocessed_paras['inci_Angle-ortho'], self.__processing_paras['Origin_META']) + self.__rows, self.__preprocessed_paras['inc_angle'], self.__processing_paras['Origin_META']) result = soilOh2004.soil_oh2004() logger.info('progress bar: 80%') @@ -287,7 +290,10 @@ class MoistureMain: shutil.copy(tif_file, product_temp_path) product_geo_path = os.path.join(tem_folder, 'SurfaceRoughnessProduct_geo.tif') - self.inter_Range2Geo(self.__preprocessed_paras['ori_sim'], product_temp_path, product_geo_path, pixelspace) + self.calInterpolation_bil_Wgs84_rc_sar_sigma(self.__processing_paras['paraMeter'], + self.__preprocessed_paras['sim_ori'], product_temp_path, + product_geo_path) + # self.inter_Range2Geo(self.__preprocessed_paras['ori_sim'], product_temp_path, product_geo_path, pixelspace) # self._tr.l1a_2_geo_int(self.__preprocessed_paras['ori_sim'], product_temp_path, product_geo_path, 'linear') # # hh_geo_path = self.__workspace_processing_path + "hh_geo.tif" @@ -298,10 +304,11 @@ class MoistureMain: bare_land_mask_path = roi().roi_process(para_names, self.__workspace_processing_path + "/roi/", self.__processing_paras, self.__preprocessed_paras) - product_path = os.path.join(self.__product_dic, 'SurfaceRoughnessProduct.tif') + SrcImageName = os.path.split(self.__input_paras["DualPolarSAR"]['ParaValue'])[1].split('.tar.gz')[0] + '-Roughness.tif' + product_path = os.path.join(self.__product_dic, SrcImageName) # 获取影像roi区域 - roi.cal_roi(product_path, product_geo_path, bare_land_mask_path, background_value=np.nan) + roi.cal_roi(product_path, product_geo_path, bare_land_mask_path, background_value=0) # 生成快视图 self.imageHandler.write_quick_view(product_path) diff --git a/surfaceRoughness_oh2004/SoilMoistureMain.spec b/surfaceRoughness_oh2004/SurfaceRoughnessMain.spec similarity index 92% rename from surfaceRoughness_oh2004/SoilMoistureMain.spec rename to surfaceRoughness_oh2004/SurfaceRoughnessMain.spec index 9478bbea..f1b5fc78 100644 --- a/surfaceRoughness_oh2004/SoilMoistureMain.spec +++ b/surfaceRoughness_oh2004/SurfaceRoughnessMain.spec @@ -3,7 +3,7 @@ block_cipher = None -a = Analysis(['SoilMoistureMain.py'], +a = Analysis(['SurfaceRoughnessMain.py'], pathex=['D:\\estar-proj\\microproduct\\surfaceRoughness_oh2004'], binaries=[], datas=[('D:/Anaconda/envs/micro/Lib/site-packages/dask/dask.yaml', './dask'), ('D:/Anaconda/envs/micro/Lib/site-packages/distributed/distributed.yaml', './distributed')], @@ -23,7 +23,7 @@ exe = EXE(pyz, a.zipfiles, a.datas, [], - name='SoilMoistureMain', + name='SurfaceRoughnessMain', debug=False, bootloader_ignore_signals=False, strip=False, diff --git a/surfaceRoughness_oh2004/SoilMoistureTool.py b/surfaceRoughness_oh2004/SurfaceRoughnessTool.py similarity index 99% rename from surfaceRoughness_oh2004/SoilMoistureTool.py rename to surfaceRoughness_oh2004/SurfaceRoughnessTool.py index 09a7fc7d..d29d0d10 100644 --- a/surfaceRoughness_oh2004/SoilMoistureTool.py +++ b/surfaceRoughness_oh2004/SurfaceRoughnessTool.py @@ -53,7 +53,7 @@ class SoilMoistureTool: leeFilter.api_lee_refined_filter_T3('', t3_path, lee_filter_path, 0, 0, atp.rows(), atp.cols()) logger.info("refine_lee filter success!") # logging.info("refine_lee filter success!") - return t3_path + return lee_filter_path # def create_incidence(self, incidence_xml): diff --git a/surfaceRoughness_oh2004/SoilMoistureXmlInfo.py b/surfaceRoughness_oh2004/SurfaceRoughnessXmlInfo.py similarity index 100% rename from surfaceRoughness_oh2004/SoilMoistureXmlInfo.py rename to surfaceRoughness_oh2004/SurfaceRoughnessXmlInfo.py