纯IDL实现ENVI格式遥感影像转(Big)GeoTIFF

缘起

• 采用ENVI接口，简单的代码ENVIRaster.Export就可以实现ENVI格式转向GeoTIFF格式，如果你愿意采用这种简单的方式，那本文到这就结束了。
• 其实IDL也有现成的WRITE_TIFF函数，可以直接写入GeoTIFF文件，但是对于大影像而言，你是无法直接构建矩阵来储存全部数据的：例如，一景GF6_PMS的单精度浮点型影像，构建整个影像的矩阵则会出现下面的内存问题：IDL> arr = FLTARR(4, 48696, 46164) % Unable to allocate memory: to make array. Cannot allocate memory % Execution halted at: $MAIN$。这个时候你就会想WRITE_TIFF函数有没有类似于append的参数，可以更新GeoTIFF，但是WRITE_TIFF的append关键字只是用于扩张image，根本不能实现逐行或者逐通道写入GeoTIFF。你开始头大，可能回头想想还是觉得ENVI接口比较香。但是！正如巨馍蘸酱所言：

我们遇到什么困难，也不要怕，微笑着面对它！消除恐惧的最好办法就是面对恐惧！坚持，才是胜利。加油！奥利给！

灵感

在IDL根目录的/lib/obsolete/下有个tiff_write函数，很显然已经被弃置了，但是其核心还是可取的，代码里面的大致思路就是从TIFF6.0的组织结构入手，采用标准I/O将各Tag和数据写入文件。个人猜想其被弃置的主要原因有：只能写入单色和三色，不支持浮点，无法写入地理标签，不能突破4G的常规TIFF储存限制。所以我花了大概一周的时间来读懂代码+学TIFF底层结构+优化到支持多通道/浮点/地理标签/BigTIFF+再优化。关于(Geo)(Big)TIFF文件结构的参考学习资料：

TIFF Tag Reference
GeoTIFF标签说明
TIFF6官方说明文档(pdf)
BigTIFF说明

不足

这个程序只是用来优化白鸽的，所以存在以下缺陷，但是如果读者自己想要做一个更加完整的，其实也不是很难。

• 只做了INT/UINT/FLOAT三种数据类型
• 一些地理标签被固定
• 无法转换BIL格式的ENVI影像

使用

输入的文件名可以是fn或fn.dat，对应的头文件可以是fn.dat.hdr或fn.hdr或hdr为大写，输出默认为同目录输出，格式为fn_TIFF.tiff，或者自定义

IDL> .compile enhancedffconvert
IDL> enhancedffConvert, inputFilename
IDL> ;or
IDL> enhancedffConvert, inputFilename, outputFilename

时间测试

精力有限，只测试了一景常规GeoTIFF和一景BigGeoTIFF

储存大小	ENVIRaster.Export方法耗时(s)	本方法耗时(s)
459M	6.64 / 5.82 / 5.80	2.84 / 2.70 / 2.75
6G	324.09 / 320.28 / 385.56	280.76 / 282.21 / 268.94

源代码

pro enhancedffConvert, i_fn, o_fn
  compile_opt idl2

  if N_PARAMS() eq 1 then o_fn = i_fn + '_TIFF.tiff'

  ;look-up-table of TIFF entry name and its code
  lut = HASH('NewSubfileType', 254, $
    'ImageWidth', 256, $
    'ImageLength', 257, $
    'BitsPerSample', 258, $
    'Compression', 259, $
    'PhotometricInterpretation', 262, $
    'StripOffsets', 273, $
    'Orientation', 274, $
    'SamplesPerPixel', 277, $
    'RowsPerStrip', 278, $
    'StripByteCounts', 279, $
    'XResolution', 282, $
    'YResolution', 283, $
    'PlanarConfiguration', 284, $
    'Software', 305, $
    'DateTime', 306, $
    'ExtraSamples', 338, $
    'SampleFormat', 339, $
    'ModelPixelScale', 33550, $
    'ModelTiepoint', 33922, $
    'GeoKeyDirectory', 34735)

  readHeader, i_fn, $
    ns = ns, nl = nl, dt = dt, $
    nb = n