利用IDL计算植被覆盖度(VFC)

最新推荐文章于 2024-05-02 00:23:45 发布

N某人

最新推荐文章于 2024-05-02 00:23:45 发布

阅读量8.2k

点赞数 27

分类专栏： IDL 文章标签： IDL VFC

本文链接：https://blog.csdn.net/weixin_43342986/article/details/108156509

版权

IDL 专栏收录该内容

8 篇文章 16 订阅

订阅专栏

0. 前言

正巧IDL实验课考核的作业是利用4个Function和主Pro过程写一个遥感图像处理的代码，要求是前一个方法的输出是另一个方法的输入。以前一直想着能不能计算NDVI和植被覆盖度(VFC)用IDL写出来，因为几乎所有老师上课一开始都是推荐用IDL写遥感图像处理的程序的。
先区分一些概念，ENVI、IDL、C++(cpp)。
ENVI(The Environment for Visualizing Images)是我们遥感专业处理中常用的一款可视化的软件，原来是RSI公司出品。现在变 Exelis Visual Information Solutions公司出品的。
而可视化软件肯定是用编程语言写出来的，而 ENVI 就是用IDL语言写出来的。
IDL是一种编程语言，它和MATLAB、Python、java等语言一样是一种解释性的语言，你在一定的编辑器甚至在 command 中输入一些基础语句，就能得到一定的运行结果。而像 c++、c 是需要进行编译之后在进行运行的，这和前面的几种语言有本质的不同。
而IDL的底层是由 c++等语言写出来的。

解释清楚这些再谈谈我自己对这个 project 的理解,即利用 MODIS 数据进行某地的植被覆盖度的变化监测

数据下载在我以前的 blog 中已经都有记录了，而实际上这个项目的真正需求和难点是长时间序列的(本人也做过，不过处理过程是用ENVI做的)，其实长时间序列只是用 MRT(MODIS Reprejection Tool) 对数据进行了重投影，再进行一个最大值滤波（最大值滤波就是把所有图像读进去用一个 '>'比较一下罢了），之后再是一些基本操作，但是我觉得这个（科研的基础入门训练）而不是科研，仅仅是掌握最基本的思维能力、文献查询、还有ENVI的最基本处理罢了，去发 paper 的话还是没有什么价值（这个项目已经被做烂了，而且GEE(Goole Earth Engine)的兴起，正在逐渐淘汰这些本地软件的处理，现在更注重遥感大数据云端处理，附以其他数据进行联合分析）这些联合分析其实就是遥感分析中的地理相关分析法，是什么导致了植被覆盖度这样变化，是城市扩张？还是自然的变换？而如果是城市扩张我又需要用什么数据去进行佐证城市在扩张（比如夜光数据，土地覆盖利用数据....），但是需要注意的是，我们在选择佐证因素的时候，需要先确定好有哪些现存的数据可以提供佐证，不可能用获得不了的卫星数据来佐证,亦或是时间序列无法满足你自身证明要求的数据。

1. 语法

1.1 函数、功能、方法

function name,var1,var2...
	process....
	balalala....
	
return var
end

例如

function pass_var,var1,var2
    c = [var1,var2]
    return,c
end

这个pass_var为函数名, $var_1$ ， $var_2$ 为输入变量， $c$ 为输出变量，需要注意的是return后面用的是逗号，并且返回值只能有一个，但是这个值可以是矩阵、单一值等等。

------------------------------------分割线------------------------------

1.2 过程

过程又称为procedure，在IDL的语法是

pro proname,var1,var2,....	
	process,...
end

例如

pro test6
    c = pass_var(1,2)
    print,c
end

这个pro就调用了上面的function，成就了一种联动。

------------------------------------分割线------------------------------

2. 正题

思路是，计算利用近红波段和红波段计算出 $N D V I$ ，再根据5%和95%的置信区间统计出 $NDVI_{soil}$ 和 $NDVI_{veg}$ (统计的时候需要去除水体，至于怎么去水，我在我的SRT中用的是MOD44W 250m全球水体掩膜的数据集)，在这里我简单的认为水体是 $N D V I < 0$ 的部分

