基于烟羽掩膜的森林火点MODIS自动探测软件

() 使用说明 编写人员朱源 编写时间 2010 年 3 月 Email:

目 录 1. 目的及意义...1 1.1 编写目的...1 1.2 作用及意义...1 2 运行环境...2 2.1 设备...2 2.2 支持系统...2 2.3 数据要求...2 3. 系统使用...3 4. 测试数据...6 5.IDL 代码...10

基于 IDL 的多图层栅格数据分区统计模块 使用说明 1. 目的及意义 1.1 编写目的 本手册提供了 基于 IDL 的多图层栅格数据分区统计模块 的基本流程 软硬件环境 安装和使用说明, 以供使用人员参考 1.2 作用及意义 图 1 ArcGIS 中的 模块分区统计在地学研究中经常用到,ArcGIS 软件提供的 软件模块 ( 如图 1 所示 ) 可以实现该功能 但是该模块仅能实现单图层数据的操作, 难以处理多图层数据 而实际应用中, 往往需要大批量地对多个图层的数据同时进行分区统计, 如统计某个区域内某种地物类型中多个年份 12 个月的 NDVI 情况 这时如果用 ArcGIS 来处理, 工作量十分巨大, 而且效率低下 基于这种情况结合我们自己的实际需要, 我们编写了这个基于 IDL 的多图层栅格数据分区统计 1

模块, 可供需要的人员使用和参考 2 运行环境 如下 应用软件系统的运行需要一定的软硬件环境支持, 本系统要求的硬件环境 2.1 设备 普通微机, 基本配置是 : CPU: 主频 1.6G 以上 显示卡 : 中等性能显卡, 显存在 32M 以上 ; 显示器 :17 纯屏显示器, 显示分辨率 1024*768 以上 ; 硬盘 : 至少 30G 以上 内存 :512M 以上, 建议使用 1G 内存 本中文环境下运行, 操作系统使用中文版 Windows XP 2.2 支持系统 操作系统环境 : Microsoft Windows XP 开发工具 : IDL6.3; 运行环境 : ENVI4.3+ IDL6.3 2.3 数据要求 输入的数据要求为 1-12 个图层待统计的数据层 ( 如 NDVI) 和 1 个类型掩膜层 ( 如地表类型或行政区 ) 叠合而成的具有 13 个波段的图像文件 ( 如图 2 所示 ),, 该步骤可以 ENVI 中通过 Layer Stacking 来实现 具体使用时可根据实际情况修改 IDL 代码来实现特定需求 2

图 2 统计模块输入数据格式示例 3. 系统使用 步骤一 : 双击启动 Zonal_Statistics IDL 程序 ; 图 3 启动 IDL 程序 步骤二 : 在 IDL 命令行输入 ENVI, 启动 ENVI 软件 ; 步骤三 : 编译和运行程序 ; 图 4 启动 ENVI 程序 3

图 5 编译和运行 IDL 程序 步骤四 : 输入和输出 ; 4

图 6 输入和输出界面 5

4. 测试数据 图 7 测试数据内容 meantestdata 压缩包中数据文件内容说明如下 :gridtest32.tif 是要统计的栅格 6

图层数据 ( 如 NDVI 等 ), 如图 8 所示 ;gridtest8.tif 是用于作类型掩膜的图层数据 ( 如地表类型等 ), 如图 9 所示 ;gridtestdata32_13band.tif 是 12 个图层的 NDVI+1 个图层的类型掩膜共 13 个图层的叠置数据, 如图 2 所示 ;gridtest32shp.shp 是统计的栅格图层 ( 如 NDVI 等 ) 转换为 POLYGON 后的矢量文件, 如图 8 所示 ; gridtest8shp.shp 是用于作类型掩膜的图层数据 ( 如地表类型等 ) 转换为 POLYGON 后的矢量文件, 如图 9 所示 ;ZonalSt_gridtes1.dbf 是 ArcGIS 所作的 结果, 如图 11 所示, 用于 IDL 计算结果的比较 ;ndvi stastic result.txt 是 IDL 的分区统计结果, 如图 12 所示 ;Zonal_Statistics 是 IDL 程序源代码文件 图 8 gridtest32.tif 是要统计的栅格图层数据 ( 如 NDVI 等 ) 7

图 9 gridtest8.tif 是用于作类型掩膜的图层数据 ( 如地表类型等 ) 图 10 用于分区统计的两个数据图层叠置显示 8

图 11 ZonalSt_gridtes1.dbf 是 ArcGIS 所作的 结果 图 12 ndvi stastic result.txt 是 IDL 的分区统计结果 9

