前提:已经配置好了python环境和GDAL库,numpy库
实现过程:
1. 读取数据:读取需要的各种数据:(这里需要用到NDVI、温度、降水、太阳辐射数据)。以读取NDVI数据为例:
filename = r'G:\graduation\code\NDVI15.tif'
NDVI = gdal.Open(filename)
2. 将读取的数据转存到numpy数组
Data = NDVI.ReadAsArray(0,0,width_NDVI,height_NDVI)
3.然后开始FPAR的计算
FPAR_NDVI = (Data-miN)*(0.95-0.001)/(maN-miN)+0.001FPAR_SR = (SR-miR)*(0.95-0.001)/(maR-miR)+0.001
FPAR = (FPAR_NDVI+FPAR_SR)/2 ##植被层对入射光合有效辐射的吸收比例
4.继而进行APAR和NPP的计算
for xi in range(height_Veg): for yj in range(width_Veg): APAR[xi][yj] = RData[xi,yj]*Data2[xi,yj]*0.5 Te2[xi][yj] = 1.184/(1+math.exp(0.2*(Topt-10-Data3[xi,yj])))*1 \ /(1+math.exp(0.3*(-Topt-10+Data3[xi,yj]))) Te[xi][yj] = Te2[xi][yj]*Te1 Ep0[xi,yj] = 16*(10*Data3[xi,yj]/I)**a Rn[xi,yj] = ((Ep0[xi,yj]*Data4[xi,yj])**0.5)*(0.369+0.589*((Ep0[xi,yj]/Data4[xi,yj])**0.5)) E[xi,yj] = Data4[xi,yj]*Rn[xi,yj]*(Data4[xi,yj]**2+Rn[xi,yj]**2+Data4[xi,yj]*Rn[xi,yj]) \ /((Data4[xi,yj]+Rn[xi,yj])*(Data4[xi,yj]**2+Rn[xi,yj]**2)) Ep[xi,yj] = (E[xi,yj]+Ep0[xi,yj])/2 W[xi,yj] = 0.5+(0.5*E[xi,yj])/Ep[xi,yj] NPP[xi,yj] = Data2[xi,yj]*FPAR[xi,yj]*0.5*Te[xi][yj]*W[xi,yj]*0.541 #NPP