CASA方法估算NPP（IDL+ENVI）

一、NPP的定义
二、估算流程
三、实现代码（IDL)

一、NPP的定义

NPP是什么：净第一性生产力，又称净初级生产力，是指从绿色植物在单位时间内单位面积上能固定的总能量中除去其自身呼吸消耗掉的部分，单位为焦耳/厘米2/年或克/米2/年。

可表示为：NPP = GPP - Ra

在上式中Ra表示自养呼吸的消耗量、GPP表示总初级生产力。

意义：NPP反映着植物固定和转化光合作用产物的效率，也决定了可供利用的物质和能量。是人类生活所需食物、原料及燃料的来源。植物通过光合作用将太阳能固定并转化为植物生物量。

以上摘自百度百科

二、估算流程

CASA模型反演NPP流程图：

1、温度胁迫因子1

Topt(x,y) 为某一区域一年内 NDVI 值到达最高时的当月平均气温（单位℃），当某个月的平均气温低于或等于-10℃，Tε1(x,y) 取值为0，此时不会发生光合作用。

2、温度胁迫因子2

当某个月平均温度 T(x,y,t) 比适宜温度 Topt(x,y) 高10℃或者低13℃时，该月的 Tε2(x,y,t) 等于月平均温度为最适宜温度 Topt(x,y) 时 Tε2(x,y,t) 的一半。

3、水分胁迫因子

EET(x,y,t) 表示(x,y)处像元t月份的实际蒸散发量（单位mm），PET(x,y,t) 表示(x,y)处像元t月份的潜在蒸散发量（单位mm）。

P(x,y,t)表示(x,y)处像元 t月份的降水（单位mm），Rn(x,y,t) 表示(x,y)处像元t月份的地表净辐射量（单位MJ/m2）。

Ep0(x,y,t) 表示局地潜在蒸发量（单位mm）。

a(x,y) 为因地而异的常数。

I(x,y) 代表12个月总和热量指标。

T(x,y,t) 为 (x,y) 处像元 t 月份的平均温度。

4、实际光能利用率

式中，Tε1(x,y) 和 Tε2(x,y) 为温度胁迫因子，Wε(x,y) 为水分胁迫因子，εMAX 为最大光能利用率。

5、植被层对入射光和有效辐射的吸收比例

NDVImax 和 NDVImin 代表草地的 NDVI 最大值和最小值，FPARmin 和 FPARmax 的取值与植被类型无关，分别为0.001 和0.95。

6、光合有效辐射

SOL(x,y,t) 表示在 t月份空间位置 (x,y) 处的像元的太阳总辐射能量（单位MJ/m2）；FAPAR(x,y,t)表示植被层对入射光和有效辐射的吸收比例。

7、净初级生产力

APAR(x,y,t) 表示在空间位置 (x,y) 处的像元在 t月份吸收的光合有效辐射（单位MJ/m2/月），ε(x,y,t) 表示在空间位置 (x,y) 处的像元在 t月份的实际光能利用率（单位g.C/MJ）。

三、实现代码（IDL)

