喜大普奔！FLAASH提供调用接口了！

请允许我先激动10分钟。。。

再平复6分钟。。。

tick tick tick…

 

今天春分。

• 春分是一年四季中阴阳平衡、昼夜均等、寒温各半的时期，所谓“春分者，阴阳相半也，故昼夜均而寒暑平”。所以我国自古就有“春分风不小，要防痛深扰”的说法，大伙别急着减衣服哦。

 

《七绝·苏醒》
　　春分雨脚落声微，柳岸斜风带客归。
　　时令北方偏向晚，可知早有绿腰肥。

 

进入主题 ^_^

 

ENVI的大气校正扩展模块提供了FLAASH和QUAC两个工具。

• 官方只提供了QUAC的调用接口，即ENVIQUACTask和ENVIQUACRaster，可以方便的进行批处理。
• FLAASH的调用接口，可以说大家已经望眼欲穿、望穿秋水、望穿秋裤。。。

 

虽然目前ENVI帮助中还没有公开FLAASH的调用接口，但是在我们真诚的打（bi）动（po）下，终于让老美提供了开发接口。

 

2017年7月24日 13:44:31更新

在ENVI54、IDL86版本调用此接口，需要在PRO代码前边加入下面几行代码才能正常调用。

 COMMON envi_acc_event_common2, wpres_toggle, $

    streams_toggle, stem_text, wa1, nb, gain, msi_arr

msi_arr = envi_get_msi_array()

 

重要！！！

如果自己在调用时报错信息类似如下，需要在执行FLAASH之前，在ENVI中关闭输入的辐射亮度数据。

 

本文将提供两个示例代码，分别处理多光谱（Landsat 8）和高光谱（AVIRIS）。

• 代码中使用的是相对路径，可直接运行。
• 输入数据在data文件夹中，输出结果放在output文件夹中。
• 虽然参数很多，但是和FLAASH面板中的参数是一一对应的。

 

源码和测试数据下载地址：http://pan.baidu.com/s/1kUETjX9

 

此次公开接口目测应该所有ENVI版本都支持！

不过示例代码中使用了新版接口，替换成对应的ENVI Classic接口即可在旧版中使用！

 

把多光谱处理的代码贴上来，敬请欣赏。

大量中文注释，让您的调用畅通无阻。

 

 

