缘起
- 采用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文件结构的参考学习资料：
不足
这个程序只是用来优化白鸽的，所以存在以下缺陷，但是如果读者自己想要做一个更加完整的，其实也不是很难。
- 只做了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 = nb, inter = inter, $
mapinfo = mapinfo
mapinfo = STRTOK(mapinfo, ',', /EXTRACT)
OPENR, ilun, i_fn, /GET_LUN
OPENW, lun, o_fn, /GET_LUN
;if file is bigger than 4G, then write BigTIFF
stat = FSTAT(ilun)
szinGb = (stat.SIZE) / (1024L^3)
if szinGb ge 4 then begin;-----begin of BigTIFF
common big, ifdB, _entryCode_
_entryCode_ = lut
;Image File Header at the head of the file
header = BYTE([73, 73, 43, 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])
WRITEU, lun, header
;initialize Image File Directory(IFD)
ifdB = {NumberOfEntry: BYTARR(8)}
ifdAddEntry_Big, 'NewSubfileType', UINT(0)
ifdAddEntry_Big, 'ImageWidth', ULONG64(ns)
ifdAddEntry_Big, 'ImageLength', ULONG64(nl)
;Bits Per Sample
POINT_LUN, -lun, pos
ifdAddEntry_Big, 'BitsPerSample', ULONG64(pos), tiff_dt = 3, n = nb
WRITEU, lun, REPLICATE(dt eq 4 or dt eq 5 ? UINT(32) : UINT(16), nb)
ifdAddEntry_Big, 'Compression', UINT(1)
ifdAddEntry_Big, 'PhotometricInterpretation', UINT(2)
;write data
case inter of
'bip': begin
POINT_LUN, -lun, pos
ifdAddEntry_Big, 'StripOffsets', ULONG64(pos)
arr = BYTARR(ULONG64(4E9)) & i = 0
while not EOF(ilun) do begin
READU, ilun, arr
CATCH, err
if err ne 0 then arr = BYTARR(stat.SIZE - ULONG64(4E9) * i)
WRITEU, lun, arr
i++
endwhile
arr = !NULL
end
'bsq': begin
POINT_LUN, -lun, pos
ifdAddEntry_Big, 'StripOffsets', ULONG64(pos), tiff_dt = 16, n = nb
WRITEU, lun, MAKE_ARRAY(nb, /UL64)
POINT_LUN, -lun, addr
POINT_LUN, lun, ULONG64(pos)
WRITEU, lun, MAKE_ARRAY(nb, start = ULONG64(addr), increment = $
ULONG64(ns) * nl * (dt eq 4 ? 4 : 2), /UL64, /INDEX)
arr = BYTARR(ULONG64(4E9)) & i = 0
while not EOF(ilun) do begin
READU, ilun, arr
CATCH, err
if err ne 0 then arr = BYTARR(stat.SIZE - ULONG64(4E9) * i)
WRITEU, lun, arr
i++
endwhile
arr = !NULL
end
else: RETURN
endcase
FREE_LUN, ilun
ifdAddEntry_Big, 'Orientation', UINT(1)
ifdAddEntry_Big, 'SamplesPerPixel', ULONG64(nb)
ifdAddEntry_Big, 'RowsPerStrip', ULONG64(nl)
;Strip Byte Counts
if inter eq 'bip' then $
ifdAddEntry_Big, 'StripByteCounts', $
ULONG64(nb) * ns * nl * (dt eq 4 or dt eq 5 ? 4 : 2) $
else begin
POINT_LUN, -lun, pos
ifdAddEntry_Big, 'StripByteCounts', ULONG64(pos), tiff_dt = 16, n = nb
WRITEU, lun, REPLICATE(ULONG64(ns) * nl * (dt eq 4 or dt eq 5 ? 4 : 2), nb)
endelse
ifdAddEntry_Big, 'XResolution', UINT(1)
ifdAddEntry_Big, 'YResolution', UINT(1)
ifdAddEntry_Big, 'PlanarConfiguration', inter eq 'bip' ? UINT(1) : UINT(2)
;Software
POINT_LUN, -lun, pos
str = 'Snowy Dove v1.3-alpha, Build Date 2020/04/17'
ifdAddEntry_Big, 'Software', ULONG64(pos), tiff_dt = 2, n = STRLEN(str) + 1
WRITEU, lun, [BYTE(str), 0B]
;Date & Time in YYYY:MM:DD HH:MM:SS format
POINT_LUN, -lun, pos
asciiTime = STRTOK(SYSTIME(0), /EXTRACT)
str = STRJOIN([asciiTime[4], $
STRING(1 + WHERE(STRUPCASE(asciiTime[1]) eq $
['JAN','FEB','MAR','APR', 'MAY', 'JUN', 'JUL', 'AUG','SEP','OCT','NOV','DEC']), $
format = '(I02)'), STRING(asciiTime[2], format = '(I02)')], $
':') + ' ' + asciiTime[3]
ifdAddEntry_Big, 'DateTime', ULONG64(pos), tiff_dt = 2, n = STRLEN(str) + 1
WRITEU, lun, [BYTE(str), 0B]
;Extra Samples
if nb gt 3 and nb le 5 then $
ifdAddEntry_Big, 'ExtraSamples', REPLICATE(UINT(0), nb-3) $
else begin
POINT_LUN, -lun, pos
ifdAddEntry_Big, 'ExtraSamples', ULONG64(pos), tiff_dt = 3, n = nb-3
WRITEU, lun, REPLICATE(UINT(0), nb-3)
endelse
;Sample Format
if inter eq 'bip' then $
ifdAddEntry_Big, 'SampleFormat', $
(HASH(2, UINT(2), 4, UINT(3), 12, UINT(1)))[dt] $
else begin
POINT_LUN, -lun, pos
ifdAddEntry_Big, 'SampleFormat', ULONG64(pos), tiff_dt = 3, n = nb
WRITEU, lun, REPLICATE((HASH(2, UINT(2), 4, UINT(3), 12, UINT(1)))[dt], nb)
endelse
;GeoTag: Model Pixel Scale
POINT_LUN, -lun, pos
ifdAddEntry_Big, 'ModelPixelScale', ULONG64(pos), tiff_dt = 12, n = 3
WRITEU, lun, [DOUBLE(mapinfo[5]), DOUBLE(mapinfo[6]), 0]
;GeoTag: Model Tiepoint
POINT_LUN, -lun, pos
ifdAddEntry_Big, 'ModelTiepoint', ULONG64(pos), tiff_dt = 12, n = 6
WRITEU, lun, [0, 0, 0, DOUBLE(mapinfo[3]), DOUBLE(mapinfo[4]), 0]
;GeoTag: GeoKey Directory
POINT_LUN, -lun, pos
ifdAddEntry_Big, 'GeoKeyDirectory', ULONG64(pos), tiff_dt = 3, n = 24
WRITEU, lun, UINT([1, 1, 2, 6, $
1024, 0, 1, 1, $
1025, 0, 1, 1, $
2054, 0, 1, 9102, $
3072, 0, 1, UINT((STRLOWCASE(STRTRIM(mapinfo[8], 2)) eq 'north' ? $
'326' : '327') + STRTRIM(mapinfo[7], 2)), $
3076, 0, 1, 9001])
;rewrite number of entries in IFD
ifdB.NumberOfEntry = BYTE(ULONG64(N_TAGS(ifdB) - 1), 0, 8)
;and add last bytes
ifdB = CREATE_STRUCT(ifdB, 'EndOfEntry', BYTARR(8))
;write IFD and header
POINT_LUN, -lun, pos
WRITEU, lun, ifdB
POINT_LUN, lun, 0
header[8] = BYTE(ULONG64(pos), 0, 8)
WRITEU, lun, header
;and then finish
FREE_LUN,lun
endif else begin;-----begin of normal TIFF
common normal, ifd, _entryCode
_entryCode = lut
header = BYTE([73, 73, 42, 0, 0, 0, 0, 0])
WRITEU, lun, header
ifd = {NumberOfEntry: BYTARR(2)}
ifdAddEntry_Normal, 'NewSubfileType', UINT(0)
ifdAddEntry_Normal, 'ImageWidth', ULONG(ns)
ifdAddEntry_Normal, 'ImageLength', ULONG(nl)
POINT_LUN, -lun, pos
ifdAddEntry_Normal, 'BitsPerSample', ULONG(pos), tiff_dt = 3, n = nb
WRITEU, lun, REPLICATE(dt eq 4 or dt eq 5 ? UINT(32) : UINT(16), nb)
ifdAddEntry_Normal, 'Compression', UINT(1)
ifdAddEntry_Normal, 'PhotometricInterpretation', UINT(2)
case inter of
'bip': begin
POINT_LUN, -lun, pos
ifdAddEntry_Normal, 'StripOffsets', ULONG(pos)
arr = MAKE_ARRAY(stat.SIZE, /BYTE)
READU, ilun, arr
WRITEU, lun, arr
arr = !NULL
end
'bsq': begin
POINT_LUN, -lun, pos
ifdAddEntry_Normal, 'StripOffsets', ULONG(pos), tiff_dt = 4, n = nb
WRITEU, lun, MAKE_ARRAY(nb, /ULONG)
POINT_LUN, -lun, addr
POINT_LUN, lun, ULONG(pos)
WRITEU, lun, MAKE_ARRAY(nb, start = ULONG(addr), increment = $
ULONG(ns) * nl * (dt eq 4 ? 4 : 2), /ULONG, /INDEX)
arr = MAKE_ARRAY(stat.SIZE, /BYTE)
READU, ilun, arr
WRITEU, lun, arr
arr = !NULL
end
else: RETURN
endcase
FREE_LUN, ilun
ifdAddEntry_Normal, 'Orientation', UINT(1)
ifdAddEntry_Normal, 'SamplesPerPixel', ULONG(nb)
ifdAddEntry_Normal, 'RowsPerStrip', ULONG(nl)
if inter eq 'bip' then $
ifdAddEntry_Normal, 'StripByteCounts', $
ULONG(nb) * ns * nl * (dt eq 4 or dt eq 5 ? 4 : 2) $
else begin
POINT_LUN, -lun, pos
ifdAddEntry_Normal, 'StripByteCounts', ULONG(pos), tiff_dt = 4, n = nb
WRITEU, lun, REPLICATE(ULONG(ns) * nl * (dt eq 4 ? 4 : 2), nb)
endelse
ifdAddEntry_Normal, 'XResolution', UINT(1)
ifdAddEntry_Normal, 'YResolution', UINT(1)
ifdAddEntry_Normal, 'PlanarConfiguration', inter eq 'bip' ? UINT(1) : UINT(2)
POINT_LUN, -lun, pos
str = 'Snowy Dove v1.3-alpha, Build Date 2020/04/17'
ifdAddEntry_Normal, 'Software', ULONG(pos), tiff_dt = 2, n = STRLEN(str) + 1
WRITEU, lun, [BYTE(str), 0B]
POINT_LUN, -lun, pos
asciiTime = STRTOK(SYSTIME(0), /EXTRACT)
str = STRJOIN([asciiTime[4], $
STRING(1 + WHERE(STRUPCASE(asciiTime[1]) eq $
['JAN','FEB','MAR','APR', 'MAY', 'JUN', 'JUL', 'AUG','SEP','OCT','NOV','DEC']), $
format = '(I02)'), STRING(asciiTime[2], format = '(I02)')], $
':') + ' ' + asciiTime[3]
ifdAddEntry_Normal, 'DateTime', ULONG(pos), tiff_dt = 2, n = STRLEN(str) + 1
WRITEU, lun, [BYTE(str), 0B]
if nb gt 3 and nb le 5 then $
ifdAddEntry_Normal, 'ExtraSamples', REPLICATE(UINT(0), nb-3) $
else begin
POINT_LUN, -lun, pos
ifdAddEntry_Normal, 'ExtraSamples', ULONG(pos), tiff_dt = 3, n = nb-3
WRITEU, lun, REPLICATE(UINT(0), nb-3)
endelse
if inter eq 'bip' then $
ifdAddEntry_Normal, 'SampleFormat', $
(HASH(2, UINT(2), 4, UINT(3), 12, UINT(1)))[dt] $
else begin
POINT_LUN, -lun, pos
ifdAddEntry_Normal, 'SampleFormat', ULONG(pos), tiff_dt = 3, n = nb
WRITEU, lun, REPLICATE((HASH(2, UINT(2), 4, UINT(3), 12, UINT(1)))[dt], nb)
endelse
POINT_LUN, -lun, pos
ifdAddEntry_Normal, 'ModelPixelScale', ULONG(pos), tiff_dt = 12, n = 3
WRITEU, lun, [DOUBLE(mapinfo[5]), DOUBLE(mapinfo[6]), 0]
POINT_LUN, -lun, pos
ifdAddEntry_Normal, 'ModelTiepoint', ULONG(pos), tiff_dt = 12, n = 6
WRITEU, lun, [0, 0, 0, DOUBLE(mapinfo[3]), DOUBLE(mapinfo[4]), 0]
POINT_LUN, -lun, pos
ifdAddEntry_Normal, 'GeoKeyDirectory', ULONG(pos), tiff_dt = 3, n = 24
WRITEU, lun, UINT([1, 1, 2, 6, $
1024, 0, 1, 1, $
1025, 0, 1, 1, $
2054, 0, 1, 9102, $
3072, 0, 1, UINT((STRLOWCASE(STRTRIM(mapinfo[8], 2)) eq 'north' ? $
'326' : '327') + STRTRIM(mapinfo[7], 2)), $
3076, 0, 1, 9001])
ifd.NumberOfEntry = BYTE(N_TAGS(ifd) - 1, 0, 2)
ifd = CREATE_STRUCT(ifd, 'EndOfEntry', BYTARR(4))
POINT_LUN, -lun, pos
WRITEU, lun, ifd
POINT_LUN, lun, 0
header[4] = BYTE(ULONG(pos), 0, 4)
WRITEU, lun, header
FREE_LUN,lun
endelse
end
;add certain entry to IFD
pro ifdAddEntry_Normal, name, data, tiff_dt = tiff_dt, n = n
compile_opt idl2
common normal, ifd, lut
newEntry = BYTARR(12)
;Tag
newEntry[0] = BYTE(lut[name], 0, 2)
;Type
if not ISA(tiff_dt) then begin
idl_dt = SIZE(data, /TYPE)
tiff_dt = (HASH($
1, 1, $; idl_byte -> TIFF_byte
4, 11, $; idl_float -> TIFF_float
5, 12, $; idl_double -> TIFF_double
12, 3, $; idl_uint -> TIFF_short
13, 4)) $; idl_ulong -> TIFF_long
[idl_dt]
end
newEntry[2] = BYTE(tiff_dt, 0, 2)
;Count
if not ISA(n) then n = N_ELEMENTS(data)
newEntry[4] = BYTE(ULONG(n), 0, 4)
;Value or Offset
newEntry[8] = BYTE(ULONG(data), 0, 4)
ifd = CREATE_STRUCT(ifd, name, newEntry)
end
;add certain entry to IFD of BigTIFF
pro ifdAddEntry_Big, name, data, tiff_dt = tiff_dt, n = n
compile_opt idl2
common big, ifdB, lut
newEntry = BYTARR(20)
;Tag
newEntry[0] = BYTE(lut[name], 0, 2)
;Type
if not ISA(tiff_dt) then begin
idl_dt = SIZE(data, /TYPE)
tiff_dt = (HASH($
1, 1, $; idl_byte -> TIFF_byte
4, 11, $; idl_float -> TIFF_float
5, 12, $; idl_double -> TIFF_double
12, 3, $; idl_uint -> TIFF_short
13, 4, $; idl_ulong -> TIFF_long
14, 17, $; idl_long64 -> TIFF_slong8(BigTIFF)
15, 16)) $; idl_ulong64 -> TIFF_long8(BigTIFF)
[idl_dt]
end
newEntry[2] = BYTE(tiff_dt, 0, 2)
;Count
if not ISA(n) then n = N_ELEMENTS(data)
newEntry[4] = BYTE(ULONG64(n), 0, 8)
;Value or Offset
newEntry[12] = BYTE(ULONG64(data), 0, 8)
ifdB = CREATE_STRUCT(ifdB, name, newEntry)
end
;read header file of ENVI format image
;and return useful information
pro readHeader, fn_bin, $
ns = ns, nl = nl, dt = dt, $
nb = nb, inter = inter, $
mapinfo = mapinfo, cs = cs, $
header = header, fn_hdr = fn_hdr
compile_opt idl2, hidden
appendHeader = ARG_PRESENT(header)
basename = FILE_BASENAME(fn_bin, '.dat')
basenameWithSuffix = FILE_BASENAME(fn_bin)
_fn_hdr = FILEPATH([basename + '.hdr', $
basenameWithSuffix + '.hdr', $
basename + '.HDR', $
basenameWithSuffix + '.HDR'], $
root_dir = FILE_DIRNAME(fn_bin))
fn_hdr = _fn_hdr[(WHERE(FILE_TEST(_fn_hdr) eq 1))[0]]
OPENR, lun, fn_hdr, /GET_LUN
tmp = ''
if appendHeader then header = []
while not EOF(lun) do begin
READF, lun, tmp
if appendHeader then header = [header, tmp]
case 1 of
STRPOS(tmp, 'samples') ne -1: $
ns = LONG(STREGEX(tmp, '[0-9]+', /EXTRACT))
STRPOS(tmp, 'lines') ne -1: $
nl = LONG(STREGEX(tmp, '[0-9]+', /EXTRACT))
STRPOS(tmp, 'data type') ne -1: $
dt = LONG(STREGEX(tmp, '[0-9]+', /EXTRACT))
STRPOS(tmp, 'bands') ne -1: $
nb = LONG(STREGEX(tmp, '[0-9]+', /EXTRACT))
STRPOS(tmp, 'interleave') ne -1: $
inter = STRMID(tmp, STRLEN(tmp) - 3, 3)
STRPOS(tmp, 'map info') ne -1: $
mapinfo = STRMID(tmp, STREGEX(tmp, '\{'))
STRPOS(tmp, 'coordinate system string') ne -1: $
cs = STRMID(tmp, STREGEX(tmp, '\{'))
else: ;pass
endcase
endwhile
FREE_LUN, lun
end
后记
- 可能用C语言写这个功能速度会更快一点
- 头秃。
2755
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