;*************************************************************************************************************
;   基于MODIS数据计算NPP（CASA)
;
;   将插值后的数据（包括温度、降水、辐射）放入同一个文件夹并命名为Rain_XX,Temp_XX,Radiate_XX,格式为.tif
;   将裁剪后的NDVI数据放入一个文件夹下，文件格式变为.dat
;
;   编写：鸽子
;   日期：2020/7/14
;*************************************************************************************************************PRO cal_NPPe=ENVI();读取数据work_dir=DIALOG_PICKFILE(title='请选择文件路径',/directory);获取文件路径CD,work_dirfns = FILE_SEARCH(work_dir,'*.tif',count = fnums);获取文件大小data=READ_TIFF(fns[0])data_size=SIZE(data);读取太阳净辐射数据data_radiate=MAKE_ARRAY(12,data_size[1],data_size[2],type=data_size[3])PRINT,SIZE(data_radiate)FOR i=0,12-1 DO BEGINENVI_OPEN_FILE,fns[i],r_fid=fidENVI_FILE_QUERY,fid,dims=dims,interleave=interleave,offset=offsetmap_info=ENVI_GET_MAP_INFO(fid=fid)data_radiate[i,*,*]=ENVI_GET_DATA(fid=fid,dims=dims,pos=0);data_radiate[i,*,*]=read_tiff(fns[i]);HELP,data_radiate[i,*,*]ENVI_FILE_MNG,id=fid,/removeENDFORdata_radiate = data_radiate/100 ;单位转换radiate_6_10=total(data_radiate[5:9,*,*],1,/NAN)o_fn=DIALOG_PICKFILE(title='radiate_6_10s计算结果保存为')ENVI_WRITE_ENVI_FILE,radiate_6_10,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=1,dsta_type=4,interleave=2,offset=offset,map_info=map_info;读取降雨量数据data_rain = MAKE_ARRAY(12,data_size[1],data_size[2],type=data_size[3]);全年总降水量rain_year = MAKE_ARRAY(1,data_size[1],data_size[2],type=data_size[3])PRINT,SIZE(data_rain)FOR i=12,24-1 DO BEGINENVI_OPEN_FILE,fns[i],r_fid=fidENVI_FILE_QUERY,fid,dims=dimsdata_rain[i-12,*,*] = ENVI_GET_DATA(fid=fid,dims=dims,pos=0)rain_year += data_rain[i-12,*,*]ENVI_FILE_MNG,id=fid,/removeENDFORo_fn=DIALOG_PICKFILE(title='rain_year计算结果保存为')ENVI_WRITE_ENVI_FILE,rain_year,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=1,dsta_type=4,interleave=2,offset=offset,map_info=map_info;读取平均气温数据并计算每个月的平均气温data_temp = MAKE_ARRAY(12,data_size[1],data_size[2],type=data_size[3])temp_ave = MAKE_ARRAY(1,data_size[1],data_size[2],type=data_size[3])PRINT,SIZE(data_temp)temp_average=FLTARR(12)FOR i=24,36-1 DO BEGINENVI_OPEN_FILE,fns[i],r_fid=fidENVI_FILE_QUERY,fid,dims=dimsdata_temp[i-24,*,*]=ENVI_GET_DATA(fid=fid,dims=dims,pos=0)temp_ave += data_temp[i-24,*,*] ;每个像素一年的总温度r=REFORM(data_temp[i-24,*,*])r1=WHERE((r GT -50)AND(r LE 50))r2=r[r1]temp_average[i-24]=mean(r2,/double,/NAN)  ;每个月的平均温度ENVI_FILE_MNG,id=fid,/removeENDFORtemp_ave/=12temp_ave_6_10=mean(data_temp[5:9,*,*],dimension=1,/NAN) ;6至10月的平均气温o_fn=DIALOG_PICKFILE(title='temp_ave计算结果保存为')ENVI_WRITE_ENVI_FILE,temp_ave_6_10,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=1,dsta_type=4,interleave=2,offset=offset,map_info=map_info;读取NDVI平均值数据,并获取NDVI最大值月份work_dir1=DIALOG_PICKFILE(title='请选择NDVI文件路径',/directory);获取文件路径CD,work_dir1fns_NDVI = FILE_SEARCH(work_dir1,'*.dat',count = fnums_NDVI);PRINT,fns_NDVI,fnums_NDVIdata_NDVI=MAKE_ARRAY(12,data_size[1],data_size[2],type=data_size[3])FVC=MAKE_ARRAY(12,data_size[1],data_size[2],type=data_size[3])PRINT,SIZE(data_NDVI)NDVI_average=FLTARR(12)NDVI_min=FLTARR(12)NDVI_max=FLTARR(12)NDVI_ave_max =0.0FOR i=0,12-1 DO BEGINENVI_OPEN_FILE,fns_NDVI[i],r_fid=fidENVI_FILE_QUERY,fid,dims=dimsdata_NDVI[i,*,*]=ENVI_GET_DATA(fid=fid,dims=dims,pos=0);    print,data_NDVI[i,500,500]r=REFORM(REFORM(data_NDVI[i,*,*],data_size[1]*data_size[2],1))a=WHERE(FINITE(r))r1=r[a]b=WHERE((r1 LT 1) AND (r1 GT 0))r2=r1[b]NDVI_average[i]=mean(r2,/NAN)NDVI_min[i]=MIN(r2,/NAN)NDVI_max[i]=MAX(r2,/NAN);计算植被覆盖度FVC[i,*,*]=(data_NDVI[i,*,*]-NDVI_min[i])/(NDVI_max[i]-NDVI_min[i])IF NDVI_average[i] GT NDVI_ave_max THEN BEGINNDVI_ave_max = NDVI_average[i]NDVI_max_mouth = i+1ENDIFENVI_FILE_MNG,id=fid,/removeENDFORPRINT,'NDVI_max_mouth=',NDVI_max_mouth;NDVI和FVC最大值合成ANDVI=max(data_NDVI,dimension=1,/NAN)ANDVI=(ANDVI < 1.1)o_fn=DIALOG_PICKFILE(title='ANDVI计算结果保存为')ENVI_WRITE_ENVI_FILE,ANDVI,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=1,dsta_type=4,interleave=2,offset=offset,map_info=map_infoAFVC=max(FVC,dimension=1,/NAN)AFVC=(ANDVI < 1)o_fn=DIALOG_PICKFILE(title='AFVC计算结果保存为')ENVI_WRITE_ENVI_FILE,AFVC,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=1,dsta_type=4,interleave=2,offset=offset,map_info=map_info;释放空间r=!NULLr1=!NULLa=!NULLr2=!NULLb=!NULL;温度胁迫因子1;temp_stress1=cal_temp_stress1(data_temp[NDVI_max_mouth-1,*,*])temp_stress1= 0.8+0.02*data_temp[NDVI_max_mouth-1,*,*]-0.0005*data_temp[NDVI_max_mouth-1,*,*]^2o_fn=DIALOG_PICKFILE(title='温度胁迫因子1计算结果保存为')ENVI_WRITE_ENVI_FILE,temp_stress1,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=1,dsta_type=4,interleave=2,offset=offset,map_info=map_info;温度胁迫因子2temp_stress2 = cal_temp_stress2(data_temp,temp_average,NDVI_max_mouth)o_fn=DIALOG_PICKFILE(title='温度胁迫因子2计算结果保存为')ENVI_WRITE_ENVI_FILE,temp_stress2,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=12,dsta_type=4,interleave=2,offset=offset,map_info=map_info;计算水分胁迫因子W=cal_water_stress(data_temp,data_rain)o_fn=DIALOG_PICKFILE(title='水分胁迫银因子计算结果保存为')ENVI_WRITE_ENVI_FILE,W,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=12,dsta_type=4,interleave=2,offset=offset,map_info=map_info;计算FPARFPAR=cal_FPAR(data_NDVI,NDVI_min,NDVI_max)o_fn=DIALOG_PICKFILE(title='FPAR计算结果保存为')ENVI_WRITE_ENVI_FILE,FPAR,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=12,dsta_type=4,interleave=2,offset=offset,map_info=map_infodata_temp=!nulldata_rain=!null;创建数组计算NPPreal_use_rate=MAKE_ARRAY(12,data_size[1],data_size[2],type=data_size[3])APAR=MAKE_ARRAY(12,data_size[1],data_size[2],type=data_size[3])NPP=MAKE_ARRAY(12,data_size[1],data_size[2],type=data_size[3])ANPP=MAKE_ARRAY(1,data_size[1],data_size[2],type=data_size[3])FOR i=0,11 DO BEGIN;实际光能利用率max_use_rate=0.389;植被层对入射光合有效辐射吸收比例real_use_rate[i,*,*] = temp_stress1*temp_stress2[i,*,*]*W[i,*,*]*max_use_rate;光和有效辐射APAR[i,*,*]=data_radiate[i,*,*]*FPAR[i,*,*]*0.5    ;净初级生产力NPP[i,*,*]=APAR[i,*,*]*real_use_rate[i,*,*];全年累计净初级生产力IF NPP[i,500,500] GE 0 THEN BEGINANPP += NPP[i,*,*]ENDIFENDFOR;保存输出结果o_fn=DIALOG_PICKFILE(title='real_use_rate计算结果保存为')ENVI_WRITE_ENVI_FILE,real_use_rate,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=12,data_type=4,interleave=2,offset=offset,map_info=map_infoo_fn=DIALOG_PICKFILE(title='APAR计算结果保存为')ENVI_WRITE_ENVI_FILE, APAR,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=12,data_type=4,interleave=2,offset=offset,map_info=map_infoo_fn=DIALOG_PICKFILE(title='NPP计算结果保存为')ENVI_WRITE_ENVI_FILE,NPP,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=12,data_type=4,interleave=2,offset=offset,map_info=map_infoo_fn=DIALOG_PICKFILE(title='ANPP计算结果保存为')ENVI_WRITE_ENVI_FILE,ANPP,out_name=o_fn,/no_copy,ns=data_size[1],$nl=data_size[2],nb=1,data_type=4,interleave=2,offset=offset,map_info=map_infoEND;计算温度胁迫因子1
FUNCTION cal_temp_stress1,datatemp_stress1 = 0.8+0.02*data-0.0005*data^2w=WHERE(data LE -10)temp_stress1[w]=0.0RETURN,temp_stress1
END;计算温度胁迫因子2
FUNCTION cal_temp_stress2,data_temp,temp_average,NDVI_max_mouthdata_size=SIZE(data_temp)temp_stress2=MAKE_ARRAY(12,data_size[2],data_size[3],type=data_size[4]);计算最大气温月份的温度胁迫因子before = 1.1814/(1+EXP(0.2*(data_temp[NDVI_max_mouth-1,*,*]-10-data_temp[NDVI_max_mouth-1,*,*])))after = 1.0/(1+EXP(0.3*(-data_temp[NDVI_max_mouth-1,*,*]-10+data_temp[NDVI_max_mouth-1,*,*])))temp_stress2[NDVI_max_mouth-1,*,*] = before*afterFOR i=0,12-1 DO BEGINIF ((temp_average[i] GT temp_average[NDVI_max_mouth-1]+10)$OR(temp_average[i] LT temp_average[NDVI_max_mouth-1]-13)) THEN BEGINtemp_stress2[i-1,*,*]=0.5*  temp_stress2[NDVI_max_mouth-1,*,*]ENDIF ELSE BEGINbefore=1.1814/(1+EXP(data_temp[NDVI_max_mouth-1,*,*]-10-data_temp[i,*,*]))after=1.0/(1+EXP(0.3*(-data_temp[NDVI_max_mouth-1,*,*]-10-data_temp[i,*,*])))temp_stress2[i,*,*]=before+afterENDELSEENDFORbefore=!NULLafter=!NULLRETURN,temp_stress2
END;计算水分胁迫因子
FUNCTION cal_water_stress,data_temp,data_rain;十二个月总和热量指标I = total((data_temp/5)^1.514,1,/NAN)data_size = SIZE(I);因地而异的常数a=(0.6751*I^3-77.1*I^2+17920*I+482390)/1000000data_size=SIZE(data_temp);开辟空间Ep0=MAKE_ARRAY(12,data_size[2],data_size[3],type=data_size[4])Rn=MAKE_ARRAY(12,data_size[2],data_size[3],type=data_size[4])EET=MAKE_ARRAY(12,data_size[2],data_size[3],type=data_size[4])PET=MAKE_ARRAY(12,data_size[2],data_size[3],type=data_size[4])W=MAKE_ARRAY(12,data_size[2],data_size[3],type=data_size[4])FOR k=0,11 DO BEGIN;局地潜在蒸发量Ep0[k,*,*] = 16*(10*data_temp[k,*,*]/I)^a;像元此月份的地表净辐射量Rn[k,*,*] = ((Ep0[k,*,*]*data_rain[k,*,*])^0.5)*(0.369+0.598*(Ep0[k,*,*]/data_rain[k,*,*])^0.5);此月份的实际蒸散发量EET[k,*,*] = (data_rain[k,*,*] * Rn[k,*,*] * (data_rain[k,*,*]^2 + Rn[k,*,*]^2 + data_rain[k,*,*] * Rn[k,*,*]))/$((data_rain[k,*,*] + Rn[k,*,*]) * (data_rain[k,*,*]^2 + Rn[k,*,*]^2));此月份的潜在蒸散量PET[k,*,*] = (EET[k,*,*] + Ep0[k,*,*])/2;水分胁迫因子W[k,*,*] = 0.5+0.5*EET[k,*,*]/PET[k,*,*]ENDFORI=!NULLa=!NULLEp0=!NULLRn=!NULLEET=!NULLPET=!NULLRETURN,W
END;计算FPAR
FUNCTION cal_FPAR,data_NDVI,NDVI_min,NDVI_maxdata_size=SIZE(data_NDVI);开辟空间SR=MAKE_ARRAY(12,data_size[2],data_size[3],type=data_size[4])FPAR_NDVI=MAKE_ARRAY(12,data_size[2],data_size[3],type=data_size[4])FPAR_SR=MAKE_ARRAY(12,data_size[2],data_size[3],type=data_size[4])FPAR=MAKE_ARRAY(12,data_size[2],data_size[3],type=data_size[4])SR_min=FLTARR(12) & SR_max=FLTARR(12);参数设置FPAR_min=0.001 & FPAR_max=0.95;循环计算FOR i=0,11 DO BEGINSR[i,*,*]=(1+data_NDVI[i,*,*])/(1-data_NDVI[i,*,*])SR_min[i] = MIN(SR[i,*,*],/NAN)SR_max[i] = MAX(SR[i,*,*],/NAN)FPAR_NDVI[i,*,*]=(data_NDVI[i,*,*]-NDVI_min[i])*(FPAR_max-FPAR_min)$/(NDVI_max[i]-NDVI_min[i])+FPAR_minFPAR_SR[i,*,*]=(SR[i,*,*]-SR_min[i])*(FPAR_max-FPAR_min)$/(SR_max[i]-SR_min[i])+FPAR_minFPAR[i,*,*]=(FPAR_NDVI[i,*,*]+FPAR_SR[i,*,*])*0.5ENDFOR;保存文件o_fn=DIALOG_PICKFILE(title='FPAR_NDVI计算结果保存为')ENVI_WRITE_ENVI_FILE,FPAR_NDVI,out_name=o_fn,/no_copy,ns=data_size[2],$nl=data_size[3],nb=12,data_type=4,interleave=2,offset=offset,map_info=map_infoo_fn=DIALOG_PICKFILE(title='FPAR_SR计算结果保存为')ENVI_WRITE_ENVI_FILE, FPAR_SR,out_name=o_fn,/no_copy,ns=data_size[2],$nl=data_size[3],nb=12,data_type=4,interleave=2,offset=offset,map_info=map_info;释放空间SR=!NULLFPAR_NDVI=!NULLFPAR_SR=!NULLRETURN,FPAR
END

