Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: MSG Native Reader
4 : * Purpose: All code for EUMETSAT Archive format reader
5 : * Author: Frans van den Bergh, fvdbergh@csir.co.za
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2005, Frans van den Bergh <fvdbergh@csir.co.za>
9 : * Copyright (c) 2008-2009, Even Rouault <even dot rouault at spatialys.com>
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 : #include "cpl_port.h"
30 : #include "cpl_error.h"
31 :
32 : #include "gdal_frmts.h"
33 : #include "gdal_priv.h"
34 : #include "ogr_spatialref.h"
35 :
36 : #include "msg_reader_core.h"
37 :
38 : #include <algorithm>
39 :
40 : using namespace msg_native_format;
41 :
42 : typedef enum
43 : {
44 : MODE_VISIR, // Visible and Infrared bands (1 through 11) in 10-bit raw mode
45 : MODE_HRV, // Pan band (band 11) only, in 10-bit raw mode
46 : MODE_RAD // Black-body temperature (K) for thermal bands only (4-10),
47 : // 64-bit float
48 : } open_mode_type;
49 :
50 : typedef enum
51 : {
52 : WHOLE_DISK,
53 : RSS, // letterbox of N 1/3 of earth
54 : SPLIT_HRV // the half-width HRV, may be sheared into two block to follow
55 : // the sun W later in the day
56 : } image_shape_type;
57 :
58 : class MSGNRasterBand;
59 :
60 : /************************************************************************/
61 : /* ==================================================================== */
62 : /* MSGNDataset */
63 : /* ==================================================================== */
64 : /************************************************************************/
65 :
66 : class MSGNDataset final : public GDALDataset
67 : {
68 : friend class MSGNRasterBand;
69 :
70 : VSILFILE *fp;
71 :
72 : Msg_reader_core *msg_reader_core;
73 : open_mode_type m_open_mode = MODE_VISIR;
74 : image_shape_type m_Shape = WHOLE_DISK;
75 : int m_nHRVSplitLine = 0;
76 : int m_nHRVLowerShiftX = 0;
77 : int m_nHRVUpperShiftX = 0;
78 : double adfGeoTransform[6];
79 : OGRSpatialReference m_oSRS{};
80 :
81 : public:
82 : MSGNDataset();
83 : ~MSGNDataset();
84 :
85 : static GDALDataset *Open(GDALOpenInfo *);
86 :
87 : CPLErr GetGeoTransform(double *padfTransform) override;
88 : const OGRSpatialReference *GetSpatialRef() const override;
89 : };
90 :
91 : /************************************************************************/
92 : /* ==================================================================== */
93 : /* MSGNRasterBand */
94 : /* ==================================================================== */
95 : /************************************************************************/
96 :
97 : class MSGNRasterBand final : public GDALRasterBand
98 : {
99 : friend class MSGNDataset;
100 :
101 : unsigned int packet_size;
102 : unsigned int bytes_per_line;
103 : unsigned int interline_spacing;
104 : unsigned int orig_band_no; // The name of the band
105 : unsigned int band_in_file; // The effective index of the band in the file
106 : open_mode_type open_mode;
107 :
108 0 : double GetNoDataValue(int *pbSuccess = nullptr) override
109 : {
110 0 : if (pbSuccess)
111 : {
112 0 : *pbSuccess = 1;
113 : }
114 0 : return MSGN_NODATA_VALUE;
115 : }
116 :
117 : double MSGN_NODATA_VALUE;
118 :
119 : char band_description[30];
120 :
121 : public:
122 : MSGNRasterBand(MSGNDataset *, int, open_mode_type mode, int orig_band_no,
123 : int band_in_file);
124 :
125 : virtual CPLErr IReadBlock(int, int, void *) override;
126 : virtual double GetMinimum(int *pbSuccess = nullptr) override;
127 : virtual double GetMaximum(int *pbSuccess = nullptr) override;
128 :
129 0 : virtual const char *GetDescription() const override
130 : {
131 0 : return band_description;
132 : }
133 : };
134 :
135 : /************************************************************************/
136 : /* MSGNRasterBand() */
137 : /************************************************************************/
138 :
139 0 : MSGNRasterBand::MSGNRasterBand(MSGNDataset *poDSIn, int nBandIn,
140 : open_mode_type mode, int orig_band_noIn,
141 0 : int band_in_fileIn)
142 : : packet_size(0), bytes_per_line(0),
143 0 : interline_spacing(poDSIn->msg_reader_core->get_interline_spacing()),
144 0 : orig_band_no(orig_band_noIn), band_in_file(band_in_fileIn),
145 0 : open_mode(mode)
146 : {
147 0 : poDS = poDSIn;
148 0 : nBand = nBandIn; // GDAL's band number, i.e. always starts at 1.
149 :
150 0 : snprintf(band_description, sizeof(band_description), "band %02u",
151 : orig_band_no);
152 :
153 0 : if (mode != MODE_RAD)
154 : {
155 0 : eDataType = GDT_UInt16;
156 0 : MSGN_NODATA_VALUE = 0;
157 : }
158 : else
159 : {
160 0 : eDataType = GDT_Float64;
161 0 : MSGN_NODATA_VALUE = -1000;
162 : }
163 :
164 0 : nBlockXSize = poDS->GetRasterXSize();
165 0 : nBlockYSize = 1;
166 :
167 0 : if (mode != MODE_HRV)
168 : {
169 0 : packet_size = poDSIn->msg_reader_core->get_visir_packet_size();
170 0 : bytes_per_line = poDSIn->msg_reader_core->get_visir_bytes_per_line();
171 : }
172 : else
173 : {
174 0 : packet_size = poDSIn->msg_reader_core->get_hrv_packet_size();
175 0 : bytes_per_line = poDSIn->msg_reader_core->get_hrv_bytes_per_line();
176 : }
177 0 : }
178 :
179 : /************************************************************************/
180 : /* IReadBlock() */
181 : /************************************************************************/
182 :
183 0 : CPLErr MSGNRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
184 : void *pImage)
185 :
186 : {
187 0 : MSGNDataset *poGDS = (MSGNDataset *)poDS;
188 :
189 : // invert y position
190 0 : const int i_nBlockYOff = poDS->GetRasterYSize() - 1 - nBlockYOff;
191 :
192 0 : const int nSamples = static_cast<int>((bytes_per_line * 8) / 10);
193 0 : if (poGDS->m_Shape == WHOLE_DISK && nRasterXSize != nSamples)
194 : {
195 0 : CPLError(CE_Failure, CPLE_AppDefined, "nRasterXSize %d != nSamples %d",
196 : nRasterXSize, nSamples);
197 0 : return CE_Failure;
198 : }
199 :
200 0 : unsigned int data_length =
201 0 : bytes_per_line + (unsigned int)sizeof(SUB_VISIRLINE);
202 0 : vsi_l_offset data_offset = 0;
203 :
204 0 : if (open_mode != MODE_HRV)
205 : {
206 0 : data_offset =
207 0 : poGDS->msg_reader_core->get_f_data_offset() +
208 0 : static_cast<vsi_l_offset>(interline_spacing) * i_nBlockYOff +
209 0 : static_cast<vsi_l_offset>(band_in_file - 1) * packet_size +
210 0 : (packet_size - data_length);
211 : }
212 : else
213 : {
214 0 : data_offset =
215 0 : poGDS->msg_reader_core->get_f_data_offset() +
216 0 : static_cast<vsi_l_offset>(interline_spacing) *
217 0 : (int(i_nBlockYOff / 3) + 1) -
218 0 : static_cast<vsi_l_offset>(packet_size) * (3 - (i_nBlockYOff % 3)) +
219 0 : (packet_size - data_length);
220 : }
221 :
222 0 : if (VSIFSeekL(poGDS->fp, data_offset, SEEK_SET) != 0)
223 0 : return CE_Failure;
224 :
225 0 : char *pszRecord = (char *)CPLMalloc(data_length);
226 0 : size_t nread = VSIFReadL(pszRecord, 1, data_length, poGDS->fp);
227 :
228 0 : SUB_VISIRLINE *p = (SUB_VISIRLINE *)pszRecord;
229 0 : to_native(*p);
230 :
231 0 : if (p->lineValidity != 1 || poGDS->m_Shape != WHOLE_DISK)
232 : { // Split lines are not full width, so NODATA all first
233 0 : for (int c = 0; c < nBlockXSize; c++)
234 : {
235 0 : if (open_mode != MODE_RAD)
236 : {
237 0 : ((GUInt16 *)pImage)[c] = (GUInt16)MSGN_NODATA_VALUE;
238 : }
239 : else
240 : {
241 0 : ((double *)pImage)[c] = MSGN_NODATA_VALUE;
242 : }
243 : }
244 : }
245 :
246 0 : if (nread != data_length ||
247 0 : (p->lineNumberInVisirGrid -
248 0 : ((open_mode == MODE_HRV && poGDS->m_Shape == RSS)
249 0 : ? (3 * poGDS->msg_reader_core->get_line_start()) - 2
250 0 : : poGDS->msg_reader_core->get_line_start())) !=
251 0 : (unsigned int)i_nBlockYOff)
252 : {
253 0 : CPLDebug("MSGN", "Shape %s",
254 0 : poGDS->m_Shape == RSS
255 : ? "RSS"
256 0 : : (poGDS->m_Shape == WHOLE_DISK ? "whole" : "split HRV"));
257 :
258 0 : CPLDebug(
259 : "MSGN", "nread = %lu, data_len %d, linenum %d, start %d, offset %d",
260 : (long unsigned int)nread, // Mingw_w64 otherwise wants %llu - MSG
261 : // read will never exceed 32 bits
262 : data_length, (p->lineNumberInVisirGrid),
263 0 : poGDS->msg_reader_core->get_line_start(), i_nBlockYOff);
264 :
265 0 : CPLFree(pszRecord);
266 :
267 0 : CPLError(CE_Failure, CPLE_AppDefined, "MSGN Scanline corrupt.");
268 :
269 0 : return CE_Failure;
270 : }
271 :
272 : // unpack the 10-bit values into 16-bit unsigned short ints
273 0 : unsigned char *cptr =
274 0 : (unsigned char *)pszRecord + (data_length - bytes_per_line);
275 0 : int bitsLeft = 8;
276 :
277 0 : if (open_mode != MODE_RAD)
278 : {
279 0 : int shift = 0;
280 0 : if (poGDS->m_Shape == SPLIT_HRV)
281 0 : shift = i_nBlockYOff < poGDS->m_nHRVSplitLine
282 0 : ? poGDS->m_nHRVLowerShiftX
283 : : poGDS->m_nHRVUpperShiftX;
284 0 : for (int c = 0; c < nSamples; c++)
285 : {
286 0 : unsigned short value = 0;
287 0 : for (int bit = 0; bit < 10; bit++)
288 : {
289 0 : value <<= 1;
290 0 : if (*cptr & 128)
291 : {
292 0 : value |= 1;
293 : }
294 0 : *cptr <<= 1;
295 0 : bitsLeft--;
296 0 : if (bitsLeft == 0)
297 : {
298 0 : cptr++;
299 0 : bitsLeft = 8;
300 : }
301 : }
302 0 : ((GUInt16 *)pImage)[nBlockXSize - 1 - c - shift] = value;
303 : }
304 : }
305 : else
306 : {
307 : // radiance mode
308 0 : for (int c = 0; c < nSamples; c++)
309 : {
310 0 : unsigned short value = 0;
311 0 : for (int bit = 0; bit < 10; bit++)
312 : {
313 0 : value <<= 1;
314 0 : if (*cptr & 128)
315 : {
316 0 : value |= 1;
317 : }
318 0 : *cptr <<= 1;
319 0 : bitsLeft--;
320 0 : if (bitsLeft == 0)
321 : {
322 0 : cptr++;
323 0 : bitsLeft = 8;
324 : }
325 : }
326 0 : double dvalue = double(value);
327 : double bbvalue =
328 0 : dvalue * poGDS->msg_reader_core
329 0 : ->get_calibration_parameters()[orig_band_no - 1]
330 0 : .cal_slope +
331 0 : poGDS->msg_reader_core
332 0 : ->get_calibration_parameters()[orig_band_no - 1]
333 0 : .cal_offset;
334 :
335 0 : ((double *)pImage)[nBlockXSize - 1 - c] = bbvalue;
336 : }
337 : }
338 0 : CPLFree(pszRecord);
339 0 : return CE_None;
340 : }
341 :
342 : /************************************************************************/
343 : /* GetMinimum() */
344 : /************************************************************************/
345 0 : double MSGNRasterBand::GetMinimum(int *pbSuccess)
346 : {
347 0 : if (pbSuccess)
348 : {
349 0 : *pbSuccess = 1;
350 : }
351 0 : return open_mode != MODE_RAD ? 1 : GDALRasterBand::GetMinimum(pbSuccess);
352 : }
353 :
354 : /************************************************************************/
355 : /* GetMaximum() */
356 : /************************************************************************/
357 0 : double MSGNRasterBand::GetMaximum(int *pbSuccess)
358 : {
359 0 : if (pbSuccess)
360 : {
361 0 : *pbSuccess = 1;
362 : }
363 0 : return open_mode != MODE_RAD ? 1023 : GDALRasterBand::GetMaximum(pbSuccess);
364 : }
365 :
366 : /************************************************************************/
367 : /* ==================================================================== */
368 : /* MSGNDataset */
369 : /* ==================================================================== */
370 : /************************************************************************/
371 :
372 0 : MSGNDataset::MSGNDataset() : fp(nullptr), msg_reader_core(nullptr)
373 : {
374 0 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
375 0 : std::fill_n(adfGeoTransform, CPL_ARRAYSIZE(adfGeoTransform), 0);
376 0 : }
377 :
378 : /************************************************************************/
379 : /* ~MSGNDataset() */
380 : /************************************************************************/
381 :
382 0 : MSGNDataset::~MSGNDataset()
383 :
384 : {
385 0 : if (fp != nullptr)
386 0 : VSIFCloseL(fp);
387 :
388 0 : if (msg_reader_core)
389 : {
390 0 : delete msg_reader_core;
391 : }
392 0 : }
393 :
394 : /************************************************************************/
395 : /* GetGeoTransform() */
396 : /************************************************************************/
397 :
398 0 : CPLErr MSGNDataset::GetGeoTransform(double *padfTransform)
399 :
400 : {
401 0 : for (int i = 0; i < 6; i++)
402 : {
403 0 : padfTransform[i] = adfGeoTransform[i];
404 : }
405 :
406 0 : return CE_None;
407 : }
408 :
409 : /************************************************************************/
410 : /* GetSpatialRef() */
411 : /************************************************************************/
412 :
413 0 : const OGRSpatialReference *MSGNDataset::GetSpatialRef() const
414 :
415 : {
416 0 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
417 : }
418 :
419 : /************************************************************************/
420 : /* Open() */
421 : /************************************************************************/
422 :
423 30371 : GDALDataset *MSGNDataset::Open(GDALOpenInfo *poOpenInfo)
424 :
425 : {
426 30371 : open_mode_type open_mode = MODE_VISIR;
427 30371 : GDALOpenInfo *open_info = poOpenInfo;
428 30363 : std::unique_ptr<GDALOpenInfo> poOpenInfoToFree;
429 :
430 30371 : if (!poOpenInfo->bStatOK)
431 : {
432 25789 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "HRV:"))
433 : {
434 0 : poOpenInfoToFree = std::make_unique<GDALOpenInfo>(
435 0 : &poOpenInfo->pszFilename[4], poOpenInfo->eAccess);
436 0 : open_info = poOpenInfoToFree.get();
437 0 : open_mode = MODE_HRV;
438 : }
439 25789 : else if (STARTS_WITH_CI(poOpenInfo->pszFilename, "RAD:"))
440 : {
441 0 : poOpenInfoToFree = std::make_unique<GDALOpenInfo>(
442 0 : &poOpenInfo->pszFilename[4], poOpenInfo->eAccess);
443 0 : open_info = poOpenInfoToFree.get();
444 0 : open_mode = MODE_RAD;
445 : }
446 : }
447 :
448 : /* -------------------------------------------------------------------- */
449 : /* Before trying MSGNOpen() we first verify that there is at */
450 : /* least one "\n#keyword" type signature in the first chunk of */
451 : /* the file. */
452 : /* -------------------------------------------------------------------- */
453 30363 : if (open_info->fpL == nullptr || open_info->nHeaderBytes < 50)
454 : {
455 26532 : return nullptr;
456 : }
457 :
458 : /* check if this is a "NATIVE" MSG format image */
459 3831 : if (!STARTS_WITH_CI((char *)open_info->pabyHeader,
460 : "FormatName : NATIVE"))
461 : {
462 3831 : return nullptr;
463 : }
464 :
465 : /* -------------------------------------------------------------------- */
466 : /* Confirm the requested access is supported. */
467 : /* -------------------------------------------------------------------- */
468 0 : if (poOpenInfo->eAccess == GA_Update)
469 : {
470 0 : CPLError(CE_Failure, CPLE_NotSupported,
471 : "The MSGN driver does not support update access to existing"
472 : " datasets.\n");
473 0 : return nullptr;
474 : }
475 :
476 : /* -------------------------------------------------------------------- */
477 : /* Create a corresponding GDALDataset. */
478 : /* -------------------------------------------------------------------- */
479 0 : VSILFILE *fp = VSIFOpenL(open_info->pszFilename, "rb");
480 0 : if (fp == nullptr)
481 : {
482 0 : return nullptr;
483 : }
484 :
485 0 : auto poDS = std::make_unique<MSGNDataset>();
486 :
487 0 : poDS->m_open_mode = open_mode;
488 0 : poDS->fp = fp;
489 :
490 : /* -------------------------------------------------------------------- */
491 : /* Read the header. */
492 : /* -------------------------------------------------------------------- */
493 : // first reset the file pointer, then hand over to the msg_reader_core
494 0 : CPL_IGNORE_RET_VAL(VSIFSeekL(poDS->fp, 0, SEEK_SET));
495 :
496 0 : poDS->msg_reader_core = new Msg_reader_core(poDS->fp);
497 :
498 0 : if (!poDS->msg_reader_core->get_open_success())
499 : {
500 0 : return nullptr;
501 : }
502 :
503 0 : poDS->nRasterXSize = poDS->msg_reader_core->get_columns();
504 0 : poDS->nRasterYSize = poDS->msg_reader_core->get_lines();
505 :
506 0 : if (open_mode == MODE_HRV)
507 : {
508 : const int nRawHRVColumns =
509 0 : (poDS->msg_reader_core->get_hrv_bytes_per_line() * 8) / 10;
510 0 : poDS->nRasterYSize *= 3;
511 0 : const auto &idr = poDS->msg_reader_core->get_image_description_record();
512 : // Check if the split layout of the HRV channel meets our expectations
513 : // to re-assemble it in a consistent way
514 0 : CPLDebug("MSGN", "HRV raw col %d raster X %d raster Y %d",
515 0 : nRawHRVColumns, poDS->nRasterXSize, poDS->nRasterYSize);
516 :
517 0 : if (idr.plannedCoverage_hrv.lowerSouthLinePlanned == 1 &&
518 0 : idr.plannedCoverage_hrv.lowerNorthLinePlanned > 1 &&
519 0 : idr.plannedCoverage_hrv.lowerNorthLinePlanned <
520 0 : poDS->nRasterYSize &&
521 0 : idr.plannedCoverage_hrv.upperSouthLinePlanned ==
522 0 : idr.plannedCoverage_hrv.lowerNorthLinePlanned + 1 &&
523 0 : idr.plannedCoverage_hrv.upperNorthLinePlanned ==
524 0 : poDS->nRasterYSize &&
525 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned >= 1 &&
526 0 : idr.plannedCoverage_hrv.lowerWestColumnPlanned ==
527 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned +
528 0 : nRawHRVColumns - 1 &&
529 0 : idr.plannedCoverage_hrv.lowerWestColumnPlanned <=
530 0 : poDS->nRasterXSize * 3 &&
531 0 : idr.plannedCoverage_hrv.upperEastColumnPlanned >= 1 &&
532 0 : idr.plannedCoverage_hrv.upperWestColumnPlanned ==
533 0 : idr.plannedCoverage_hrv.upperEastColumnPlanned +
534 0 : nRawHRVColumns - 1 &&
535 0 : idr.plannedCoverage_hrv.upperWestColumnPlanned <=
536 0 : poDS->nRasterXSize * 3)
537 : {
538 0 : poDS->nRasterXSize *= 3;
539 0 : poDS->m_Shape = SPLIT_HRV;
540 0 : poDS->m_nHRVSplitLine =
541 0 : idr.plannedCoverage_hrv.upperSouthLinePlanned;
542 0 : poDS->m_nHRVLowerShiftX =
543 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned - 1;
544 0 : poDS->m_nHRVUpperShiftX =
545 0 : idr.plannedCoverage_hrv.upperEastColumnPlanned - 1;
546 : }
547 0 : else if (idr.plannedCoverage_hrv.upperNorthLinePlanned == 0 &&
548 0 : idr.plannedCoverage_hrv.upperSouthLinePlanned == 0 &&
549 0 : idr.plannedCoverage_hrv.upperWestColumnPlanned == 0 &&
550 0 : idr.plannedCoverage_hrv.upperEastColumnPlanned ==
551 0 : 0 && // RSS only uses the lower section
552 0 : idr.plannedCoverage_hrv.lowerNorthLinePlanned ==
553 0 : idr.referencegrid_hrv.numberOfLines && // start at max N
554 : // full expected width
555 0 : idr.plannedCoverage_hrv.lowerWestColumnPlanned ==
556 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned +
557 0 : nRawHRVColumns - 1 &&
558 0 : idr.plannedCoverage_hrv.lowerSouthLinePlanned > 1 &&
559 0 : idr.plannedCoverage_hrv.lowerSouthLinePlanned <
560 0 : idr.referencegrid_hrv.numberOfLines &&
561 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned >= 1 &&
562 0 : idr.plannedCoverage_hrv.lowerWestColumnPlanned <=
563 0 : poDS->nRasterXSize * 3 &&
564 : // full height
565 0 : idr.plannedCoverage_hrv.lowerNorthLinePlanned ==
566 0 : idr.plannedCoverage_hrv.lowerSouthLinePlanned +
567 0 : poDS->nRasterYSize - 1)
568 : {
569 0 : poDS->nRasterXSize *= 3;
570 0 : poDS->m_Shape = RSS;
571 : }
572 : else
573 : {
574 0 : CPLError(
575 : CE_Failure, CPLE_AppDefined,
576 : "HRV neither Whole Disk nor RSS - don't know how to handle");
577 0 : return nullptr;
578 : }
579 : }
580 : else
581 : {
582 : const int nRawVisIRColumns =
583 0 : (poDS->msg_reader_core->get_visir_bytes_per_line() * 8) / 10;
584 :
585 0 : const auto &idr = poDS->msg_reader_core->get_image_description_record();
586 : // Check if the VisIR channel is RSS or not, and if it meets our
587 : // expectations to re-assemble it in a consistent way
588 0 : CPLDebug("MSGN", "raw col %d raster X %d raster Y %d", nRawVisIRColumns,
589 0 : poDS->nRasterXSize, poDS->nRasterYSize);
590 :
591 0 : if (idr.plannedCoverage_visir.southernLinePlanned == 1 &&
592 0 : idr.plannedCoverage_visir.northernLinePlanned ==
593 0 : poDS->nRasterYSize &&
594 0 : idr.plannedCoverage_visir.easternColumnPlanned >= 1 &&
595 0 : idr.plannedCoverage_visir.westernColumnPlanned ==
596 0 : idr.plannedCoverage_visir.easternColumnPlanned +
597 0 : nRawVisIRColumns - 1 &&
598 0 : idr.plannedCoverage_visir.westernColumnPlanned <=
599 0 : poDS->nRasterXSize)
600 : {
601 0 : poDS->m_Shape = WHOLE_DISK;
602 : }
603 0 : else if (idr.plannedCoverage_visir.northernLinePlanned ==
604 0 : idr.referencegrid_visir.numberOfLines && // start at max N
605 : // full expected width
606 0 : idr.plannedCoverage_visir.westernColumnPlanned ==
607 0 : idr.plannedCoverage_visir.easternColumnPlanned +
608 0 : nRawVisIRColumns - 1 &&
609 0 : idr.plannedCoverage_visir.southernLinePlanned > 1 &&
610 0 : idr.plannedCoverage_visir.easternColumnPlanned >= 1 &&
611 0 : idr.plannedCoverage_visir.westernColumnPlanned <=
612 0 : poDS->nRasterXSize &&
613 : // full height
614 0 : idr.plannedCoverage_visir.northernLinePlanned ==
615 0 : idr.plannedCoverage_visir.southernLinePlanned +
616 0 : poDS->nRasterYSize - 1)
617 : {
618 0 : poDS->m_Shape = RSS;
619 : }
620 : else
621 : {
622 0 : CPLError(CE_Failure, CPLE_AppDefined,
623 : "Neither Whole Disk nor RSS - don't know how to handle");
624 0 : return nullptr;
625 : }
626 : }
627 :
628 0 : CPLDebug("MSGN", "Shape %s",
629 0 : poDS->m_Shape == RSS
630 : ? "RSS"
631 0 : : (poDS->m_Shape == WHOLE_DISK ? "whole" : "split HRV"));
632 :
633 : /* -------------------------------------------------------------------- */
634 : /* Create band information objects. */
635 : /* -------------------------------------------------------------------- */
636 : unsigned int i;
637 0 : unsigned int band_count = 1;
638 0 : unsigned int missing_band_count = 0;
639 0 : const unsigned char *bands = poDS->msg_reader_core->get_band_map();
640 0 : unsigned char band_map[MSG_NUM_CHANNELS + 1] = {
641 : 0}; // map GDAL band numbers to MSG channels
642 0 : for (i = 0; i < MSG_NUM_CHANNELS; i++)
643 : {
644 0 : if (bands[i])
645 : {
646 0 : bool ok_to_add = false;
647 0 : switch (open_mode)
648 : {
649 0 : case MODE_VISIR:
650 0 : ok_to_add = i < MSG_NUM_CHANNELS - 1;
651 0 : break;
652 0 : case MODE_RAD:
653 0 : ok_to_add = (i <= 2) ||
654 0 : (Msg_reader_core::Blackbody_LUT[i + 1].B != 0);
655 0 : break;
656 0 : case MODE_HRV:
657 0 : ok_to_add = i == MSG_NUM_CHANNELS - 1;
658 0 : break;
659 : }
660 0 : if (ok_to_add)
661 : {
662 0 : poDS->SetBand(band_count,
663 0 : new MSGNRasterBand(poDS.get(), band_count,
664 0 : open_mode, i + 1,
665 0 : i + 1 - missing_band_count));
666 0 : band_map[band_count] = (unsigned char)(i + 1);
667 0 : band_count++;
668 : }
669 : }
670 : else
671 : {
672 0 : missing_band_count++;
673 : }
674 : }
675 :
676 : double pixel_gsd_x;
677 : double pixel_gsd_y;
678 : double origin_x;
679 : double origin_y;
680 :
681 : {
682 0 : const auto &idr = poDS->msg_reader_core->get_image_description_record();
683 : /* there are a number of 'magic' constants below
684 : I trimmed them to get registration for MSG4, MSG3, MSG2 with country
685 : outlines from
686 : http://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units
687 :
688 : Adjust in two phases P1, P2.
689 : I describe direction as outline being NSEW of coast shape when number
690 : is changed
691 : */
692 :
693 0 : if (open_mode != MODE_HRV)
694 : {
695 0 : pixel_gsd_x =
696 0 : 1000.0 * poDS->msg_reader_core
697 0 : ->get_col_dir_step(); // convert from km to m
698 0 : pixel_gsd_y =
699 0 : 1000.0 * poDS->msg_reader_core
700 0 : ->get_line_dir_step(); // convert from km to m
701 0 : origin_x = -pixel_gsd_x * (-(Conversions::nlines / 2.0) +
702 0 : poDS->msg_reader_core->get_col_start() -
703 : 1); // all vis/NIR E-W -ve E
704 0 : origin_y = -pixel_gsd_y * ((Conversions::nlines / 2.0) -
705 0 : poDS->msg_reader_core->get_line_start() +
706 : 1.5); // set with 4 N-S +ve S
707 : }
708 : else
709 : {
710 0 : pixel_gsd_x =
711 0 : 1000.0 * poDS->msg_reader_core
712 0 : ->get_hrv_col_dir_step(); // convert from km to m
713 0 : pixel_gsd_y =
714 0 : 1000.0 * poDS->msg_reader_core
715 0 : ->get_hrv_line_dir_step(); // convert from km to m
716 0 : if (poDS->m_Shape == RSS)
717 : {
718 0 : origin_x = -pixel_gsd_x *
719 0 : (-(3 * Conversions::nlines / 2.0) -
720 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned -
721 : 1); // MSG3 HRV E-W -ve E
722 0 : origin_y = -pixel_gsd_y *
723 0 : ((3 * Conversions::nlines / 2.0) -
724 0 : idr.plannedCoverage_hrv.lowerSouthLinePlanned +
725 : 2); // N-S -ve S
726 : }
727 : else
728 : {
729 0 : origin_x =
730 0 : -pixel_gsd_x * (-(3 * Conversions::nlines / 2.0) +
731 0 : 1 * poDS->msg_reader_core->get_col_start() -
732 : 3); // MSG4, MSG2 HRV E-W -ve E
733 0 : origin_y = -pixel_gsd_y *
734 0 : ((3 * Conversions::nlines / 2.0) -
735 0 : 1 * poDS->msg_reader_core->get_line_start() +
736 : 4); // N-S +ve S
737 : }
738 : }
739 :
740 : /* the conversion to lat/long is in two parts:
741 : pixels to m (around imaginary circle r=sta height) in the geo
742 : projection (affine transformation) geo to lat/long via the GEOS
743 : projection (in WKT) and the ellipsoid
744 :
745 : CGMS/DOC/12/0017 section 4.4.2
746 : */
747 :
748 0 : poDS->adfGeoTransform[0] = -origin_x;
749 0 : poDS->adfGeoTransform[1] = pixel_gsd_x;
750 0 : poDS->adfGeoTransform[2] = 0.0;
751 :
752 0 : poDS->adfGeoTransform[3] = -origin_y;
753 0 : poDS->adfGeoTransform[4] = 0.0;
754 0 : poDS->adfGeoTransform[5] = -pixel_gsd_y;
755 :
756 0 : poDS->m_oSRS.SetProjCS("Geostationary projection (MSG)");
757 :
758 0 : poDS->m_oSRS.SetGeogCS(
759 : "MSG Ellipsoid", "MSG_DATUM", "MSG_SPHEROID",
760 0 : Conversions::req *
761 : 1000.0, // SetGeogCS doesn't specify length units, so all must
762 : // be the same - here m.
763 0 : 1.0 / Conversions::oblate // 1 / ( 1 -
764 : // Conversions::rpol/Conversions::req)
765 : );
766 :
767 0 : poDS->m_oSRS.SetGEOS(
768 0 : idr.longitudeOfSSP,
769 0 : (Conversions::altitude - Conversions::req) *
770 : 1000.0, // we're using meters as length unit
771 : 0.0,
772 : // false northing to handle the fact RSS is only 1/3 disk
773 : pixel_gsd_y *
774 0 : ((poDS->m_Shape == RSS)
775 0 : ? ((open_mode != MODE_HRV)
776 0 : ? -(idr.plannedCoverage_visir.southernLinePlanned -
777 : 1)
778 : : // MSG-3 vis/NIR N-S P2
779 0 : -(idr.plannedCoverage_hrv.lowerSouthLinePlanned +
780 : 1))
781 : : // MSG-3 HRV N-S P2 -ve N
782 : 0.0));
783 : }
784 :
785 : const CALIBRATION *cal =
786 0 : poDS->msg_reader_core->get_calibration_parameters();
787 : char tagname[30];
788 : char field[300];
789 :
790 0 : poDS->SetMetadataItem("Radiometric parameters format", "offset slope");
791 0 : for (i = 1; i < band_count; i++)
792 : {
793 0 : snprintf(tagname, sizeof(tagname), "ch%02u_cal", band_map[i]);
794 0 : CPLsnprintf(field, sizeof(field), "%.12e %.12e",
795 0 : cal[band_map[i] - 1].cal_offset,
796 0 : cal[band_map[i] - 1].cal_slope);
797 0 : poDS->SetMetadataItem(tagname, field);
798 : }
799 :
800 0 : snprintf(
801 : field, sizeof(field), "%04u%02u%02u/%02u:%02u",
802 0 : poDS->msg_reader_core->get_year(), poDS->msg_reader_core->get_month(),
803 0 : poDS->msg_reader_core->get_day(), poDS->msg_reader_core->get_hour(),
804 0 : poDS->msg_reader_core->get_minute());
805 0 : poDS->SetMetadataItem("Date/Time", field);
806 :
807 0 : snprintf(field, sizeof(field), "%u %u",
808 0 : poDS->msg_reader_core->get_line_start(),
809 0 : poDS->msg_reader_core->get_col_start());
810 0 : poDS->SetMetadataItem("Origin", field);
811 :
812 0 : return poDS.release();
813 : }
814 :
815 : /************************************************************************/
816 : /* GDALRegister_MSGN() */
817 : /************************************************************************/
818 :
819 1509 : void GDALRegister_MSGN()
820 :
821 : {
822 1509 : if (GDALGetDriverByName("MSGN") != nullptr)
823 295 : return;
824 :
825 1214 : GDALDriver *poDriver = new GDALDriver();
826 :
827 1214 : poDriver->SetDescription("MSGN");
828 1214 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
829 1214 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
830 1214 : "EUMETSAT Archive native (.nat)");
831 1214 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/msgn.html");
832 1214 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "nat");
833 1214 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
834 :
835 1214 : poDriver->pfnOpen = MSGNDataset::Open;
836 :
837 1214 : GetGDALDriverManager()->RegisterDriver(poDriver);
838 : }
|
