# # 模型计算的库 # import cython cimport cython # 必须导入 import numpy as np cimport numpy as np from libc.math cimport pi from scipy.optimize import leastsq import random import logging logger = logging.getLogger("mylog") def WMCModel(param_arr,sample_lai,sample_soil,sample_inc,sample_sigma): """ WMC模型 增加 归一化植被指数 Args: param_arr (np.ndarray): 参数数组 sample_lai (double): 叶面积指数 sample_soil (double): 土壤含水量 sample_inc (double): 入射角（弧度值） sample_sigma (double): 后向散射系数（线性值） Returns: double: 方程值 """ # 映射参数，方便修改模型 A,B,C,D,M,N=param_arr # 在这里修改模型 V_lai=sample_lai #V_lai=E*sample_lai+F exp_gamma=np.exp(-2*B*((V_lai*D+C))*(1/np.cos(sample_inc))) sigma_soil=M*sample_soil+N sigma_veg=A*((V_lai))*np.cos(sample_inc) f_veg=1 result=sigma_veg*(1-exp_gamma)+sigma_soil*exp_gamma-sample_sigma return result def train_WMCmodel(lai_water_inc_sigma_list,params_X0,train_err_image_path,draw_flag=True): """ 训练模型参数 Args: lai_waiter_inc_sigma_list (list): 训练模型使用的样本呢 """ def f(X): eqs=[] for lai_water_inc_sigma_item in lai_water_inc_sigma_list: sample_lai=lai_water_inc_sigma_item[4] sample_sigma=lai_water_inc_sigma_item[5] # 5: csv_sigma, 8:tiff_sigma sample_soil=lai_water_inc_sigma_item[6] sample_inc=lai_water_inc_sigma_item[7] FVC=lai_water_inc_sigma_item[8] eqs.append(WMCModel(X,sample_lai,sample_soil,sample_inc,sample_sigma)) return eqs X0 = params_X0 # 初始值 # logger.info(str(X0)) h = leastsq(f, X0) # logger.info(h[0],h[1]) err_f=f(h[0]) x_arr=[lai_waiter_inc_sigma_item[4] for lai_waiter_inc_sigma_item in lai_water_inc_sigma_list] # 根据误差大小进行排序 # logger.info("训练集：\n根据误差输出点序\n数量：{}\n点序\t误差值\t 样点信息".format(str(np.array(err_f).shape))) # for i in np.argsort(np.array(err_f)): # logger.info('{}\t{}\t{}'.format(i,err_f[i],str(lai_water_inc_sigma_list[i]))) # logger.info("\n误差点序输出结束\n") if draw_flag: # logger.info(err_f) # logger.info(np.where(np.abs(err_f)<10)) from matplotlib import pyplot as plt plt.scatter(x_arr,err_f) plt.title("equation-err") plt.savefig(train_err_image_path,dpi=600) plt.show() return h[0] def test_WMCModel(lai_waiter_inc_sigma_list,param_arr,lai_X0,test_err_image_path,draw_flag=True): """ 测试模型训练结果 Args: lai_waiter_inc_sigma_list (list): 测试使用的样本集 A (_type_): 参数A B (_type_): 参数B C (_type_): 参数C D (_type_): 参数D M (_type_): 参数M N (_type_): 参数N lai_X0 (_type_): 初始值 Returns: list: 误差列表 [sample_lai,err,predict] """ err=[] err_f=[] x_arr=[] err_lai=[] for lai_waiter_inc_sigma_item in lai_waiter_inc_sigma_list: sample_time,sample_code,sample_lon,sample_lat,sample_lai,csv_sigma,sample_soil,sample_inc,sample_sigma=lai_waiter_inc_sigma_item def f(X): lai=X[0] eqs=[WMCModel(param_arr,lai,sample_soil,sample_inc,csv_sigma)] return eqs X0=lai_X0 h = leastsq(f, X0) temp_err=h[0]-sample_lai err_lai.append(temp_err[0]) # lai预测的插值 err.append([sample_lai,temp_err[0],h[0][0],sample_code]) err_f.append(f(h[0])[0]) # 方程差 x_arr.append(sample_lai) # 根据误差大小进行排序 # logger.info("测试集：\n根据误差输出点序\n数量：{}\n点序\t误差值\t 方程差\t样点信息".format(str(np.array(err_lai).shape))) # for i in np.argsort(np.array(err_lai)): # logger.info('{}\t{}\t{}\t{}'.format(i,err_lai[i],err_f[i],str(lai_waiter_inc_sigma_list[i]))) # logger.info("\n误差点序输出结束\n") if draw_flag: from matplotlib import pyplot as plt plt.scatter(x_arr,err_lai) plt.title("equation-err") plt.savefig(test_err_image_path,dpi=600) plt.show() return err def processs_WMCModel(param_arr,lai_X0,sigma,inc_angle,soil_water): if(sigma<0 ): return np.nan def f(X): lai=X[0] eqs=[WMCModel(param_arr,lai,soil_water,inc_angle,sigma )] return eqs h = leastsq(f, [lai_X0]) return h[0][0] # Cython 的扩展地址 cpdef np.ndarray[double,ndim=2] process_tiff(np.ndarray[double,ndim=2] sigma_tiff, np.ndarray[double,ndim=2] inc_tiff, np.ndarray[double,ndim=2] soil_water_tiff, np.ndarray[double,ndim=1] param_arr, double lai_X0): cdef np.ndarray[double,ndim=2] result=sigma_tiff cdef int param_arr_length=param_arr.shape[0] cdef int height=sigma_tiff.shape[0] cdef int width=sigma_tiff.shape[1] cdef int i=0 cdef int j=0 cdef double temp=0 while i=0 else np.nan result[i,j]=temp j=j+1 i=i+1 return result