CASA方法估算NPP（IDL+ENVI）相关推荐

Python实现CASA模型估算植被净初级生产力——最适温度篇
在估算植被的NPP时,相信很多人都会选择CASA模型,因为它实在是一个参数少易实现且应用广泛的光能利用率模型啊. CASA模型估算NPP,公式就是NPP = APAR * E,APAR是植被吸收的光合 ...
【300】◀▶ IDL - ENVI API
参考:ENVI API 参考:ENVI Classic Display 序号类名称功能说明语法 & 举例 01 ENVI 函数 ====<<<< De ...
IDL(ENVI/IDL) 简（jian）明（lou）教程：二、ENVI/IDL批处理入门（以投影转换为例）
二.ENVI/IDL批处理入门 ENVI/IDL集成了ENVI软件的高级功能,比如打开文件直接使用envi_open_file, File ,r_fid=fid即可,不用考虑什么格式等,再比如做文件投 ...
08.用户故事与敏捷方法——估算用户故事笔记
00.估算故事最好方法: *无论什么时候获得有关故事的新信息,都允许我们改变之前的想法 *适用于史诗故事和小故事 *不需要花很多时间 *提供进度和剩余工作的有用信息 *不太精确的估算也不会有太大问题 ...
怀旧服服务器在线人数查询,民间大神用土方法估算出魔兽怀旧服在线人数震动了官方...
随着<魔兽怀旧服>的正式发布,许多老玩家都在第一时间聚集到当初所怀念的"世界"当中.即便官方在开服前紧急新增了几组服务器,也难以抵挡住玩家所带来的热情,排队是正常的事情 ...
如何用数学方法估算一个女生前男友的数量？
全世界只有3.14 % 的人关注了爆炸吧知识如果一个女生说,她集齐了十二个星座的前男友,我们应该如何估计她前男友的数量? 小学生:这个问题相当简单公式: 数据: {白羊座, 金牛座, 双子座, ...
CASA——估算陆地生态系统植被净初级生产力（NPP）的经典模型
文章目录第一讲:CASA模型介绍第二讲:CASA初步操作第三讲:CASA数据制备(一) 第四讲:CASA数据制备(二) 第五讲:CASA数据制备(三) 第六讲:CASA数据制备(四) 第七讲:土 ...
基于改进的 CASA 模型反演30m分辨率NPP
应太多粉丝要求,想学习下npp计算怎么上手,本文从数据获取.处理.以及计算这三个方面来展开教学,希望大家能够点个关注点个赞加收藏! NPP软件下载百度云链接: 基于改进的 CASA 模型反演 NP ...
IDL学习记录和Java调用IDL方法
IDL学习记录和Java调用IDL方法 2018年02月06日 08:32:02 回首1949 阅读数:385更多个人分类: 随想版权声明:乐呵乐呵得了 https://blog.csdn.net ...
ENVI 丢失idl.dll（基于win10 安全中心的错误识别检测）
最近使用ENVI 5.3 64bit ,打开时候弹出对话框丢失idl.dll, 提示重新安装可以解决此问题.由于重新安装过于复杂,所以这里找了一下丢失原因:原来是由于win10 安全中心将idl.d ...

CASA方法估算NPP（IDL+ENVI）

CASA方法估算NPP（IDL+ENVI）

一、NPP的定义

二、估算流程

三、实现代码（IDL)

CASA方法估算NPP（IDL+ENVI）相关推荐

最新文章

热门文章

CASA方法估算NPP（IDL+ENVI）

CASA方法估算NPP（IDL+ENVI）

一、NPP的定义

二、估算流程

三、实现代码 （IDL)

CASA方法估算NPP（IDL+ENVI）相关推荐

最新文章

热门文章

三、实现代码（IDL)