|
|
|
/* ---------------------------------------------------------------- */
/* Transforms a GeoTiff file as produced by ASTER GDEM, see: */
/* */
/* http://gdem.ersdac.jspacesystems.or.jp */
/* */
/* to 'hgt' format as used by SRTM, see: */
/* */
/* https://lpdaac.usgs.gov/products/measures_products_table/srtmgl3 */
/* */
/* Foa a nice explanation of the (geo)tiff format see: */
/* */
/* http://www.awaresystems.be/imaging/tiff.html */
/* */
/* ---------------------------------------------------------------- */
/* */
/* Originally by Ruurd J. Idenburg */
/* */
/* No copyright, no licence, no guarantees or warrantees, be it */
/* explicit, implicit or whatever. Usage is totally and completely */
/* at the users own risk, the author shall not be liable for any */
/* damages whatsoever, for any reason whatsoever. */
/* */
/* Please keep this comment block intact when modifying this code */
/* and add a note with date and a description. */
/* */
/* ---------------------------------------------------------------- */
/* Parameter(s): */
/* */
/* file - the (geo)tiff path and filename to be converted */
/* */
/* */
/* Result: */
/* */
/* A file on the same path and filename as the input tiff file */
/* but with file extension 'hgt' and in the SRTM format */
/* */
/* ---------------------------------------------------------------- */
/* 2014/02/01 - Initial version */
/* ---------------------------------------------------------------- */
use arg geoTiffFile
geoTiff = .stream~new(geoTiffFile)
geoTiff~open("read")
/* ---------------------------------------------------------------- */
/* The Image File Directory(IFD) is always at offset 4 in the file. */
/* ---------------------------------------------------------------- */
ifdOffset = (geoTiff~charin(5,4)~reverse~c2d)+1 -- cater for base 1
say "Found the IFD at offset" ifdOffset "('"ifdOffset~d2x(8)"'X)."
/* ---------------------------------------------------------------- */
/* Get the number of tags present in this IFD. */
/* ---------------------------------------------------------------- */
tagCount = geoTiff~charin(ifdOffset,2)~reverse~c2d
say "IFD contains" tagCount "tags."
/* ---------------------------------------------------------------- */
/* The tags follow directly after the count and each 12 bytes long. */
/* ---------------------------------------------------------------- */
tags = geoTiff~charin(ifdOffset+2,tagCount*12)
nextIfdOffset = geoTiff~charin(ifdOffset+2+tagCount*12,4)~reverse~c2d
-- say nextIfdOffset -- for now assume no following IFD's
tagDir = .directory~new
/* ---------------------------------------------------------------- */
/* Each tag consists of Key, Datatype, Valuecount and Value(s). */
/* If the values don't fit in 4 bytes, the value is an offset to */
/* where the values can be found. */
/* ---------------------------------------------------------------- */
do t=0 to tagCount-1 --cater for base 1
aTag = tags~substr(1+12*t,12)
tagKey = aTag~substr(1,2)~reverse~c2x
dataType = aTag~substr(3,2)~reverse~c2x
valueCount = aTag~substr(5,4)~reverse~c2d
valuePointer = aTag~substr(9,4)~reverse~c2d -- can be a value only
tag = .directory~new
tagDir~put(dataType valueCount valuePointer,tagKey)
end
columns = tagDir~0100~word(3)
say "The elevation information is" columns "columns wide."
rows = tagDir~0101~word(3)
say "The elevation information is" rows "rows deep."
--rowsPerStrip = tagDir~0116~word(3) -- not needed for further processing
--say rowsPerStrip
stripByteCountValues = tagDir~0117~word(2)
say "There are" stripByteCountValues "byte counts for the row strips to read."
stripByteCountsOffset = tagDir~0117~word(3)+1 -- cater for base 1
--say stripByteCountsOffset
stripByteCounts = geoTiff~charin(stripByteCountsOffset,stripByteCountValues*4)
--say stripByteCounts~c2x
stripCount = tagDir~0111~word(2)
say "The strip count" stripCount "should be equal to previous value."
stripOffsetsPointer = tagDir~0111~word(3)+1 -- cater for base 1
--say stripOffsetsPointer
stripPointers = geoTiff~charin(stripOffsetsPointer,stripCount*4)
heightFile = geoTiffFile~substr(1,geoTiffFile~lastPos(".")-1)".hgt"
srtm = .stream~new(heightFile)
srtm~open("write" "replace")
do a=0 to stripCount-1 -- cater for base 1
charinLength = stripByteCounts~substr(1+4*a,4)~reverse~c2d
stripOffset = stripPointers~substr(1+4*a,4)~reverse~c2d+1
strip = geoTiff~charin(stripOffset,charinLength)
-- make strip contents big-endian from little-endian
test = makeBigEndian(strip)
--say test~length
srtm~charout(test)
end
geoTiff~close
srtm~close
exit
::routine makeBigEndian
use arg string
big = .mutablebuffer~new(,string~length)
do i=1 by 2 to string~length
big~append(string~substr(i,2)~reverse)
end
return big
exit
| If you feel inclined to make corrections, suggestions etc.,
please mail me any.
| |
|
All content © Ruurd Idenburg, 2007–,
except where marked otherwise. All rights reserved. This page is primarily for non-commercial use only.
The Idenburg website records no personal information and sets no ‘cookies’.
This site is hosted on a VPS(Virtual Private System) rented from Transip.nl, a Dutch company, falling under Dutch (privacy) laws (I think).
This page updated on by Ruurd Idenburg.
|
|