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