基于 IDL 的遥感要素序列提取代码

好文速递：基于 IDL 的遥感要素序列提取代码

 

    

  遥感要素序列是遥感应用中经常需要完成的工作，手动提取效率低下且不能满足大数据处理的要求。为了解决这些问题，我们使用 IDL 编写了自动批量提取遥感要素的代码。验证结果表明，我们基于 IDL编写的代码可有效提取遥感要素序列，且可以适应大数据量处理。

 

1 数据源

 

使用 MODIS 的 MOD13Q1 产品作为数据源。首先从 NASA 官网（https://ladsweb.nascom.nasa.gov/）下载 2000 至 2018 年覆盖西藏自治区范围的所有 MOD13Q1 数据，它是采用 hdf 数据格式存储的。然后使用 MRT（MODIS Reprojection Tool）工具和 ENVI5.1 的“小熊工具箱”插件完成图像的拼接和提取 NDVI 图层工作，并且存储为 tif 格式（见下图） 

                       拼接完成的 NDVI 遥感影像  

最后将每年的 NDVI 图像文件放在同一个文件夹里作为待处理数据。

                         分组存储的 NDVI 遥感影像 

2  代码实现 

 

2.1 遥感要素序列提取的原理 

 

        遥感要素一般存储在遥感图像中，我们往往需要连续提取多幅图像对应位置的图像信息。因此，遥感要素提取的一种思路就是将需要提取要素的位置转化为掩膜图像，给图像中需要提取要素的位置特殊的编码。然后就可以通过对比栅格编号提取需要的值并对其进行运算。由于遥感图像数据量巨大，必须使用循环来完成多次的数据提取，如需要提取要素的位置较多，还可以加入另一循环来完成。

 

2.2 程序框架 

 

       这个程序由运行环境设置模块、搜索 MODIS图像基本信息模块、掩膜数据读取模块、MODIS 图像读取模块、NDVI 序列提取模块和数据输出模块六个模块组成。在进行正式的数据提取之前需要配置 IDL 环境，使得它可以调用 ENVI5.1 的内置函数。环境配置完毕后，需要先搜索存储遥感图像的文件夹下遥感图像的数量和各自的文件名。搜索完毕后就可以开始 NDVI 数据提取的循环了。提取循环由内中外三层循环组成，外层循环负责遍历存储MODIS 图像的文件夹，中层循环逐个提取位置掩膜，内层循环逐个提取 MODIS 遥感图像特定位置的取值。

 

                                           代码流程图

 

2.3 运行环境设置模块 

 

IDL 代码编译有时候会出现所调用的 ENVI 函数找不到的情况，所以我们使用“COMPILE_OPT idl2”语句避免这种情况。要运行 IDL 的批处理程序还需要加载 ENVI 的部分核心 sav 文件，完成这个功能的代码是“ENVI,  /RESTORE_BASE_ SAVE_FILES”。IDL 的批处理模式需要一段代码开启，这个代码是“ENVI_BATCH_INIT, LOG_FILE = 'batch_log.txt'”。完成了以上步骤也就完成了 IDL 批处理模式的环境配置。

 

2.4 搜索 MODIS 图像基本信息

 

要想处理文件夹下的 MODIS 遥感影像，首先要搜索该文件夹包含的 MODIS 影像文件的个数，这个功能使用 file_search 函数完成。除此之外，还需要从遥感影像的文件名中读取它的时间信息，这个功能靠 strmid 函数从文件名字符串中截取完成。

 

 

2.5 图像读取 

 

要想提取图像中信息就得先打开图像，决定NDVI 提取范围的掩膜数据本质上也是图像，所以自动读取图像就成为这个程序中最为重要的部分。IDL 的图像读取过方法有很多，我们这里采用比较常规的方案。首先使用 ENVI_OPEN_FILE 函数返回文件名对应文件的 fid 参数，这个参数相当图像文件的编号。然后使用 ENVI_OPEN_FILE 根据返回的 fid 读取图像的波段数、行列数和文件类型等信息,最后使用 ENVI_GET_DATA 函数结合波段数等信息将 MODIS 图像读入为一个矩阵。这样就完成了图像的读取。

 

2.6 遥感要素提取核心算法 

 

 

遥感要素（NDVI）提取的核心算法本质上是比对遥感图像矩阵和掩膜矩阵，将特定行列的遥感图像矩阵记录在输出矩阵中。首先要设置一个足够大提取矩阵 a1。然后读入遥感图像矩阵和掩膜矩阵，遥感图像矩阵是一个每个单元格都记录着对应位置 10000*NDVI 值的矩阵，掩膜矩阵是一个和遥感图像矩阵具有相同行列数的矩阵，它的单元格里只包含 0、1 值，需要提取数据的位置取值为 1，其他位置取值为 0。然后逐栅格对比掩膜矩阵的取值，如果单元格取值为 1，就将遥感图像矩阵对应位置单元格的取值传给 a1 矩阵。最后对 a1 矩阵进行取平均值、最大值和最小值并传输给中间输出结果矩阵 dat_temp 就完成了一次 NDVI 的提取。

 

2.7 文件输出 

 

考虑到数据提取成果需要便于被多种软件读取，我们采用 dat 格式输出成果。在输出成果前还需要设置输出文件名，我们采用“掩膜名称+'static'+数据年份”的形式命名每个输出文件。数据输出采用 openw 函数完成。

 

2.8 变量说明

 

以下是源代码中使用的主要变量的变量说明: 

Root_path 为存放数据文件的目录； 

Mask_file 为存放掩膜数据的文件目录； 

outsum_path 为结果输出路径； 

outsum 为提取结果的输出文件名； 

data_tif 为掩膜图像矩阵； 

data_ndvi1 为遥感影像矩阵； 

a1 为单幅图像提取结果矩阵； 

data_temp 为输出结果矩阵。 

 

2.9 源代码

 

pro ROI_MODIS_NOVAA 

   COMPILE_OPT idl2 

   ENVI, /RESTORE_BASE_SAVE_FILES 

   ENVI_BATCH_INIT, LOG_FILE = 

'batch_log.txt' 

   ;配置运行环境 

   Data_tpe=’*.tif’ 

   ;全局变量，代表输入图像的格式 

   Root_path='H:\PROJECT\yajiang\MOD13Q1\' 

   ;存放数据文件的目录 

   Mask_file='H:\me\个人使用材料\雅江物候论文

\xinchuli\' 

   ;存放掩膜数据的文件目录 

   outsum_path="H:\me\个人使用材料\雅江物候论

文\NDVI_serisws\" 

   ;结果输出路径 

    for cc=2005,2018 do begin 

    data_temp=fltarr(5,50) 

  i=fix(50) 

    base_path='NDVI_'+strmid(string(cc),8,4) 

    search_path= Root_path +base_path 

    filename=file_search(search_path,'*.tif') 

    filename1=file_basename(filename) 

    filename2= strmid(filename1,9,4) 

    filename3=strmid(filename1,13,3) 

    ;读取MODIS数据的基本信息 

    mask_path=file_search(Mask_file,'*.tif') 

    ;搜索掩膜数据 

    n = N_ELEMENTS(filename) 

    b = N_ELEMENTS(mask_path) 

    for j=0,b-1 do begin 

ENVI_OPEN_FILE,mask_path[j],r_fid=tempfid 

ENVI_OPEN_FILE,temp-

fid,NB=NB,DIMS=DIMS,BNAMES=BNAMES,

DATA_TYPE=DATA_TYPE,NL=NL,NS=NS 

       data_tif = ENVI_GET_DATA(fid=tempfid, 

dims=dims, pos=0) 

       c=file_basename(mask_path[j]) 

       ;读取掩膜数据 

       FOR i=0,n-1 DO BEGIN 

ENVI_OPEN_FILE,filename[i],r_fid=tempfid 

ENVI_FILE_QUERY,temp-

fid,NB=NB,DIMS=DIMS,BNAMES=BNA

MES,DATA_TYPE=DATA_TYPE,NL=NL

,NS=NS 

       data_ndvi1 = ENVI_GET_DATA(fid=temp-

fid, dims=dims, pos=0) 

       index1=where((data_tif eq 1) and (data_ndvi1 

gt -10000)) 

       ;读取MODIS数据 

a1=data_ndvi1[index1] 

 data_temp[0,i]=filename2[i] 

       data_temp[1,i]=filename3[i] 

       data_temp[2,i]=mean(a1) 

       data_temp[3,i]=max(a1) 

       data_temp[4,i]=min(a1) 

       ;提取淹没范围内NDVI的均值 

       ENDFOR 

outsum=outsum_path+c+'static'+filename2[1]+'.dat'

;提取结果的输出文件名 

     openw,lun,outsum,/get_lun,/append      

     printf,lun,data_temp 

free_lun,lun 

;输出提取结果

endfor

endfor 

end 

 

3 NDVI 数据提取实例 

 

现在为了提取西藏那曲的草地的物候信息，需要逐年提取那曲地区 2000 至 2018 年的 NDVI 序列。我们下载了相应的遥感影像并完成了拼接和分文件夹存储。然后采用我们编写的代码提取这个序列（图 4）。提取结果如图 5，代码提取结果具有不错的周期性，能够满足物候提取的需要。

 

                                                 代码运行界面

                                        代码提取的 NDVI 序列 

 

资料来源：《高原农业 》，作者：司国新

更多遥感和GIS知识，关注我的微信公众号：遥感加油站