1 : /******************************************************************************
2 : * $Id: msgndataset.cpp 17664 2009-09-21 21:16:45Z rouault $
3 : *
4 : * Project: MSG Native Reader
5 : * Purpose: All code for EUMETSAT Archive format reader
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 "gdal_priv.h"
31 : #include "ogr_spatialref.h"
32 :
33 : #include "msg_reader_core.h"
34 : using namespace msg_native_format;
35 :
36 : CPL_CVSID("$Id: msgndataset.cpp 17664 2009-09-21 21:16:45Z rouault $");
37 :
38 : CPL_C_START
39 : void GDALRegister_MSGN(void);
40 : CPL_C_END
41 :
42 : typedef enum {
43 : MODE_VISIR, // Visible and Infrared bands (1 through 11) in 10-bit raw mode
44 : MODE_HRV, // Pan band (band 11) only, in 10-bit raw mode
45 : MODE_RAD // Black-body temperature (K) for thermal bands only (4-10), 64-bit float
46 : } open_mode_type;
47 :
48 : class MSGNRasterBand;
49 :
50 : /************************************************************************/
51 : /* ==================================================================== */
52 : /* MSGNDataset */
53 : /* ==================================================================== */
54 : /************************************************************************/
55 :
56 : class MSGNDataset : public GDALDataset
57 : {
58 : friend class MSGNRasterBand;
59 :
60 : FILE *fp;
61 : GByte abyHeader[1012];
62 :
63 : Msg_reader_core* msg_reader_core;
64 : double adfGeoTransform[6];
65 : char *pszProjection;
66 :
67 : public:
68 : MSGNDataset();
69 : ~MSGNDataset();
70 :
71 : static GDALDataset *Open( GDALOpenInfo * );
72 :
73 : CPLErr GetGeoTransform( double * padfTransform );
74 : const char *GetProjectionRef();
75 :
76 : };
77 :
78 : /************************************************************************/
79 : /* ==================================================================== */
80 : /* MSGNRasterBand */
81 : /* ==================================================================== */
82 : /************************************************************************/
83 :
84 : class MSGNRasterBand : public GDALRasterBand
85 0 : {
86 : friend class MSGNDataset;
87 :
88 : unsigned int packet_size;
89 : unsigned int bytes_per_line;
90 : unsigned int interline_spacing;
91 : unsigned int orig_band_no; // The name of the band
92 : unsigned int band_in_file; // The effective index of the band in the file
93 : open_mode_type open_mode;
94 :
95 0 : double GetNoDataValue (int *pbSuccess=NULL) {
96 0 : if (pbSuccess) {
97 0 : *pbSuccess = 1;
98 : }
99 0 : return MSGN_NODATA_VALUE;
100 : }
101 :
102 : double MSGN_NODATA_VALUE;
103 :
104 : char band_description[30];
105 :
106 : public:
107 :
108 : MSGNRasterBand( MSGNDataset *, int , open_mode_type mode, int orig_band_no, int band_in_file);
109 :
110 : virtual CPLErr IReadBlock( int, int, void * );
111 : virtual double GetMinimum( int *pbSuccess = NULL );
112 : virtual double GetMaximum(int *pbSuccess = NULL );
113 0 : virtual const char* GetDescription() const { return band_description; }
114 : };
115 :
116 :
117 : /************************************************************************/
118 : /* MSGNRasterBand() */
119 : /************************************************************************/
120 :
121 0 : MSGNRasterBand::MSGNRasterBand( MSGNDataset *poDS, int nBand , open_mode_type mode, int orig_band_no, int band_in_file)
122 :
123 : {
124 0 : this->poDS = poDS;
125 0 : this->nBand = nBand; // GDAL's band number, i.e. always starts at 1
126 0 : this->open_mode = mode;
127 0 : this->orig_band_no = orig_band_no;
128 0 : this->band_in_file = band_in_file;
129 :
130 0 : sprintf(band_description, "band %02d", orig_band_no);
131 :
132 0 : if (mode != MODE_RAD) {
133 0 : eDataType = GDT_UInt16;
134 0 : MSGN_NODATA_VALUE = 0;
135 : } else {
136 0 : eDataType = GDT_Float64;
137 0 : MSGN_NODATA_VALUE = -1000;
138 : }
139 :
140 0 : nBlockXSize = poDS->GetRasterXSize();
141 0 : nBlockYSize = 1;
142 :
143 0 : if (mode != MODE_HRV) {
144 0 : packet_size = poDS->msg_reader_core->get_visir_packet_size();
145 0 : bytes_per_line = poDS->msg_reader_core->get_visir_bytes_per_line();
146 : } else {
147 0 : packet_size = poDS->msg_reader_core->get_hrv_packet_size();
148 0 : bytes_per_line = poDS->msg_reader_core->get_hrv_bytes_per_line();
149 : }
150 :
151 0 : interline_spacing = poDS->msg_reader_core->get_interline_spacing();
152 0 : }
153 :
154 : /************************************************************************/
155 : /* IReadBlock() */
156 : /************************************************************************/
157 :
158 0 : CPLErr MSGNRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
159 : void * pImage )
160 :
161 : {
162 0 : MSGNDataset *poGDS = (MSGNDataset *) poDS;
163 :
164 : // invert y position
165 0 : int i_nBlockYOff = poDS->GetRasterYSize() - 1 - nBlockYOff;
166 :
167 : char *pszRecord;
168 :
169 0 : unsigned int data_length = bytes_per_line + sizeof(SUB_VISIRLINE);
170 0 : unsigned int data_offset = 0;
171 :
172 0 : if (open_mode != MODE_HRV) {
173 : data_offset = poGDS->msg_reader_core->get_f_data_offset() +
174 : interline_spacing*i_nBlockYOff + (band_in_file-1)*packet_size +
175 0 : (packet_size - data_length);
176 : } else {
177 : data_offset = poGDS->msg_reader_core->get_f_data_offset() +
178 : interline_spacing*(int(i_nBlockYOff/3) + 1) -
179 0 : packet_size*(3 - (i_nBlockYOff % 3)) + (packet_size - data_length);
180 : }
181 :
182 0 : VSIFSeek( poGDS->fp, data_offset, SEEK_SET );
183 :
184 0 : pszRecord = (char *) CPLMalloc(data_length);
185 0 : size_t nread = VSIFRead( pszRecord, 1, data_length, poGDS->fp );
186 :
187 0 : SUB_VISIRLINE* p = (SUB_VISIRLINE*) pszRecord;
188 0 : to_native(*p);
189 :
190 0 : if (p->lineValidity != 1) {
191 0 : for (int c=0; c < nBlockXSize; c++) {
192 0 : if (open_mode != MODE_RAD) {
193 0 : ((GUInt16 *)pImage)[c] = (GUInt16)MSGN_NODATA_VALUE;
194 : } else {
195 0 : ((double *)pImage)[c] = MSGN_NODATA_VALUE;
196 : }
197 : }
198 : }
199 :
200 0 : if ( nread != data_length ||
201 : ( open_mode != MODE_HRV && (p->lineNumberInVisirGrid - poGDS->msg_reader_core->get_line_start()) != (unsigned int)i_nBlockYOff )
202 : ) { // no sophisticated checking for HRV at the moment
203 0 : CPLFree( pszRecord );
204 :
205 0 : CPLError( CE_Failure, CPLE_AppDefined, "MSGN Scanline corrupt." );
206 :
207 0 : return CE_Failure;
208 : }
209 :
210 : // unpack the 10-bit values into 16-bit unsigned short ints
211 : unsigned char *cptr = (unsigned char*)pszRecord +
212 0 : (data_length - bytes_per_line);
213 0 : int bitsLeft = 8;
214 0 : unsigned short value = 0;
215 :
216 0 : if (open_mode != MODE_RAD) {
217 0 : for (int c=0; c < nBlockXSize; c++) {
218 0 : value = 0;
219 0 : for (int bit=0; bit < 10; bit++) {
220 0 : value <<= 1;
221 0 : if (*cptr & 128) {
222 0 : value |= 1;
223 : }
224 0 : *cptr <<= 1;
225 0 : bitsLeft--;
226 0 : if (bitsLeft == 0) {
227 0 : cptr++;
228 0 : bitsLeft = 8;
229 : }
230 : }
231 0 : ((GUInt16 *)pImage)[nBlockXSize-1 - c] = value;
232 : }
233 : } else {
234 : // radiance mode
235 0 : for (int c=0; c < nBlockXSize; c++) {
236 0 : value = 0;
237 0 : for (int bit=0; bit < 10; bit++) {
238 0 : value <<= 1;
239 0 : if (*cptr & 128) {
240 0 : value |= 1;
241 : }
242 0 : *cptr <<= 1;
243 0 : bitsLeft--;
244 0 : if (bitsLeft == 0) {
245 0 : cptr++;
246 0 : bitsLeft = 8;
247 : }
248 : }
249 0 : double dvalue = double(value);
250 0 : double bbvalue = MSGN_NODATA_VALUE;
251 0 : bbvalue = dvalue * poGDS->msg_reader_core->get_calibration_parameters()[orig_band_no-1].cal_slope +
252 0 : poGDS->msg_reader_core->get_calibration_parameters()[orig_band_no-1].cal_offset;
253 :
254 0 : ((double *)pImage)[nBlockXSize-1 -c] = bbvalue;
255 : }
256 : }
257 0 : CPLFree( pszRecord );
258 0 : return CE_None;
259 : }
260 :
261 : /************************************************************************/
262 : /* GetMinimum() */
263 : /************************************************************************/
264 0 : double MSGNRasterBand::GetMinimum( int *pbSuccess ) {
265 0 : if (pbSuccess) {
266 0 : *pbSuccess = 1;
267 : }
268 0 : return open_mode != MODE_RAD ? 1 : GDALRasterBand::GetMinimum(pbSuccess);
269 : }
270 :
271 : /************************************************************************/
272 : /* GetMaximum() */
273 : /************************************************************************/
274 0 : double MSGNRasterBand::GetMaximum(int *pbSuccess ) {
275 0 : if (pbSuccess) {
276 0 : *pbSuccess = 1;
277 : }
278 0 : return open_mode != MODE_RAD ? 1023 : GDALRasterBand::GetMaximum(pbSuccess);
279 : }
280 :
281 : /************************************************************************/
282 : /* ==================================================================== */
283 : /* MSGNDataset */
284 : /* ==================================================================== */
285 : /************************************************************************/
286 :
287 0 : MSGNDataset::MSGNDataset() {
288 0 : pszProjection = CPLStrdup("");
289 0 : msg_reader_core = 0;
290 0 : }
291 :
292 : /************************************************************************/
293 : /* ~MSGNDataset() */
294 : /************************************************************************/
295 :
296 0 : MSGNDataset::~MSGNDataset()
297 :
298 : {
299 0 : if( fp != NULL )
300 0 : VSIFClose( fp );
301 :
302 0 : if (msg_reader_core) {
303 0 : delete msg_reader_core;
304 : }
305 :
306 0 : CPLFree(pszProjection);
307 0 : }
308 :
309 : /************************************************************************/
310 : /* GetGeoTransform() */
311 : /************************************************************************/
312 :
313 0 : CPLErr MSGNDataset::GetGeoTransform( double * padfTransform )
314 :
315 : {
316 :
317 0 : for (int i=0; i < 6; i++) {
318 0 : padfTransform[i] = adfGeoTransform[i];
319 : }
320 :
321 0 : return CE_None;
322 : }
323 :
324 : /************************************************************************/
325 : /* GetProjectionRef() */
326 : /************************************************************************/
327 :
328 0 : const char *MSGNDataset::GetProjectionRef()
329 :
330 : {
331 0 : return ( pszProjection );
332 : }
333 :
334 : /************************************************************************/
335 : /* Open() */
336 : /************************************************************************/
337 :
338 8884 : GDALDataset *MSGNDataset::Open( GDALOpenInfo * poOpenInfo )
339 :
340 : {
341 8884 : open_mode_type open_mode = MODE_VISIR;
342 8884 : GDALOpenInfo* open_info = poOpenInfo;
343 :
344 8884 : if (!poOpenInfo->bStatOK) {
345 8245 : if ( EQUALN(poOpenInfo->pszFilename, "HRV:", 4) ) {
346 0 : open_info = new GDALOpenInfo(&poOpenInfo->pszFilename[4], poOpenInfo->eAccess);
347 0 : open_mode = MODE_HRV;
348 : } else
349 8245 : if ( EQUALN(poOpenInfo->pszFilename, "RAD:", 4 ) ) {
350 0 : open_info = new GDALOpenInfo(&poOpenInfo->pszFilename[4], poOpenInfo->eAccess);
351 0 : open_mode = MODE_RAD;
352 : }
353 : }
354 :
355 : /* -------------------------------------------------------------------- */
356 : /* Before trying MSGNOpen() we first verify that there is at */
357 : /* least one "\n#keyword" type signature in the first chunk of */
358 : /* the file. */
359 : /* -------------------------------------------------------------------- */
360 8884 : if( open_info->fp == NULL || open_info->nHeaderBytes < 50 )
361 8463 : return NULL;
362 :
363 : /* check if this is a "NATIVE" MSG format image */
364 421 : if( !EQUALN((char *)open_info->pabyHeader,
365 : "FormatName : NATIVE", 36) )
366 : {
367 421 : return NULL;
368 : }
369 :
370 : /* -------------------------------------------------------------------- */
371 : /* Confirm the requested access is supported. */
372 : /* -------------------------------------------------------------------- */
373 0 : if( poOpenInfo->eAccess == GA_Update )
374 : {
375 : CPLError( CE_Failure, CPLE_NotSupported,
376 : "The MSGN driver does not support update access to existing"
377 0 : " datasets.\n" );
378 0 : return NULL;
379 : }
380 :
381 : /* -------------------------------------------------------------------- */
382 : /* Create a corresponding GDALDataset. */
383 : /* -------------------------------------------------------------------- */
384 : MSGNDataset *poDS;
385 :
386 0 : poDS = new MSGNDataset();
387 :
388 0 : poDS->fp = open_info->fp;
389 0 : open_info->fp = NULL;
390 :
391 : /* -------------------------------------------------------------------- */
392 : /* Read the header. */
393 : /* -------------------------------------------------------------------- */
394 : // first reset the file pointer, then hand over to the msg_reader_core
395 0 : VSIFSeek( poDS->fp, 0, SEEK_SET );
396 :
397 0 : poDS->msg_reader_core = new Msg_reader_core(poDS->fp);
398 :
399 0 : if (!poDS->msg_reader_core->get_open_success()) {
400 0 : delete poDS;
401 0 : return NULL;
402 : }
403 :
404 0 : poDS->nRasterXSize = poDS->msg_reader_core->get_columns();
405 0 : poDS->nRasterYSize = poDS->msg_reader_core->get_lines();
406 :
407 0 : if (open_mode == MODE_HRV) {
408 0 : poDS->nRasterXSize *= 3;
409 0 : poDS->nRasterYSize *= 3;
410 : }
411 :
412 :
413 : /* -------------------------------------------------------------------- */
414 : /* Create band information objects. */
415 : /* -------------------------------------------------------------------- */
416 : unsigned int i;
417 0 : unsigned int band_count = 1;
418 0 : unsigned int missing_band_count = 0;
419 0 : unsigned char* bands = poDS->msg_reader_core->get_band_map();
420 : unsigned char band_map[MSG_NUM_CHANNELS+1]; // map GDAL band numbers to MSG channels
421 0 : for (i=0; i < MSG_NUM_CHANNELS; i++) {
422 0 : if (bands[i]) {
423 0 : bool ok_to_add = false;
424 0 : switch (open_mode) {
425 : case MODE_VISIR:
426 0 : ok_to_add = i < MSG_NUM_CHANNELS - 1;
427 0 : break;
428 : case MODE_RAD:
429 0 : ok_to_add = (i <= 2) || (Msg_reader_core::Blackbody_LUT[i+1].B != 0);
430 0 : break;
431 : case MODE_HRV:
432 0 : ok_to_add = i == MSG_NUM_CHANNELS - 1;
433 : break;
434 : }
435 0 : if (ok_to_add) {
436 0 : poDS->SetBand( band_count, new MSGNRasterBand( poDS, band_count, open_mode, i+1, i+1 - missing_band_count));
437 0 : band_map[band_count] = i+1;
438 0 : band_count++;
439 : }
440 : } else {
441 0 : missing_band_count++;
442 : }
443 : }
444 :
445 : double pixel_gsd_x;
446 : double pixel_gsd_y;
447 : double origin_x;
448 : double origin_y;
449 :
450 0 : if (open_mode != MODE_HRV) {
451 0 : pixel_gsd_x = 1000 * poDS->msg_reader_core->get_col_dir_step(); // convert from km to m
452 0 : pixel_gsd_y = 1000 * poDS->msg_reader_core->get_line_dir_step(); // convert from km to m
453 0 : origin_x = -pixel_gsd_x * (-(Conversions::nlines / 2.0) + poDS->msg_reader_core->get_col_start());
454 0 : origin_y = -pixel_gsd_y * ((Conversions::nlines / 2.0) - poDS->msg_reader_core->get_line_start());
455 : } else {
456 0 : pixel_gsd_x = 1000 * poDS->msg_reader_core->get_col_dir_step() / 3.0; // convert from km to m, approximate for HRV
457 0 : pixel_gsd_y = 1000 * poDS->msg_reader_core->get_line_dir_step() / 3.0; // convert from km to m, approximate for HRV
458 0 : origin_x = -pixel_gsd_x * (-(3*Conversions::nlines / 2.0) + 3*poDS->msg_reader_core->get_col_start());
459 0 : origin_y = -pixel_gsd_y * ((3*Conversions::nlines / 2.0) - 3*poDS->msg_reader_core->get_line_start());
460 : }
461 :
462 0 : poDS->adfGeoTransform[0] = origin_x;
463 0 : poDS->adfGeoTransform[1] = pixel_gsd_x;
464 0 : poDS->adfGeoTransform[2] = 0.0;
465 :
466 0 : poDS->adfGeoTransform[3] = origin_y;
467 0 : poDS->adfGeoTransform[4] = 0.0;
468 0 : poDS->adfGeoTransform[5] = -pixel_gsd_y;
469 :
470 0 : OGRSpatialReference oSRS;
471 :
472 0 : oSRS.SetProjCS("Geostationary projection (MSG)");
473 0 : oSRS.SetGEOS( 0, 35785831, 0, 0 );
474 : oSRS.SetGeogCS(
475 : "MSG Ellipsoid",
476 : "MSG_DATUM",
477 : "MSG_SPHEROID",
478 : Conversions::rpol * 1000.0,
479 : 1 / ( 1 - Conversions::rpol/Conversions::req)
480 0 : );
481 :
482 0 : oSRS.exportToWkt( &(poDS->pszProjection) );
483 :
484 0 : CALIBRATION* cal = poDS->msg_reader_core->get_calibration_parameters();
485 : char tagname[30];
486 : char field[300];
487 :
488 0 : poDS->SetMetadataItem("Radiometric parameters format", "offset slope");
489 0 : for (i=1; i < band_count; i++) {
490 0 : sprintf(tagname, "ch%02d_cal", band_map[i]);
491 0 : sprintf(field, "%.12e %.12e", cal[band_map[i]-1].cal_offset, cal[band_map[i]-1].cal_slope);
492 0 : poDS->SetMetadataItem(tagname, field);
493 : }
494 :
495 : sprintf(field, "%04d%02d%02d/%02d:%02d",
496 : poDS->msg_reader_core->get_year(),
497 : poDS->msg_reader_core->get_month(),
498 : poDS->msg_reader_core->get_day(),
499 : poDS->msg_reader_core->get_hour(),
500 : poDS->msg_reader_core->get_minute()
501 0 : );
502 0 : poDS->SetMetadataItem("Date/Time", field);
503 :
504 : sprintf(field, "%d %d",
505 : poDS->msg_reader_core->get_line_start(),
506 : poDS->msg_reader_core->get_col_start()
507 0 : );
508 0 : poDS->SetMetadataItem("Origin", field);
509 :
510 :
511 0 : if (open_info != poOpenInfo) {
512 0 : delete open_info;
513 : }
514 :
515 0 : return( poDS );
516 : }
517 :
518 : /************************************************************************/
519 : /* GDALRegister_MSGN() */
520 : /************************************************************************/
521 :
522 338 : void GDALRegister_MSGN()
523 :
524 : {
525 : GDALDriver *poDriver;
526 :
527 338 : if( GDALGetDriverByName( "MSGN" ) == NULL )
528 : {
529 336 : poDriver = new GDALDriver();
530 :
531 336 : poDriver->SetDescription( "MSGN" );
532 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
533 336 : "EUMETSAT Archive native (.nat)" );
534 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
535 336 : "frmt_msgn.html" );
536 336 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "nat" );
537 :
538 336 : poDriver->pfnOpen = MSGNDataset::Open;
539 :
540 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
541 : }
542 338 : }
|