输入：近红波段图像和红波段图像。
输出：VFC植被覆盖度

需要注意的是，我在前期已经将图像的背景赋值为 NAN 了，且我用的输入图像是MOD13Q1经过MRT处理输出的近红和红波段。

------------------------------------分割线------------------------------

2.1 第一步

选取图像，这里我用的dialog_pickfile.可以自己选取图像，这里是主过程。

pro caculate_VFC
   ;coding = GBK
    ;version 1.0
    ;Calculate VFC with red band and near red band
    ;整个需要选择的图像主要为红波段图像、近红波段图像，计算NDVI，并根据NDVI计算直方图得到NDVI_min,NDVI_max
    ;Neverland,Aug 21st 2020
    ;测试中发现在计算VFC计算步骤中会导致浮点数非法操作，但是必须使用这种方法，否则会导致整幅图变为2值图像 
    filters = ['*.jpg;*.jpeg', '*.tif;*.tiff', '*.png'];      ;过滤器
    result1 = dialog_message('请选择红波段的图像',/information)
    red_file = dialog_pickfile(/READ, FILTER = filters,GET_PATH=work_dir,title='请选择红波段图像')       ;选择红波段图像
    result2 = dialog_message('请选择近红波段的图像',/information)
    nir_file = dialog_pickfile(/READ, FILTER = filters,path=work_dir,title='请选择近红波段图像')       ;选择近红波段图像
    b1 = read_tiff(red_file,Geotiff=GeoKeys)                  ;读取红波段
    b2 = read_tiff(nir_file,Geotiff=GeoKeys)                  ;读取近红波段
    mn_size = size(b1,/dimensions)                            ;图像的大小
    m = mn_size[0]                                            ;行
    n = mn_size[1]
    ；-------------------计算NDVI-----------------------
    NDVI = cal_NDVI(b1,b2)
    
    ;--------------------统计直方图----------------------
    NDVI_stastic = cal_maxmin(NDVI)
    NDVI_soil = NDVI_stastic[0]                ;5%的置信区间
    NDVI_veg = NDVI_stastic[1]                 ;95的置信区间
    print,'NDVI_soil: ',NDVI_soil, '; NDVI_veg:', NDVI_veg

    ----------------------计算VFC------------------------
    ;第四步，计算VFC并写出
    VFC = cal_vegfraction(NDVI,NDVI_soil,NDVI_veg)

    VFC_out = dialog_pickfile(/write,path=work_dir,title='请选择VFC保存位置的图像')          ;自定义写出的名称
    VFC_outname = VFC_out + '.tif'
    im1 = image(VFC,dimensions=[n,m],margin=0,window_title='VFC图像',order=1)        ;图像
    write_tiff,VFC_outname,VFC,/float,geotiff=Geokeys             
end

dialog_message弹出对话框，类型有warning，information…等等，然后dialog_pickfile交互选择图像，根据提示选就行。这个也是主 pro ，里面会调用之后写的函数。

2.2 第二步

计算NDVI
$\dfrac{NIR-Red}{NIR+Red}$
这应该是每个RSer必知的公式！

function cal_NDVI,b1,b2
    ;输入参数为b1红波段,b2近红波段并计算NDVI,最后写出
    NDVI = (float(b2)-b1)/(float(b2)+b1)
    NDVI_out = dialog_pickfile(/write,path=work_dir,title='请选择NDVI输出的图像路径与名称')
    write_tiff,NDVI_outname,NDVI,/float,geotiff=Geokeys
    return,NDVI
end

于是在主过程中补上去，就是代码中分割线下面的部分。

2.3第三步

根据直方图统计 $NDVI_{soil}$ ， $NDVI_{veg}$ ，也有人称为 $NDVI_{max}$ ， $NDVI_{min}$ , 我习惯前者,于是第三个function。
输入是 $N D V I$ ，输出是 $NDVI _{soil}$ 和 $NDVI_{veg}$

function cal_maxmin,NDVI

