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