PRO test_flaash_object_class_multi

 COMPILE_OPT idl2

 e=envi()

 

 ;获取源码所在文件

  Path = FILE_DIRNAME(ROUTINE_FILEPATH())+PATH_SEP()

 

 ;输入多光谱辐射亮度数据

  radiance_file = FILEPATH('Landsat_8_OLI_Rad.dat', $

   root_dir=Path, subdirectory=['data'])

 

 ;大气校正输出结果文件路径

  reflect_file = FILEPATH('Landsat_8_OLI_Rad_flaash.dat', $

   root_dir=Path, subdirectory=['output'])

 

 ;打开栅格数据，获取元数据信息

  raster = e.OpenRaster(radiance_file)

 

 ;数据信息

  nspatial = raster.nColumns  ;列数

  nlines = raster.nRows       ;行数

  fid = ENVIRasterToFID(raster)

 ENVI_FILE_QUERY, fid, data_type=data_type ;数据类型

 

 ;获取输入文件后缀 （.dat）

  exten = stregex(radiance_file, '\..+$', /extract)

 ;对应 Rootname for FLAASH files 参数，设置为 Landsat_8_OLI_Rad_

  user_stem_name = FILE_BASENAME(radiance_file, exten)+'_'

 

 ;对应 Output Directory for FLAASH Files 参数

 ;使用输出路径作为临时路径（不建议使用系统临时路径）

 modtran_directory = Path+'output'

 

 ;获取光谱响应函数路径

 filter_func_filename = FILEPATH('landsat8_oli.sli', $

   root_dir=e.root_dir, subdirectory=['resource','filterfuncs'])

 

 ;获取时间信息

 IF OBJ_VALID(raster.Time) THEN BEGIN

   ;如果元数据中有时间信息，则自动获取

   tmpTimes = STRSPLIT(raster.Time.ACQUISITION, '-T:Z', /extract)

   year = FIX(tmpTimes[0])

   month = FIX(tmpTimes[1])

   day = FIX(tmpTimes[2])

   gmt = DOUBLE(tmpTimes[3]) + $

     DOUBLE(tmpTimes[4])/60D + DOUBLE(tmpTimes[5])/60D^2

 ENDIF ELSE BEGIN

   ;如果元数据中没有，则手动设置

   year = 2013

   month = 10

   day = 3

   gmt = 2.923418

 ENDELSE

 

 ;坐标信息

  ref = raster.SPATIALREF

 IF ref NE !NULL THEN BEGIN

   ;如果有坐标系，则自动获取经纬度、分辨率

   pixel_size = (ref.pixel_size)[0]

   ref.ConvertFileToMap, nspatial/2, nlines/2, MapX, MapY

   ref.ConvertMapToLonLat, MapX, MapY, longitude, latitude

 ENDIF ELSE BEGIN

   ;如果没有坐标系，则手动设置

   pixel_size = 30.0

   longitude = 117.08846

   latitude = 40.506906

 ENDELSE

 

 ;获取波长信息

  metadata = Raster.Metadata

 wavelength_units = metadata['WAVELENGTH UNITS']

  lambda = metadata['WAVELENGTH']

 ;fwhm如果没有，可设置值全部为-1，

 ;例如4个波段的多光谱数据，设置为[-1.0, -1.0, -1.0, -1.0]

 IF metadata.HasTag('FWHM') THEN $

   fwhm = metadata['FWHM'] $

 ELSE fwhm = DBLARR(raster.nBands)-1.0

 ;缩放系数，如果定标时设置了FLAASH Setting，则设置value=1.0即可。

  input_scale = MAKE_ARRAY(raster.nbands, value=1.0, /double)

 

 ;初始化FLAASH对象

 ;可选关键字如下：

 ; rad_remove FLAASH执行完毕后，自动关闭输入文件

 ; anc_remove FLAASH执行完毕后，自动关闭生成的辅助数据

 ; anc_delete FLAASH执行完毕后，自动关闭并删除辅助数据

  flaash_obj = obj_new('flaash_batch', /anc_delete)

 

 ;设置大量的输入参数

 flaash_obj->SetProperty, $

   hyper = 0, $ ;设置为1，表示高光谱；设置为0，表示多光谱

   ;

   ; FLAASH工程参数----

   radiance_file = radiance_file, $

   reflect_file = reflect_file, $

   filter_func_filename = filter_func_filename, $

   filter_func_file_index = 0, $

   water_band_choice = 1.13, $

   red_channel = 4, $  ;0表示undefined，LC8红波段为第4波段

   green_channel = 3, $ ;0表示undefined，LC8绿波段为第3波段

   blue_channel = 2, $ ;0表示undefined，LC8蓝波段为第2波段

 

   ;水汽反演，没有所需波段，所以设置为0，表示undefined

   ;分别对应Multispectral Setting中Water Retrieval选项卡中的两个参数

   water_retrieval = 0, $ ;Water Retrieval参数。0表示No，1表示Yes

   water_abs_channel = 0, $

   water_ref_channel = 0, $

 

   ;气溶胶反演----

   ;对应Multispectral Setting中Kaufman-Tanre Aerosol Retrieval选项卡中的参数

   kt_upper_channel = 6, $ ;设置短波红外2（SWIR 2）

   kt_lower_channel = 4, $ ;设置红波段（Red）

   kt_cutoff = 0.08, $ ;Maximum Upper Channel Reflectance

   kt_ratio = 0.500, $ ;Reflectance Ratio

   cirrus_channel = 0, $ ;0表示undefined

 

   ;前边已经定义

   user_stem_name = user_stem_name, $

   modtran_directory = modtran_directory, $

   ;

   ; MODTRAN参数---

   visvalue = 40.0000, $ ;能见度，默认40km

 

   ;为了进行水汽反演，需要如下3个波段范围中的一个：

   ; 1050-1210nm, 770-870nm, 870-1020nm

   ; 而且要求此范围的波段光谱分辨率最低为15nm

   f_resolution = 15.0000, $

 

   ;时间信息----

   day = day, $

   month = month, $

   year = year, $

   gmt = gmt, $

   latitude = latitude, $

   longitude = longitude, $

   sensor_altitude = 705.0000, $ ;传感器高度

   ground_elevation = 0.043, $ ;平均海拔，单位km

 

   ;分别对应 Advanced Setting 中的同名参数，默认即可

   view_zenith_angle = 180.0000, $

   view_azimuth = 0.0000, $

 

   ;大气模型：0-SAW；1-MLW；2-U.S. Standard；3-SAS；4-MLS；5-T

   atmosphere_model = 4, $

   ;气溶胶模型：0-No Aerosol；1-Rural；2-Maritime；3-Urban；4-Tropospheric

   aerosol_model = 3, $

 

   ;如下几个参数对应 Advanced Setting同名参数，默认即可。

   multiscatter_model = 0, $

   disort_streams = 8, $

   co2mix = 390.0000, $

   water_column_multiplier = 1.0000, $

   ;

   ;图像参数----

   nspatial = nspatial, $

   nlines = nlines, $

   data_type = data_type, $

   margin1 = 0, $

   margin2 = 0, $

   nskip = 0, $

   pixel_size = pixel_size, $

   sensor_name = 'Landsat-8 OLI', $

 

   ;分析参数----

   ;对应Advanced Setting中的 Aerosol Scale Height

   aerosol_scaleht = 2.0000, $

   ;对应Advanced Setting中的 Use Adjacency Correction

   ;中高分辨率设置为1，低分辨率（如Modis）设置为0

   use_adjacency = 1, $

 

   ;输出缩放系数，输出结果放大了10000倍，变为UINT数据类型。

   ;对应Advanced Setting中的Output Reflectance Scale Factor

   output_scale = 10000.0000, $ ;输出结果缩放系数

 

   ;对应 Width (number of bands) 参数，多光谱设置0即可。

   polishing_res = 0, $

 

   ;对应 Aerosol Retrieval 参数。

   ; 0 表示 None；1 表示 2-Band (K-T)；2 表示 2-Band Over Water

   aerosol_retrieval = 1, $

 

   ;对应FLAASH面板中的 Wavelength Recalibration，多光谱一般为0

   calc_wl_correction = 0, $

   reuse_modtran_calcs = 0, $

   use_square_slit_function = 0, $

   convolution_method = 'fft', $

 

   ;对应Advanced Setting中的 Use Tiled Processing

   ;1-Yes；0-No

   use_tiling = 0, $

   tile_size = 1024.0000, $

 

   ; Spectral Parameters

   wavelength_units = wavelength_units, $

   lambda = lambda, $

   fwhm = fwhm, $

   input_scale = input_scale

 

 ;重要！！！！！！重要！！！！！！重要！！！！！

 ;执行FLAASH之前，必须在ENVI中把输入文件关闭

 Raster.Close

 

 ;开始执行FLAASH

 flaash_obj->processImage

 

 ;获取输入输出文件的FID

 flaash_obj->getResults, rad_fid=rad_fid, reflect_fid=reflect_fid

END