;输入为NDVI,根据NDVI统计5%和95%的置信区间

 w = where(finite(NDVI) gt 0,count)                 ;由于背景是NAN，如果不适用finite函数会导致NAN无法与float比较
 ht = histogram(NDVI[w],nbins=200,locations=locations)
    ht_acc = total(ht,/cumulative)/count
    
;统计累积频率为5%和95%的NDVI作为NDVI_soil和NDVI_veg
    w1 = where(ht_acc gt 0.05)                 ;5%的累计区间
    NDVI_soil = locations[w1[0]-1]
    
    w2 = where(ht_acc ge 0.95)                 ;95%的累计区间
    NDVI_veg = locations(w2[0])
    NDVI_stastic = [NDVI_soil,NDVI_veg]          ;合并以便写出
    return,NDVI_stastic
end

再补这个function进主pro中，这个function返回值有2个，所以返回了第0个和1个值，作为 $NDVI_{soil}$ 和 $NDVI_{veg}$ 。histogram是直方图，然后统计什么时候到了5%什么时候到了95%，就是 $NDVI_{\min}$ , $NDVI_{\max}$ 。每个地方有不同的植被特性，以及地区环境不同需要选取合适的置信区间。

2.4 第四步

计算VFC,VFC的计算公式
$\frac{NDVI- NDVI_{soil}}{NDVI_{Veg}-NDVI_{soil}}$ 当然会有计算出 $V F C < 0$ 或者 $V F C > 1$ 的情形，需要自己进行赋值。

;-----------------------------------------------计算植被覆盖度-----------------------------------------------------------
function cal_vegfraction,NDVI,NDVI_soil,NDVI_veg
    ;NDVI为归一化植被指数，NDVI_soil为5%裸土的NDVI值， 
     NDVI_veg为95%纯植被像元的NDVI值
      
     ;像元二分法
     VFC = (NDVI-NDVI_soil)/(NDVI_veg-NDVI_soil)
     
     ;将小于NDVI_soil的值变为0
     w = where(NDVI le NDVI_soil or finite(NDVI) le 0)
     VFC[w] = 0
      
      ;将NDVI大于NDVI_veg的像元变为1
      w = where(NDVI ge NDVI_veg)
      VFC[w] = 1

核心的也就一句话，之后都是异常值处理，where就是寻访出下标
然后再处理，最后再把这个function写入主程序之中。于是就写完了，其实后面还写了个密度分割的，但是想想太麻烦了就不贴了。
对了如果出现浮点数的警告也没关系，我当时也是警告，但照样能读和写计算完VFC的图像，如果真的要处理背景值，就需要看我后面总结的NAN处理的办法，把背景值赋值成空值了，当然用ENVI或者用GIS都是可以的，但是可视化都是“流氓”🐦🐦🐦。

4. 总结

以前一直觉得用IDL做这个Project很困难，后来今天花了2-3个小时写完了，发现也挺简单，其实多尝试就行，自己一直顺着思路写Debug就行了。

我也知道有些人觉得代码过长，看到一堆英文就不看了，从来无所谓，写blog是记录自己思考过程和收获的过程，你看不看又跟我没关系，我记录了，那我就收获了，我就成功了。

一世随缘，随缘一世才能活得自在。

最后，我将在另外一个账号写博客，有兴趣的可以关注。

N某人

关注

27
点赞
踩
139

收藏

觉得还不错? 一键收藏
10
评论
利用IDL计算植被覆盖度(VFC)

前言正巧IDL实验课考核的作业是利用4个Function和主Pro过程写一个遥感图像处理的代码，以前一直想着能不能计算NDVI和植被覆盖度(VFC)用IDL写出来，因为老是总强调处理遥感影像比较正统，一定程度上是这样说，因为IDL是Harris公司的，然后envi的底层是IDL，IDL的底层是cpp，写完了整个代码大概是140行，当然这个是仅仅利用2幅影像计算的，其实看大多数人的需求是长时间序列的，其实长时间序列只是进行了一个最大值滤波，之后再是几幅图像的计算。语法function name,..
复制链接

扫一扫