gma中计算CWDI（作物水分亏缺指数）的源代码

这次是干货作物水分亏缺指数作物水分亏缺指数Crop Water Deficit IndexCWDI%从农田水分平衡出发引入了作物系数考虑了作物需水特性能很好好的反应作物缺水状况。计算公式如下源代码可以自由的适配和改造代码构建过程采用了向量化的思路进而提高计算效率并减少内存消耗对于巨量数据有效fromgma._algos.arrmtimportto_num_arrays_with_same_shapeimportnumpyasnp### to_num_arrays_with_same_shape 用于将数据数据格式化为统一的形状且保证为数字型整数、浮点数等数组。### 可选择性移除这个过程因为不是必须的。### 复杂度高是应为适应了n维数据按照不同轴axis计算classMoistureIndex: 降水-蒸散 构建的相关指数 def__init__(self,pre,et0,axisNone):## 初始化计算数据self._pre,self._et0to_num_arrays_with_same_shape(pre,et0)## 脱离 gma 不验证和处理原始数据直接使用下行代码#### self._pre, self._et0 pre, et0ifaxisisNone:self._preself._pre.flatten()self._et0self._et0.flatten()self._base_axis0else:self._base_axisaxis self._preself._pre.swapaxes(0,self._base_axis)self._et0self._et0.swapaxes(0,self._base_axis)self._shapeself._pre.shapedefCWDI(self,weights[0.3,0.25,0.2,0.15,0.1],duration_per_weight10):# etc 在这里是 self._et0## 确定每个计算单元所需的数据长度0轴上x_lenlen(weights)*duration_per_weight# 注意axis上前x_len-1是无数据的为np.nan因为不满足累积量weightsnp.array(weights)## 输出数组resnp.full(self._shape,np.nan)## 循环 x_len 次计算每个位置的结果foriinrange(x_len):## 每 x_len 个数据为一组末尾不足 x_len 的数据本轮计算时被丢弃d_num(self._shape[0]-i)%x_len d_endself._shape[0]-d_num c_pre,c_et0self._pre[i:d_end],self._et0[i:d_end]## 每 duration_per_weight 个数据为一组计算累加值以及i_cwdinew_shape(c_pre.shape[0]//duration_per_weight,duration_per_weight,*c_pre.shape[1:])c_prec_pre.reshape(new_shape).sum(axis1)c_et0c_et0.reshape(new_shape).sum(axis1)i_cwdi(1-c_pre/c_et0)*100i_cwdi[i_cwdi0]0## 每 len(weights) 个数据为一组计算 cwdii_cwdii_cwdi.reshape(i_cwdi.shape[0]//len(weights),len(weights),*i_cwdi.shape[1:])d_weightsweights[None,:,*([None]*(i_cwdi.ndim-2))]cwdi(i_cwdi*d_weights).sum(axis1)res[(ix_len-1)::x_len]cwdireturnres.swapaxes(0,self._base_axis)示例从 excel 开始计算in_filePRE_ET0.xlsx# 使用pandasimportpandasaspd dfpd.read_excel(in_file)pre,etcdf[PRE],df[ET0]##这里仅做演示实际使用时请使用真实的ETcmiMoistureIndex(pre,etc)cwdimi.CWDI()## 使用 gma3.0.0a12# from gma import gio, climet# xlsx_ly gio.open_vector(in_file)# df xlsx_ly.to_pandas()# pre, etc df[PRE], df[ET0] ##这里仅做演示实际使用时请使用真实的ETc# cwdi climet.index.CWDI(pre, etc)从 tif 开始计算## 使用 gma3.0.0a12fromgmaimportgio,climet ds_pregio.open_raster(rD:\BaiduNetdiskDownload\ne\第5章\PRE_Luoyang_1981-2020.tif)ds_etcgio.open_raster(rD:\BaiduNetdiskDownload\ne\第5章\ET0_Luoyang_1981-2020.tif)##这里仅做演示实际使用时请使用真实的ETc## 将两个 tif 组合为虚拟文件方便整体计算dsgio.VirtualRasterDataset([ds_pre,ds_etc])defcal(in_ar):# 一共 960 个波段其中前480个是 pre后480个是 etcprein_ar[:480]etcin_ar[480:]miMoistureIndex(pre,etc,axis0)cwdimi.CWDI()## 或者直接使用 gma 函数# cwdi climet.index.CWDI(pre, etc)cwdi[pre.mask]ds.nodata# 掩膜掉 nodatreturncwdi# algebraic 方法直接返回一个栅格数据集默认存储在内存中ndsds.algebraic(cal)### 或者直接保存到文件# out_file cwdi.tif# nds ds.algebraic(cal, out_dst out_file)获取本文用到的数据私信作者