1: /* ---------------------------------------------------------------- */
2: /* Transforms a GeoTiff file as produced by ASTER GDEM, see: */
3: /* */
4: /* http://gdem.ersdac.jspacesystems.or.jp */
5: /* */
6: /* to 'hgt' format as used by SRTM, see: */
7: /* */
8: /* https://lpdaac.usgs.gov/products/measures_products_table/srtmgl3 */
9: /* */
10: /* Foa a nice explanation of the (geo)tiff format see: */
11: /* */
12: /* http://www.awaresystems.be/imaging/tiff.html */
13: /* */
14: /* ---------------------------------------------------------------- */
15: /* */
16: /* Originally by Ruurd J. Idenburg */
17: /* */
18: /* No copyright, no licence, no guarantees or warrantees, be it */
19: /* explicit, implicit or whatever. Usage is totally and completely */
20: /* at the users own risk, the author shall not be liable for any */
21: /* damages whatsoever, for any reason whatsoever. */
22: /* */
23: /* Please keep this comment block intact when modifying this code */
24: /* and add a note with date and a description. */
25: /* */
26: /* ---------------------------------------------------------------- */
27: /* Parameter(s): */
28: /* */
29: /* file - the (geo)tiff path and filename to be converted */
30: /* */
31: /* */
32: /* Result: */
33: /* */
34: /* A file on the same path and filename as the input tiff file */
35: /* but with file extension 'hgt' and in the SRTM format */
36: /* */
37: /* ---------------------------------------------------------------- */
38: /* 2014/02/01 - Initial version */
39: /* ---------------------------------------------------------------- */
40:
41: use arg geoTiffFile
42: geoTiff = .stream~new(geoTiffFile)
43: geoTiff~open("read")
44: /* ---------------------------------------------------------------- */
45: /* The Image File Directory(IFD) is always at offset 4 in the file. */
46: /* ---------------------------------------------------------------- */
47: ifdOffset = (geoTiff~charin(5,4)~reverse~c2d)+1 -- cater for base 1
48: say "Found the IFD at offset" ifdOffset "('"ifdOffset~d2x(8)"'X)."
49: /* ---------------------------------------------------------------- */
50: /* Get the number of tags present in this IFD. */
51: /* ---------------------------------------------------------------- */
52: tagCount = geoTiff~charin(ifdOffset,2)~reverse~c2d
53: say "IFD contains" tagCount "tags."
54: /* ---------------------------------------------------------------- */
55: /* The tags follow directly after the count and each 12 bytes long. */
56: /* ---------------------------------------------------------------- */
57: tags = geoTiff~charin(ifdOffset+2,tagCount*12)
58: nextIfdOffset = geoTiff~charin(ifdOffset+2+tagCount*12,4)~reverse~c2d
59: -- say nextIfdOffset -- for now assume no following IFD's
60: tagDir = .directory~new
61: /* ---------------------------------------------------------------- */
62: /* Each tag consists of Key, Datatype, Valuecount and Value(s). */
63: /* If the values don't fit in 4 bytes, the value is an offset to */
64: /* where the values can be found. */
65: /* ---------------------------------------------------------------- */
66: do t=0 to tagCount-1 --cater for base 1
67: aTag = tags~substr(1+12*t,12)
68: tagKey = aTag~substr(1,2)~reverse~c2x
69: dataType = aTag~substr(3,2)~reverse~c2x
70: valueCount = aTag~substr(5,4)~reverse~c2d
71: valuePointer = aTag~substr(9,4)~reverse~c2d -- can be a value only
72: tag = .directory~new
73: tagDir~put(dataType valueCount valuePointer,tagKey)
74: end
75: columns = tagDir~0100~word(3)
76: say "The elevation information is" columns "columns wide."
77: rows = tagDir~0101~word(3)
78: say "The elevation information is" rows "rows deep."
79: --rowsPerStrip = tagDir~0116~word(3) -- not needed for further processing
80: --say rowsPerStrip
81: stripByteCountValues = tagDir~0117~word(2)
82: say "There are" stripByteCountValues "byte counts for the row strips to read."
83: stripByteCountsOffset = tagDir~0117~word(3)+1 -- cater for base 1
84: --say stripByteCountsOffset
85: stripByteCounts = geoTiff~charin(stripByteCountsOffset,stripByteCountValues*4)
86: --say stripByteCounts~c2x
87: stripCount = tagDir~0111~word(2)
88: say "The strip count" stripCount "should be equal to previous value."
89: stripOffsetsPointer = tagDir~0111~word(3)+1 -- cater for base 1
90: --say stripOffsetsPointer
91: stripPointers = geoTiff~charin(stripOffsetsPointer,stripCount*4)
92: heightFile = geoTiffFile~substr(1,geoTiffFile~lastPos(".")-1)".hgt"
93: srtm = .stream~new(heightFile)
94: srtm~open("write" "replace")
95: do a=0 to stripCount-1 -- cater for base 1
96: charinLength = stripByteCounts~substr(1+4*a,4)~reverse~c2d
97: stripOffset = stripPointers~substr(1+4*a,4)~reverse~c2d+1
98: strip = geoTiff~charin(stripOffset,charinLength)
99: -- make strip contents big-endian from little-endian
100: test = makeBigEndian(strip)
101: --say test~length
102: srtm~charout(test)
103: end
104: geoTiff~close
105: srtm~close
106: exit
107:
108: ::routine makeBigEndian
109: use arg string
110: big = .mutablebuffer~new(,string~length)
111: do i=1 by 2 to string~length
112: big~append(string~substr(i,2)~reverse)
113: end
114: return big
115: exit
116:
117:
