像元二分模型计算植被覆盖度

最新推荐文章于 2024-09-13 10:05:52 发布

江北20190411

最新推荐文章于 2024-09-13 10:05:52 发布

阅读量1.1w

点赞数 8

分类专栏：遥感 IDL

本文链接：https://blog.csdn.net/tk20190411/article/details/105861955

版权

遥感同时被 2 个专栏收录

40 篇文章 41 订阅

订阅专栏

IDL

15 篇文章 9 订阅

订阅专栏

前言

关于改进的像元二分模型第一次出现是在李苗苗老师的博士论文中，但她后来发表的一篇期刊文章把这一部分单独摘了出来并做了应用案例，期刊链接在下方。

密云水库上游植被覆盖度的遥感估算

论文要点

计算公式为：Fc =(NDVI -NDVIsoil) (NDVIveg -NDVIsoil)
问题的关键在于NDVIveg 与NDVIsoil如何取值。

其中 ,NDVIsoil 为完全是裸土或无植被覆盖区域的NDVI 值 , NDVIveg 则代表完全被植被所覆盖的像元的NDVI 值, 即纯植被像元的NDVI 值。NDVIveg 代表着全植被覆盖像元的最大值, 由于植被类型的影响, NDVIveg 值也会随着时间和空间而改变。因此 , 计算植被覆盖度时, 既使同一景影像, 对于NDVIsoil 和NDVIveg 值不能取固定值。在实地应用中,植被覆盖类型随土地利用类型而变化, 对于某一土地利用类型, 由于植被类型近似, 其 NDVIveg 值也相近;而对于特定的土壤类型 ,其 NDVIsoil 值也应该是一定的。因此, 土地利用图和土壤图应作为计算 NDVIveg 和NDVIsoil 值的基础

计算过程

辅助数据介绍

全球土壤数据
HWSD全球土壤数据下载处理
全球土地分类数据
MCD12Q1
关于MCDQ1，由官网介绍可知数据中包含5种分类方法得到的结果。
在这里插入图片描述
数据的说明手册目前网址404了。但我们可以查看数据的头文件获取分类代码信息。
这里贴出了type3的分类代码信息，比较符合我们计算植被覆盖度的目的。

Dataset #3: Land_Cover_Type_3  
  Dims: BYTE (2400 x 2400)    
Attribute 3-1: "long_name"  
  "Land_Cover_Type_3"   
Attribute 3-2: "Water Bodies"  
  0   
Attribute 3-3: "Grasslands"  
  1   
Attribute 3-4: "Shrublands"  
  2   
Attribute 3-5: "Broadleaf Croplands"  
  3   
Attribute 3-6: "Savannas"  
  4   
Attribute 3-7: "Evergreen Broadleaf Forests"  
  5   
Attribute 3-8: "Deciduous Broadleaf Forests"  
  6   
Attribute 3-9: "Evergreen Needleleaf Forests"  
  7  
Attribute 3-10: "Deciduous Needleleaf Forests"  
  8   
Attribute 3-11: "Unvegetated"  
  9   
Attribute 3-12: "Urban and Built-up Lands"  
  10   
Attribute 3-13: "Unclassified"  
  255   
Attribute 3-14: "valid_range"  
  0, 10    
Attribute 3-15: "_FillValue"  
  255

使用MCTK对type3进行投影，重采样，至此LULC和HWSD，NDVI数据都准备好了。

IDL实现

pro vfc_cal
  compile_opt idl2
  e = envi()
  workpath = 'F:\FY3D_MERSI_data\Subset'
  cd,workpath
  ndvifile = workpath + '\NDVI'
  lulcfile = workpath + '\LULC'
  hwsdfile = workpath + '\HWSD'
  
  ndvir = e.openraster(ndvifile)
  hwsdr = e.openraster(hwsdfile)
  lulcr = e.openraster(lulcfile)
  
  spatialref = ndvir.spatialRef
  
  ndvi = ndvir.getdata(bands=0)
  hwsd = hwsdr.getdata(bands=0)
  lulc = lulcr.getdata(bands=0)
  
  ndvi[where(finite(ndvi,/nan))]  = mean(ndvi,/nan)
  lulc[where(finite(lulc,/nan))]  = 4
  
  hwsd_class = hwsd[uniq(hwsd,sort(hwsd))]
  lulc_class = lulc[uniq(lulc,sort(lulc))]
  print,hwsd_class,lulc_class
   
  lutmin = make_array(2,n_elements(hwsd_class),type=5)
  for i=0,n_elements(hwsd_class)-1 do begin 
    lutmin[0,i] = hwsd_class[i]
    ndvi_hwsd = ndvi[where(hwsd eq lutmin[0,i])] 
    ndvi_hwsd = ndvi_hwsd[sort(ndvi_hwsd)]
    ndvi_hwsd = ndvi_hwsd[where(hwsd ge -0.1)]
    number = n_elements(ndvi_hwsd)
    lutmin[1,i] =  ndvi_hwsd[round(0.05*number)]  
  endfor

  lutmax = make_array(2,n_elements(lulc_class),type=5)
  for i=0,n_elements(lulc_class)-1 do begin
    lutmax[0,i] = lulc_class[i]
    ndvi_lulc = ndvi[where(lulc eq lutmax[0,i])]
    ndvi_lulc = ndvi_lulc[sort(ndvi_lulc)]
    number = n_elements(ndvi_lulc)
    lutmax[1,i] =  ndvi_lulc[floor(0.95*number)]
  endfor
 
  dim = size(ndvi)
  n = n_elements(ndvi)
  vfc = make_array(dim[1],dim[2],type=5)
  for i=0,n-1 do begin
    for j=0,n_elements(hwsd_class)-1 do begin
      for k=0,n_elements(lulc_class)-1 do begin
        if lulc[i] eq lutmax[0,k] and hwsd[i] eq lutmin[0,j] then begin
          case 1 of 
            ndvi[i] ge lutmax[1,k] : vfc[i] = 1
            ndvi[i] le lutmin[1,j] : vfc[i] = 0
            else: vfc[i] = (ndvi[i]-lutmin[1,j]) / (lutmax[1,k]-lutmin[1,j])
          endcase
        endif
      endfor
    endfor
  endfor
  
  tempfn = e.GetTemporaryFilename()
  vfcr = e.CreateRaster(tempfn, vfc, SpatialRef = spatialref)
  vfcr.save
 
end