1 : /******************************************************************************
2 : * $Id: msg_reader_core.cpp 16666 2009-03-28 12:46:49Z rouault $
3 : *
4 : * Project: MSG Native Reader
5 : * Purpose: Base class for reading in the headers of MSG native images
6 : * Author: Frans van den Bergh, fvdbergh@csir.co.za
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2005, Frans van den Bergh <fvdbergh@csir.co.za>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "msg_reader_core.h"
31 : #include "msg_basic_types.h"
32 : #include <stdio.h>
33 : #include <string.h>
34 : #include <math.h>
35 :
36 : #ifdef DEBUG
37 : #ifdef GDAL_SUPPORT
38 : #undef DEBUG
39 : #endif
40 : #endif
41 :
42 : #ifdef GDAL_SUPPORT
43 : #include "cpl_vsi.h"
44 :
45 : CPL_CVSID("$Id: msg_reader_core.cpp 16666 2009-03-28 12:46:49Z rouault $");
46 :
47 : #else
48 : #define VSIFSeek(fp, pos, ref) fseek(fp, pos, ref)
49 : #define VSIFRead(p, bs, nb, fp) fread(p, bs, nb, ref)
50 : #endif
51 :
52 : namespace msg_native_format {
53 :
54 : const Blackbody_lut_type Msg_reader_core::Blackbody_LUT[MSG_NUM_CHANNELS+1] = {
55 : {0,0,0}, // dummy channel
56 : {0,0,0}, // N/A
57 : {0,0,0}, // N/A
58 : {0,0,0}, // N/A
59 : {2569.094, 0.9959, 3.471},
60 : {1598.566, 0.9963, 2.219},
61 : {1362.142, 0.9991, 0.485},
62 : {1149.083, 0.9996, 0.181},
63 : {1034.345, 0.9999, 0.060},
64 : { 930.659, 0.9983, 0.627},
65 : { 839.661, 0.9988, 0.397},
66 : { 752.381, 0.9981, 0.576},
67 : {0,0,0} // N/A
68 : };
69 :
70 :
71 0 : Msg_reader_core::Msg_reader_core(const char* fname) {
72 :
73 0 : FILE* fin = fopen(fname, "rb");
74 0 : if (!fin) {
75 0 : fprintf(stderr, "Could not open file %s\n", fname);
76 0 : return;
77 : }
78 0 : read_metadata_block(fin);
79 : }
80 :
81 0 : Msg_reader_core::Msg_reader_core(FILE* fp) {
82 0 : read_metadata_block(fp);
83 0 : }
84 :
85 :
86 0 : void Msg_reader_core::read_metadata_block(FILE* fin) {
87 0 : _open_success = true;
88 :
89 : unsigned int i;
90 :
91 0 : VSIFRead(&_main_header, sizeof(_main_header), 1, fin);
92 0 : VSIFRead(&_sec_header, sizeof(_sec_header), 1, fin);
93 :
94 : #ifdef DEBUG
95 : // print out all the fields in the header
96 : PH_DATA* hd = (PH_DATA*)&_main_header;
97 : for (int i=0; i < 6; i++) {
98 : to_string(*hd);
99 : printf("[%02d] %s %s", i, hd->name, hd->value);
100 : hd++;
101 : }
102 : PH_DATA_ID* hdi = (PH_DATA_ID*)&_main_header.dataSetIdentification;
103 :
104 : for (i=0; i < 5; i++) {
105 : printf("%s %s %s", hdi->name, hdi->size, hdi->address);
106 : hdi++;
107 : }
108 : hd = (PH_DATA*)(&_main_header.totalFileSize);
109 : for (int i=0; i < 19; i++) {
110 : to_string(*hd);
111 : printf("[%02d] %s %s", i, hd->name, hd->value);
112 : hd++;
113 : }
114 : #endif // DEBUG
115 :
116 : // extract data & header positions
117 :
118 0 : for (i=0; i < 5; i++) {
119 0 : PH_DATA_ID* hdi = (PH_DATA_ID*)&_main_header.dataSetIdentification[i];
120 0 : if (strncmp(hdi->name, "15Header", strlen("15Header")) == 0) {
121 0 : sscanf(hdi->size, "%d", &_f_header_size);
122 0 : sscanf(hdi->address, "%d", &_f_header_offset);
123 : } else
124 0 : if (strncmp(hdi->name, "15Data", strlen("15Data")) == 0) {
125 0 : sscanf(hdi->size, "%d", &_f_data_size);
126 0 : sscanf(hdi->address, "%d", &_f_data_offset);
127 : }
128 : }
129 : #ifdef DEBUG
130 : printf("Data: %d %d\n", _f_data_offset, _f_data_size);
131 : printf("Header: %d %d\n", _f_header_offset, _f_header_size);
132 : #endif // DEBUG
133 :
134 : unsigned int lines;
135 0 : sscanf(_sec_header.northLineSelectedRectangle.value, "%d", &_lines);
136 0 : sscanf(_sec_header.southLineSelectedRectangle.value, "%d", &lines);
137 0 : _line_start = lines;
138 0 : _lines -= lines - 1;
139 :
140 : unsigned int cols;
141 0 : sscanf(_sec_header.westColumnSelectedRectangle.value, "%d", &_columns);
142 0 : sscanf(_sec_header.eastColumnSelectedRectangle.value, "%d", &cols);
143 0 : _col_start = cols;
144 0 : _columns -= cols - 1;
145 :
146 : #ifdef DEBUG
147 : printf("lines = %d, cols = %d\n", _lines, _columns);
148 : #endif // DEBUG
149 :
150 0 : int records_per_line = 0;
151 0 : for (i=0; i < MSG_NUM_CHANNELS; i++) {
152 0 : if (_sec_header.selectedBandIds.value[i] == 'X') {
153 0 : _bands[i] = 1;
154 0 : records_per_line += (i == (MSG_NUM_CHANNELS-1)) ? 3 : 1;
155 : } else {
156 0 : _bands[i] = 0;
157 : }
158 : }
159 :
160 : #ifdef DEBUG
161 : printf("reading a total of %d records per line\n", records_per_line);
162 : #endif // DEBUG
163 :
164 : // extract time fields, assume that SNIT is the correct field:
165 0 : sscanf(_main_header.snit.value + 0, "%04d", &_year);
166 0 : sscanf(_main_header.snit.value + 4, "%02d", &_month);
167 0 : sscanf(_main_header.snit.value + 6, "%02d", &_day);
168 0 : sscanf(_main_header.snit.value + 8, "%02d", &_hour);
169 0 : sscanf(_main_header.snit.value + 10, "%02d", &_minute);
170 :
171 : // read radiometric block
172 : RADIOMETRIC_PROCCESSING_RECORD rad;
173 0 : off_t offset = RADIOMETRICPROCESSING_RECORD_OFFSET + _f_header_offset + sizeof(GP_PK_HEADER) + sizeof(GP_PK_SH1) + 1;
174 0 : VSIFSeek(fin, offset, SEEK_SET);
175 0 : VSIFRead(&rad, sizeof(RADIOMETRIC_PROCCESSING_RECORD), 1, fin);
176 0 : to_native(rad);
177 0 : memcpy((void*)_calibration, (void*)&rad.level1_5ImageCalibration,sizeof(_calibration));
178 :
179 : #ifdef DEBUG
180 : for (unsigned int i=0; i < MSG_NUM_CHANNELS; i++) {
181 : if (_calibration[i].cal_slope < 0 || _calibration[i].cal_slope > 0.4) {
182 : printf("Warning: calibration slope (%lf) out of nominal range. MSG reader probably broken\n", _calibration[i].cal_slope);
183 :
184 : }
185 : if (_calibration[i].cal_offset > 0 || _calibration[i].cal_offset < -20) {
186 : printf("Warning: calibration offset (%lf) out of nominal range. MSG reader probably broken\n", _calibration[i].cal_offset);
187 : }
188 : }
189 : #endif
190 :
191 : // read image description block
192 : IMAGE_DESCRIPTION_RECORD idr;
193 0 : offset = RADIOMETRICPROCESSING_RECORD_OFFSET - IMAGEDESCRIPTION_RECORD_LENGTH + _f_header_offset + sizeof(GP_PK_HEADER) + sizeof(GP_PK_SH1) + 1;
194 0 : VSIFSeek(fin, offset, SEEK_SET);
195 0 : VSIFRead(&idr, sizeof(IMAGE_DESCRIPTION_RECORD), 1, fin);
196 0 : to_native(idr);
197 0 : _line_dir_step = idr.referencegrid_visir.lineDirGridStep;
198 0 : _col_dir_step = idr.referencegrid_visir.columnDirGridStep;
199 :
200 :
201 : // Rather convoluted, but this code is required to compute the real data block sizes
202 : // It does this by reading in the first line of every band, to get to the packet size field
203 : GP_PK_HEADER gp_header;
204 : GP_PK_SH1 sub_header;
205 : SUB_VISIRLINE visir_line;
206 :
207 0 : VSIFSeek(fin, _f_data_offset, SEEK_SET);
208 :
209 0 : _hrv_packet_size = 0;
210 0 : _interline_spacing = 0;
211 0 : visir_line.channelId = 0;
212 :
213 : int scanned_bands[MSG_NUM_CHANNELS];
214 0 : int band_count = 0;
215 0 : for (i=0; i < MSG_NUM_CHANNELS; i++) {
216 0 : scanned_bands[i] = _bands[i];
217 0 : band_count += _bands[i];
218 : }
219 :
220 0 : do {
221 0 : VSIFRead(&gp_header, sizeof(GP_PK_HEADER), 1, fin);
222 0 : VSIFRead(&sub_header, sizeof(GP_PK_SH1), 1, fin);
223 0 : VSIFRead(&visir_line, sizeof(SUB_VISIRLINE), 1, fin);
224 0 : to_native(visir_line);
225 0 : to_native(gp_header);
226 :
227 : // skip over the actual line data
228 : VSIFSeek(fin,
229 : gp_header.packetLength - (sizeof(GP_PK_SH1) + sizeof(SUB_VISIRLINE) - 1),
230 : SEEK_CUR
231 0 : );
232 :
233 0 : if (visir_line.channelId == 0 || visir_line.channelId > MSG_NUM_CHANNELS) {
234 0 : _open_success = false;
235 0 : break;
236 : }
237 :
238 0 : if (scanned_bands[visir_line.channelId - 1]) {
239 0 : scanned_bands[visir_line.channelId - 1] = 0;
240 0 : band_count--;
241 :
242 0 : if (visir_line.channelId != 12) { // not the HRV channel
243 0 : _visir_bytes_per_line = gp_header.packetLength - (sizeof(GP_PK_SH1) + sizeof(SUB_VISIRLINE) - 1);
244 0 : _visir_packet_size = gp_header.packetLength + sizeof(GP_PK_HEADER) + 1;
245 0 : _interline_spacing += _visir_packet_size;
246 : } else {
247 0 : _hrv_bytes_per_line = gp_header.packetLength - (sizeof(GP_PK_SH1) + sizeof(SUB_VISIRLINE) - 1);
248 0 : _hrv_packet_size = gp_header.packetLength + sizeof(GP_PK_HEADER) + 1;
249 0 : _interline_spacing += 3*_hrv_packet_size;
250 0 : VSIFSeek(fin, 2*gp_header.packetLength, SEEK_CUR );
251 : }
252 : }
253 : } while (band_count > 0);
254 0 : }
255 :
256 : #ifndef GDAL_SUPPORT
257 :
258 : int Msg_reader_core::_chan_to_idx(Msg_channel_names channel) {
259 : unsigned int idx = 0;
260 : while (idx < MSG_NUM_CHANNELS) {
261 : if ( (1 << (idx+1)) == (int)channel ) {
262 : return idx;
263 : }
264 : idx++;
265 : }
266 : return 0;
267 : }
268 :
269 : void Msg_reader_core::get_pixel_geo_coordinates(unsigned int line, unsigned int column, double& longitude, double& latitude) {
270 : Conversions::convert_pixel_to_geo((unsigned int)(line + _line_start), (unsigned int)(column + _col_start), longitude, latitude);
271 : }
272 :
273 : void Msg_reader_core::get_pixel_geo_coordinates(double line, double column, double& longitude, double& latitude) {
274 : Conversions::convert_pixel_to_geo(line + _line_start, column + _col_start, longitude, latitude);
275 : }
276 :
277 : double Msg_reader_core::compute_pixel_area_sqkm(double line, double column) {
278 : return Conversions::compute_pixel_area_sqkm(line + _line_start, column + _col_start);
279 : }
280 :
281 : #endif // GDAL_SUPPORT
282 :
283 : } // namespace msg_native_format
284 :
|
