gma中计算CWDI（作物水分亏缺指数）的源代码
这次是干货
作物水分亏缺指数
作物水分亏缺指数（Crop Water Deficit Index，CWDI，%）从农田水分平衡出发，引入了作物系数，考虑了作物需水特性，能很好好的反应作物缺水状况。计算公式如下：
源代码（可以自由的适配和改造）
代码构建过程采用了向量化的思路，进而提高计算效率并减少内存消耗（对于巨量数据有效）
fromgma._algos.arrmtimportto_num_arrays_with_same_shapeimportnumpyasnp### to_num_arrays_with_same_shape 用于将数据数据格式化为统一的形状，且保证为数字型（整数、浮点数等）数组。### 可选择性移除这个过程，因为不是必须的。### 复杂度高是应为适应了n维数据按照不同轴（axis）计算classMoistureIndex:''' 降水-蒸散 构建的相关指数 '''def__init__(self,pre,et0,axis=None):## 初始化计算数据self._pre,self._et0=to_num_arrays_with_same_shape(pre,et0)## 脱离 gma （不验证和处理原始数据）直接使用下行代码#### self._pre, self._et0 = pre, et0ifaxisisNone:self._pre=self._pre.flatten()self._et0=self._et0.flatten()self._base_axis=0else:self._base_axis=axis self._pre=self._pre.swapaxes(0,self._base_axis)self._et0=self._et0.swapaxes(0,self._base_axis)self._shape=self._pre.shapedefCWDI(self,weights=[0.3,0.25,0.2,0.15,0.1],duration_per_weight=10):# etc 在这里是 self._et0## 确定每个计算单元所需的数据长度（0轴上）x_len=len(weights)*duration_per_weight# 注意：axis上前x_len-1是无数据的（为np.nan），因为不满足累积量weights=np.array(weights)## 输出数组res=np.full(self._shape,np.nan)## 循环 x_len 次，计算每个位置的结果foriinrange(x_len):## 每 x_len 个数据为一组，末尾不足 x_len 的数据本轮计算时被丢弃d_num=(self._shape[0]-i)%x_len d_end=self._shape[0]-d_num c_pre,c_et0=self._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_pre=c_pre.reshape(new_shape).sum(axis=1)c_et0=c_et0.reshape(new_shape).sum(axis=1)i_cwdi=(1-c_pre/c_et0)*100i_cwdi[i_cwdi<0]=0## 每 len(weights) 个数据为一组，计算 cwdii_cwdi=i_cwdi.reshape(i_cwdi.shape[0]//len(weights),len(weights),*i_cwdi.shape[1:])d_weights=weights[None,:,*([None]*(i_cwdi.ndim-2))]cwdi=(i_cwdi*d_weights).sum(axis=1)res[(i+x_len-1)::x_len]=cwdireturnres.swapaxes(0,self._base_axis)示例
从 excel 开始计算
in_file="PRE_ET0.xlsx"# 使用pandasimportpandasaspd df=pd.read_excel(in_file)pre,etc=df['PRE'],df['ET0']##这里仅做演示，实际使用时请使用真实的ETcmi=MoistureIndex(pre,etc)cwdi=mi.CWDI()## 使用 gma==3.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 开始计算
## 使用 gma==3.0.0a12fromgmaimportgio,climet ds_pre=gio.open_raster(r"D:\BaiduNetdiskDownload\ne\第5章\PRE_Luoyang_1981-2020.tif")ds_etc=gio.open_raster(r"D:\BaiduNetdiskDownload\ne\第5章\ET0_Luoyang_1981-2020.tif")##这里仅做演示，实际使用时请使用真实的ETc## 将两个 tif 组合为虚拟文件，方便整体计算ds=gio.VirtualRasterDataset([ds_pre,ds_etc])defcal(in_ar):# 一共 960 个波段，其中前480个是 pre，后480个是 etcpre=in_ar[:480]etc=in_ar[480:]mi=MoistureIndex(pre,etc,axis=0)cwdi=mi.CWDI()## 或者直接使用 gma 函数# cwdi = climet.index.CWDI(pre, etc)cwdi[pre.mask]=ds.nodata# 掩膜掉 nodatreturncwdi# algebraic 方法直接返回一个栅格数据集，默认存储在内存中nds=ds.algebraic(cal)### 或者直接保存到文件# out_file = 'cwdi.tif'# nds = ds.algebraic(cal, out_dst = out_file)