首页 > 分享 > 利用IDL计算植被覆盖度(VFC)

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

花匠小妙招
2024-12-26 18:10

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 123456

例如

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

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

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

1.2 过程

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

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

例如

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

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

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

2. 正题

思路是，计算利用近红波段和红波段计算出 N D V I NDVI NDVI ，再根据5%和95%的置信区间统计出 N D V I s o i l NDVI_{soil} NDVIsoil 和 N D V I v e g NDVI_{veg} NDVIveg (统计的时候需要去除水体，至于怎么去水，我在我的SRT中用的是MOD44W 250m全球水体掩膜的数据集)，在这里我简单的认为水体是 N D V I < 0 NDVI<0 NDVI<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

1234567891011121314151617181920212223242526272829303132333435

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

2.2 第二步

计算NDVI
N D V I = N I R − R e d N I R + R e d NDVI = dfrac{NIR-Red}{NIR+Red} NDVI=NIR+RedNIR−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 1234567

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

2.3第三步

根据直方图统计 N D V I s o i l NDVI_{soil} NDVIsoil， N D V I v e g NDVI_{veg} NDVIveg，也有人称为 N D V I m a x NDVI_{max} NDVImax， N D V I m i n NDVI_{min} NDVImin, 我习惯前者,于是第三个function。
输入是 N D V I NDVI NDVI，输出是 N D V I s o i l NDVI _{soil} NDVIsoil和 N D V I v e g NDVI_{veg} NDVIveg

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

1234567891011121314151617

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

2.4 第四步

计算VFC,VFC的计算公式
V F C = N D V I − N D V I s o i l N D V I V e g − N D V I s o i l VFC = frac{NDVI- NDVI_{soil}}{NDVI_{Veg}-NDVI_{soil}} VFC=NDVIVeg−NDVIsoilNDVI−NDVIsoil当然会有计算出 V F C < 0 VFC<0 VFC<0 或者 V F C > 1 VFC>1 VFC>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 123456789101112131415

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