CASA方法估算NPP（IDL+ENVI）

最新推荐文章于 2025-02-12 17:12:19 发布

南鱼木舟

最新推荐文章于 2025-02-12 17:12:19 发布

阅读量2.8w

点赞数 57

分类专栏： IDL 文章标签： idl 算法

本文链接：https://blog.csdn.net/qq_42315165/article/details/111240678

版权

IDL 专栏收录该内容

1 篇文章

订阅专栏

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ε2(x,y,t)=1.18141/(1+e^(0.2*[Topt(x,y)-10-T(x,y,t)]) )*1/(1+e^(0.3×[-Topt(x,y)-10+T(x,y,t)] ) )*** (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_NPP

  e=ENVI()
  ;读取数据
  work_dir=DIALOG_PICKFILE(title='请选择文件路径',/directory);获取文件路径
  CD,work_dir
  fns = 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 BEGIN
    ENVI_OPEN_FILE,fns[i],r_fid=fid
    ENVI_FILE_QUERY,fid,dims=dims,interleave=interleave,offset=offset
    map_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,/remove
  ENDFOR
  data_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 BEGIN
    ENVI_OPEN_FILE,fns[i],r_fid=fid
    ENVI_FILE_QUERY,fid,dims=dims
    data_rain[i-12,*,*] = ENVI_GET_DATA(fid=fid,dims=dims,pos=0)
    rain_year += data_rain[i-12,*,*]
    ENVI_FILE_MNG,id=fid,/remove
  ENDFOR
  o_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 BEGIN
    ENVI_OPEN_FILE,fns[i],r_fid=fid
    ENVI_FILE_QUERY,fid,dims=dims
    data_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,/remove
  ENDFOR
  temp_ave/=12
  temp_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_dir1
  fns_NDVI = FILE_SEARCH(work_dir1,'*.dat',count = fnums_NDVI)
  ;PRINT,fns_NDVI,fnums_NDVI

  data_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.0
  FOR i=0,12-1 DO BEGIN
    ENVI_OPEN_FILE,fns_NDVI[i],r_fid=fid
    ENVI_FILE_QUERY,fid,dims=dims
    data_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 BEGIN
      NDVI_ave_max = NDVI_average[i]
      NDVI_max_mouth = i+1
    ENDIF
    ENVI_FILE_MNG,id=fid,/remove
  ENDFOR
  PRINT,'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_info
    
  AFVC=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=!NULL
  r1=!NULL
  a=!NULL
  r2=!NULL
  b=!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,*,*]^2
  o_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

  ;温度胁迫因子2
  temp_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

  ;计算FPAR
  FPAR=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_info

  data_temp=!null
  data_rain=!null
  
  ;创建数组计算NPP
  real_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 BEGIN
      ANPP += NPP[i,*,*]
    ENDIF
    
  ENDFOR

  ;保存输出结果
  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_info
    
  o_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_info
    
  o_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_info
    
  o_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_info
 
END

;计算温度胁迫因子1
FUNCTION cal_temp_stress1,data
  temp_stress1 = 0.8+0.02*data-0.0005*data^2
  w=WHERE(data LE -10)
  temp_stress1[w]=0.0
  RETURN,temp_stress1
END

;计算温度胁迫因子2
FUNCTION cal_temp_stress2,data_temp,temp_average,NDVI_max_mouth

  data_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*after

  FOR i=0,12-1 DO BEGIN
    IF ((temp_average[i] GT temp_average[NDVI_max_mouth-1]+10)$
      OR(temp_average[i] LT temp_average[NDVI_max_mouth-1]-13)) THEN BEGIN
      temp_stress2[i-1,*,*]=0.5*  temp_stress2[NDVI_max_mouth-1,*,*]
    ENDIF ELSE BEGIN
      before=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+after
    ENDELSE
  ENDFOR

  before=!NULL
  after=!NULL
  RETURN,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)/1000000
  data_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,*,*]
  ENDFOR

  I=!NULL
  a=!NULL
  Ep0=!NULL
  Rn=!NULL
  EET=!NULL
  PET=!NULL
  RETURN,W
END

;计算FPAR
FUNCTION cal_FPAR,data_NDVI,NDVI_min,NDVI_max

  data_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 BEGIN
    SR[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_min
    FPAR_SR[i,*,*]=(SR[i,*,*]-SR_min[i])*(FPAR_max-FPAR_min)$
      /(SR_max[i]-SR_min[i])+FPAR_min
    FPAR[i,*,*]=(FPAR_NDVI[i,*,*]+FPAR_SR[i,*,*])*0.5
  ENDFOR
    ;保存文件
  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_info

  o_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=!NULL
  FPAR_NDVI=!NULL
  FPAR_SR=!NULL

  RETURN,FPAR
END