5.IDL 代码 Zonal_Statistics 是 IDL 程序源代码文件 : PRO Zonal_Statistics;Zonal_Statistics ; input 12+1 bands image of ndvi ENVI_SELECT, fid=ref_fid, dims=ref_dims, title='select ndvi 12 months+classtype Image' if ref_fid eq -1 then return ;defien the output text file base1 = widget_auto_base(title='enter the name of output text file') ;wo = widget_outf(base, uvalue='outf', /auto) wo1 = widget_outf(base1,auto_manage=1,default='ndvi stastic', uvalue='outf') result1 = auto_wid_mng(base1) if (result1.accept eq 0) then return print, 'Selected File path and name:', result1.outf outfile=result1.outf start_time = systime(1) ; Show status bar envi_report_init,'in progress,please wait!', title='zonal_statistics Calculation',base=base ; Update status bar envi_report_stat,base,1,10 ; get the reflectance of band 1-12 NDVI1=envi_get_data(fid=ref_fid, pos=0, dims=ref_dims); NDVI2=envi_get_data(fid=ref_fid, pos=1, dims=ref_dims); NDVI3=envi_get_data(fid=ref_fid, pos=2, dims=ref_dims); NDVI4=envi_get_data(fid=ref_fid, pos=3, dims=ref_dims); NDVI5=envi_get_data(fid=ref_fid, pos=4, dims=ref_dims); NDVI6=envi_get_data(fid=ref_fid, pos=5, dims=ref_dims); NDVI7=envi_get_data(fid=ref_fid, pos=6, dims=ref_dims); NDVI8=envi_get_data(fid=ref_fid, pos=7, dims=ref_dims); NDVI9=envi_get_data(fid=ref_fid, pos=8, dims=ref_dims); NDVI10=envi_get_data(fid=ref_fid, pos=9, dims=ref_dims); NDVI11=envi_get_data(fid=ref_fid, pos=10, dims=ref_dims); NDVI12=envi_get_data(fid=ref_fid, pos=11, dims=ref_dims); class=envi_get_data(fid=ref_fid, pos=12, dims=ref_dims); ; Update status bar envi_report_stat,base,1,10 ; Write to output text file 10

outfile=strmid(outfile,0,22);the first character position is 0.,22 is Length outfile=string(outfile,'.txt') print,'outfile=',outfile openw,unit,outfile,/get_lun ; 定义植被种类数 classnumber=18;18 是植被种类数量 printf,unit,'ndvi stastic result' ; 定义数组 classcount=intarr(classnumber+1) classcode=intarr(classnumber+1) printf,unit,' month, classtype, pixelscount, mean, max, min' ;Begin cycle of all image per pixels for i=1,classnumber do begin ;ndvi month1 if count gt 0 then ndvimean1 = mean(ndvi1[index]) if count gt 0 then ndvimax1 =max(ndvi1[index]) if count gt 0 then ndvimin1 =min(ndvi1[index]) if count gt 0 then classcount[i] =n_elements(ndvi1[index]); 统计 class[i] 类型的像元个数 if count gt 0 then printf,unit,1,',',i,',',classcount[i],',',ndvimean1,',',ndvimax1,',',ndvimin1 ;ndvi month2 if count gt 0 then ndvimean2 = mean(ndvi2[index]) if count gt 0 then ndvimax2 =max(ndvi2[index]) if count gt 0 then ndvimin2 =min(ndvi2[index]) if count gt 0 then classcount[i] =n_elements(ndvi2[index]); 统计 class[i] 类型的像元个数 if count gt 0 then printf,unit,2,',',i,',',classcount[i],',',ndvimean2,',',ndvimax2,',',ndvimin2 ;ndvi month3 if count gt 0 then ndvimean3 = mean(ndvi3[index]) if count gt 0 then ndvimax3 =max(ndvi3[index]) if count gt 0 then ndvimin3 =min(ndvi3[index]) if count gt 0 then classcount[i] =n_elements(ndvi3[index]); 统计 class[i] 类型的像元个数 11

