-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathread_hrsc_dtm.m
174 lines (142 loc) · 7.04 KB
/
read_hrsc_dtm.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
function [id, core, corelonvec, corelatvec]=read_hrsc_dtm(cubfile)
% This function reads in DTMs derived from the stereo channels
% on the High Resolution Stereo Camera aboard Mars Express.
%
%function [id, core, corelonvec, corelatvec]=read_hrsc_dtm(cubfile)
%
%Input:
% -'cubfile' is the HRSC DT4 or DA4 file to extract the data
% from. These files usually have a '.IMG' extension if downloaded
% from the Planetary Science Archive at ESA.
%Output:
% -'id' is the identifier of the image being processed
% -'core' is the actual dtm data
% -'corelonvec' is the longitude vector for the 'core' data
% -'corelatvec' is the latitude vector for the 'core' data
%
% Copyright 2008-2016 Elliot Sefton-Nash
%---------FIRST OPEN FILE AS TEXT TO READ LABEL INFO---------
%
% Define the set of image parameters found in the label/header
% which are needed for Matlab reading and processing of image data.
%
% The structure of a the keys2get array addressing is as follows:
% 'Parameter Parent1 Parent2 Parent3'
%
keys2get = {'PRODUCT_ID'; % 1 Image identifier
'RECORD_BYTES'; % 2 Number of bytes per record
'^IMAGE'; % 3 Number of the record that the image data starts at
'INTERCHANGE_FORMAT'; % 4
'MINIMUM_LATITUDE IMAGE_MAP_PROJECTION'; % 5
'MAXIMUM_LATITUDE IMAGE_MAP_PROJECTION'; % 6
'WESTERNMOST_LONGITUDE IMAGE_MAP_PROJECTION'; % 7
'EASTERNMOST_LONGITUDE IMAGE_MAP_PROJECTION'; % 8
'LINES IMAGE'; % 9
'LINE_SAMPLES IMAGE'; % 10
'SAMPLE_TYPE IMAGE'; % 11 This should be MSB_INTEGER for HRSC DTMs
'SAMPLE_BITS IMAGE'; % 12 This should be 16
'MEX:DTM_MISSING_DN MEX:DTM'}; % 13 The DN for null pixels.
% Create the keyword data structure with attributes addresses and
% values.
sizekeys2get = size(keys2get);
for i = 1:sizekeys2get(1)
keys(i).address = keys2get(i);
keys(i).value = 0;
end
inheader = 1;
infooter = 0;
level = 0; % Holds the level of indentation that the cursor is at. Allows us to address the keywords based on thier 'address' held in keys2get.
footerstartbyte = -1; % To begin with, assume there is no EOL. Some files may not have them.
address = '';
fid = fopen(cubfile, 'rt');
% If we find an 'End' in the header, we proceed to the start of the EOL,
% whose location may have been read. If it has not then we close the file.
% Sometimes in the footer there is an 'End' before the end of the EOL. So
% we keep reading regardless until the EOF if we are in the EOL.
while ~feof(fid) && (inheader || infooter)
hline = fgetl(fid);
% Separate hline into keyword and value either side of the '=' delimiter.
[keyword,value] = strtok(hline,'=');
if ~isempty(keyword)
% Remove spaces and the '=' assignment operator.
keyword = sscanf(keyword,'%s');
value = strtrim(strtok(value, '='));
% Test if the keyword is Group OR Object OR End_Group OR End_Object. If so we need to adjust
% the parent tree.
switch upper(keyword)
case {'OBJECT', 'GROUP'}
level = level + 1;
address = [value, ' ', address]; %#ok<AGROW> % Push the new parent at the front of the address.
case {'END_OBJECT', 'END_GROUP'}
level = level - 1;
[reject, address] = strtok(address, ' '); % Take off the front value from the address string
address = strtrim(address);
case 'END'
% If there is an END is the header, it means we should
% proceed to the footer, ONLY if footerstartbyte != -1.
% Else we should quit this while loop, done by setting
% both infooter and inheader = 0.
if inheader && (footerstartbyte > -1)
fseek(fid, footerstartbyte, 'bof');
infooter = 1;
end
inheader = 0;
otherwise
% Now perform an action based on the current address with
% the current keyword added to it.
fulladdress = strtrim([keyword, ' ', address]);
if strcmp(upper(fulladdress), 'STARTBYTE ORIGINALLABEL');
footerstartbyte = str2double(value);
else
for i = 1:sizekeys2get
if strcmpi(keys(i).address, fulladdress)
keys(i).value = value;
end
end
end
fulladdress = '';
end
end
end
fclose(fid);
% Should test here to see if any of the keys are empty and return an error
% if so, as we cannot load the image unless they are all obtained from the file.
%------------------------- NOW READ CORE IMAGE DATA -----------------------
% IMG FILE PARAMETERS
id = keys(1).value;
record_bytes = str2double(keys(2).value); % Label bytes
image_pointer = str2double(keys(3).value); % Record pointer to the start of the image data
samples = str2double(keys(10).value); % Image samples
lines = str2double(keys(9).value); % Image lines
% The start byte of the data is (^IMAGE - 1) * RECORD_BYTES
data_start_byte = (image_pointer - 1)*record_bytes;
% Find out the precision based on passing, in this case, the value
% associated with the PDS keyword 'SAMPLE_BITS'.
[precision, pixel_bytes] = get_precision(keys(12).value);
% Find out what byte-ordering the data is. The keyword SAMPLE_TYPE tells us
% this. For HRSC it's usually MSB_INTEGER, which is big endian.
endian = get_endian(keys(11).value);
% Open the file as binary read-only.
fid = fopen(cubfile, 'r', endian);
% Skip the header.
fseek(fid, data_start_byte, 'bof');
% Read the image data (fread fills column-wise, so must transpose later)
core = fread(fid, [samples, lines], precision);
% Rotate by 90 degrees.
core = core';
fclose(fid);
% Set all pixels equal to the null value equal to zero.
null_dn = str2double(keys(13).value);
core(find(core == null_dn)) = 0.0;
%-----LAT LON-----
maxlat = str2double(keys(6).value); % MAXIMUM_LATITUDE IMAGE_MAP_PROJECTION
minlat = str2double(keys(5).value); % MINIMUM_LATITUDE IMAGE_MAP_PROJECTION
maxlon = str2double(keys(8).value); % EASTERNMOST_LONGITUDE IMAGE_MAP_PROJECTION
minlon = str2double(keys(7).value); % WESTERNMOST_LONGITUDE IMAGE_MAP_PROJECTION
% Extract degrees per pixel in lat and lon.
dlat = (maxlat-minlat)/(lines-1);
dlon = (maxlon-minlon)/(samples-1);
% Create the lat/lon vectors.
corelatvec = maxlat:-1*dlat:minlat;
corelonvec = minlon:dlon:maxlon;
end