if count gt 0 then printf,unit,3,',',i,',',classcount[i],',',ndvimean3,',',ndvimax3,',',ndvimin3 ;ndvi month4 if count gt 0 then ndvimean4 = mean(ndvi4[index]) if count gt 0 then ndvimax4 =max(ndvi4[index]) if count gt 0 then ndvimin4 =min(ndvi4[index]) if count gt 0 then classcount[i] =n_elements(ndvi4[index]); 统计 class[i] 类型的像元个数 if count gt 0 then printf,unit,4,',',i,',',classcount[i],',',ndvimean4,',',ndvimax4,',',ndvimin4 ;ndvi month5 if count gt 0 then ndvimean5 = mean(ndvi5[index]) if count gt 0 then ndvimax5 =max(ndvi5[index]) if count gt 0 then ndvimin5 =min(ndvi5[index]) if count gt 0 then classcount[i] =n_elements(ndvi5[index]); 统计 class[i] 类型的像元个数 if count gt 0 then printf,unit,5,',',i,',',classcount[i],',',ndvimean5,',',ndvimax5,',',ndvimin5 ;ndvi month6 if count gt 0 then ndvimean6 = mean(ndvi6[index]) if count gt 0 then ndvimax6 =max(ndvi6[index]) if count gt 0 then ndvimin6 =min(ndvi6[index]) if count gt 0 then classcount[i] =n_elements(ndvi6[index]); 统计 class[i] 类型的像元个数 if count gt 0 then printf,unit,6,',',i,',',classcount[i],',',ndvimean6,',',ndvimax6,',',ndvimin6 ; Update status bar envi_report_stat,base,5,10 ;ndvi month7 if count gt 0 then ndvimean7 = mean(ndvi7[index]) if count gt 0 then ndvimax7 =max(ndvi7[index]) if count gt 0 then ndvimin7 =min(ndvi7[index]) if count gt 0 then classcount[i] =n_elements(ndvi7[index]); 统计 class[i] 类型的像元个数 if count gt 0 then printf,unit,7,',',i,',',classcount[i],',',ndvimean7,',',ndvimax7,',',ndvimin7 ;ndvi month8 12

if count gt 0 then ndvimean8 = mean(ndvi8[index]) if count gt 0 then ndvimax8 =max(ndvi8[index]) if count gt 0 then ndvimin8 =min(ndvi8[index]) if count gt 0 then classcount[i] =n_elements(ndvi8[index]); 统计 class[i] 类型的像元个数 if count gt 0 then printf,unit,8,',',i,',',classcount[i],',',ndvimean8,',',ndvimax8,',',ndvimin8 ;ndvi month9 if count gt 0 then ndvimean9 = mean(ndvi9[index]) if count gt 0 then ndvimax9 =max(ndvi9[index]) if count gt 0 then ndvimin9 =min(ndvi9[index]) if count gt 0 then classcount[i] =n_elements(ndvi9[index]); 统计 class[i] 类型的像元个数 if count gt 0 then printf,unit,9,',',i,',',classcount[i],',',ndvimean9,',',ndvimax9,',',ndvimin9 ;ndvi month10 if count gt 0 then ndvimean10 = mean(ndvi10[index]) if count gt 0 then ndvimax10 =max(ndvi10[index]) if count gt 0 then ndvimin10 =min(ndvi10[index]) if count gt 0 then classcount[i] =n_elements(ndvi10[index]); 统计 class[i] 类型的像元个数 if count gt 0 then printf,unit,10,',',i,',',classcount[i],',',ndvimean10,',',ndvimax10,',',ndvimin10 ;ndvi month11 if count gt 0 then ndvimean11 = mean(ndvi11[index]) if count gt 0 then ndvimax11 =max(ndvi11[index]) if count gt 0 then ndvimin11 =min(ndvi11[index]) if count gt 0 then classcount[i] =n_elements(ndvi11[index]); 统计 class[i] 类型的像元个数 if count gt 0 then printf,unit,11,',',i,',',classcount[i],',',ndvimean11,',',ndvimax11,',',ndvimin11 ;ndvi month12 if count gt 0 then ndvimean12 = mean(ndvi12[index]) if count gt 0 then ndvimax12 =max(ndvi12[index]) if count gt 0 then ndvimin12 =min(ndvi12[index]) if count gt 0 then classcount[i] =n_elements(ndvi12[index]); 统计 class[i] 类型的像元个数 if count gt 0 then printf,unit,12,',',i,',',classcount[i],',',ndvimean12,',',ndvimax12,',',ndvimin12 endfor;end cycle of all image per pixels 13

free_lun,unit;end of GET_LUN, ; Update status bar envi_report_stat,base,10,10 ;out put image to memory ;out_map_info = envi_get_map_info(fid=ref_fid) ;envi_enter_data,ndvi1,bnames='ndvi1', map_info=out_map_info print,'process Time :', systime(1) - start_time ; Close status bar envi_report_init, base=base, /finish END Result = DIALOG_MESSAGE('Process Done!',/INFORMATION, TITLE='Zonal_Statistics' ) 14