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 : double adfGeoTransform[6];
63 : OGRSpatialReference m_oSRS{};
64 :
65 : public:
66 : MSGNDataset();
67 : ~MSGNDataset();
68 :
69 : static GDALDataset *Open(GDALOpenInfo *);
70 :
71 : CPLErr GetGeoTransform(double *padfTransform) 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 = (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 = (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 : std::fill_n(adfGeoTransform, CPL_ARRAYSIZE(adfGeoTransform), 0);
360 0 : }
361 :
362 : /************************************************************************/
363 : /* ~MSGNDataset() */
364 : /************************************************************************/
365 :
366 0 : MSGNDataset::~MSGNDataset()
367 :
368 : {
369 0 : if (fp != nullptr)
370 0 : VSIFCloseL(fp);
371 :
372 0 : if (msg_reader_core)
373 : {
374 0 : delete msg_reader_core;
375 : }
376 0 : }
377 :
378 : /************************************************************************/
379 : /* GetGeoTransform() */
380 : /************************************************************************/
381 :
382 0 : CPLErr MSGNDataset::GetGeoTransform(double *padfTransform)
383 :
384 : {
385 0 : for (int i = 0; i < 6; i++)
386 : {
387 0 : padfTransform[i] = adfGeoTransform[i];
388 : }
389 :
390 0 : return CE_None;
391 : }
392 :
393 : /************************************************************************/
394 : /* GetSpatialRef() */
395 : /************************************************************************/
396 :
397 0 : const OGRSpatialReference *MSGNDataset::GetSpatialRef() const
398 :
399 : {
400 0 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
401 : }
402 :
403 : /************************************************************************/
404 : /* Open() */
405 : /************************************************************************/
406 :
407 32431 : GDALDataset *MSGNDataset::Open(GDALOpenInfo *poOpenInfo)
408 :
409 : {
410 32431 : open_mode_type open_mode = MODE_VISIR;
411 32431 : GDALOpenInfo *open_info = poOpenInfo;
412 32429 : std::unique_ptr<GDALOpenInfo> poOpenInfoToFree;
413 :
414 32431 : if (!poOpenInfo->bStatOK)
415 : {
416 27621 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "HRV:"))
417 : {
418 0 : poOpenInfoToFree = std::make_unique<GDALOpenInfo>(
419 0 : &poOpenInfo->pszFilename[4], poOpenInfo->eAccess);
420 0 : open_info = poOpenInfoToFree.get();
421 0 : open_mode = MODE_HRV;
422 : }
423 27621 : else if (STARTS_WITH_CI(poOpenInfo->pszFilename, "RAD:"))
424 : {
425 0 : poOpenInfoToFree = std::make_unique<GDALOpenInfo>(
426 0 : &poOpenInfo->pszFilename[4], poOpenInfo->eAccess);
427 0 : open_info = poOpenInfoToFree.get();
428 0 : open_mode = MODE_RAD;
429 : }
430 : }
431 :
432 : /* -------------------------------------------------------------------- */
433 : /* Before trying MSGNOpen() we first verify that there is at */
434 : /* least one "\n#keyword" type signature in the first chunk of */
435 : /* the file. */
436 : /* -------------------------------------------------------------------- */
437 32429 : if (open_info->fpL == nullptr || open_info->nHeaderBytes < 50)
438 : {
439 28469 : return nullptr;
440 : }
441 :
442 : /* check if this is a "NATIVE" MSG format image */
443 3960 : if (!STARTS_WITH_CI((char *)open_info->pabyHeader,
444 : "FormatName : NATIVE"))
445 : {
446 3960 : return nullptr;
447 : }
448 :
449 : /* -------------------------------------------------------------------- */
450 : /* Confirm the requested access is supported. */
451 : /* -------------------------------------------------------------------- */
452 0 : if (poOpenInfo->eAccess == GA_Update)
453 : {
454 0 : ReportUpdateNotSupportedByDriver("MSGN");
455 0 : return nullptr;
456 : }
457 :
458 : /* -------------------------------------------------------------------- */
459 : /* Create a corresponding GDALDataset. */
460 : /* -------------------------------------------------------------------- */
461 0 : VSILFILE *fp = VSIFOpenL(open_info->pszFilename, "rb");
462 0 : if (fp == nullptr)
463 : {
464 0 : return nullptr;
465 : }
466 :
467 0 : auto poDS = std::make_unique<MSGNDataset>();
468 :
469 0 : poDS->m_open_mode = open_mode;
470 0 : poDS->fp = fp;
471 :
472 : /* -------------------------------------------------------------------- */
473 : /* Read the header. */
474 : /* -------------------------------------------------------------------- */
475 : // first reset the file pointer, then hand over to the msg_reader_core
476 0 : CPL_IGNORE_RET_VAL(VSIFSeekL(poDS->fp, 0, SEEK_SET));
477 :
478 0 : poDS->msg_reader_core = new Msg_reader_core(poDS->fp);
479 :
480 0 : if (!poDS->msg_reader_core->get_open_success())
481 : {
482 0 : return nullptr;
483 : }
484 :
485 0 : poDS->nRasterXSize = poDS->msg_reader_core->get_columns();
486 0 : poDS->nRasterYSize = poDS->msg_reader_core->get_lines();
487 :
488 0 : if (open_mode == MODE_HRV)
489 : {
490 : const int nRawHRVColumns =
491 0 : (poDS->msg_reader_core->get_hrv_bytes_per_line() * 8) / 10;
492 0 : poDS->nRasterYSize *= 3;
493 0 : const auto &idr = poDS->msg_reader_core->get_image_description_record();
494 : // Check if the split layout of the HRV channel meets our expectations
495 : // to re-assemble it in a consistent way
496 0 : CPLDebug("MSGN", "HRV raw col %d raster X %d raster Y %d",
497 0 : nRawHRVColumns, poDS->nRasterXSize, poDS->nRasterYSize);
498 :
499 0 : if (idr.plannedCoverage_hrv.lowerSouthLinePlanned == 1 &&
500 0 : idr.plannedCoverage_hrv.lowerNorthLinePlanned > 1 &&
501 0 : idr.plannedCoverage_hrv.lowerNorthLinePlanned <
502 0 : poDS->nRasterYSize &&
503 0 : idr.plannedCoverage_hrv.upperSouthLinePlanned ==
504 0 : idr.plannedCoverage_hrv.lowerNorthLinePlanned + 1 &&
505 0 : idr.plannedCoverage_hrv.upperNorthLinePlanned ==
506 0 : poDS->nRasterYSize &&
507 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned >= 1 &&
508 0 : idr.plannedCoverage_hrv.lowerWestColumnPlanned ==
509 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned +
510 0 : nRawHRVColumns - 1 &&
511 0 : idr.plannedCoverage_hrv.lowerWestColumnPlanned <=
512 0 : poDS->nRasterXSize * 3 &&
513 0 : idr.plannedCoverage_hrv.upperEastColumnPlanned >= 1 &&
514 0 : idr.plannedCoverage_hrv.upperWestColumnPlanned ==
515 0 : idr.plannedCoverage_hrv.upperEastColumnPlanned +
516 0 : nRawHRVColumns - 1 &&
517 0 : idr.plannedCoverage_hrv.upperWestColumnPlanned <=
518 0 : poDS->nRasterXSize * 3)
519 : {
520 0 : poDS->nRasterXSize *= 3;
521 0 : poDS->m_Shape = SPLIT_HRV;
522 0 : poDS->m_nHRVSplitLine =
523 0 : idr.plannedCoverage_hrv.upperSouthLinePlanned;
524 0 : poDS->m_nHRVLowerShiftX =
525 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned - 1;
526 0 : poDS->m_nHRVUpperShiftX =
527 0 : idr.plannedCoverage_hrv.upperEastColumnPlanned - 1;
528 : }
529 0 : else if (idr.plannedCoverage_hrv.upperNorthLinePlanned == 0 &&
530 0 : idr.plannedCoverage_hrv.upperSouthLinePlanned == 0 &&
531 0 : idr.plannedCoverage_hrv.upperWestColumnPlanned == 0 &&
532 0 : idr.plannedCoverage_hrv.upperEastColumnPlanned ==
533 0 : 0 && // RSS only uses the lower section
534 0 : idr.plannedCoverage_hrv.lowerNorthLinePlanned ==
535 0 : idr.referencegrid_hrv.numberOfLines && // start at max N
536 : // full expected width
537 0 : idr.plannedCoverage_hrv.lowerWestColumnPlanned ==
538 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned +
539 0 : nRawHRVColumns - 1 &&
540 0 : idr.plannedCoverage_hrv.lowerSouthLinePlanned > 1 &&
541 0 : idr.plannedCoverage_hrv.lowerSouthLinePlanned <
542 0 : idr.referencegrid_hrv.numberOfLines &&
543 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned >= 1 &&
544 0 : idr.plannedCoverage_hrv.lowerWestColumnPlanned <=
545 0 : poDS->nRasterXSize * 3 &&
546 : // full height
547 0 : idr.plannedCoverage_hrv.lowerNorthLinePlanned ==
548 0 : idr.plannedCoverage_hrv.lowerSouthLinePlanned +
549 0 : poDS->nRasterYSize - 1)
550 : {
551 0 : poDS->nRasterXSize *= 3;
552 0 : poDS->m_Shape = RSS;
553 : }
554 : else
555 : {
556 0 : CPLError(
557 : CE_Failure, CPLE_AppDefined,
558 : "HRV neither Whole Disk nor RSS - don't know how to handle");
559 0 : return nullptr;
560 : }
561 : }
562 : else
563 : {
564 : const int nRawVisIRColumns =
565 0 : (poDS->msg_reader_core->get_visir_bytes_per_line() * 8) / 10;
566 :
567 0 : const auto &idr = poDS->msg_reader_core->get_image_description_record();
568 : // Check if the VisIR channel is RSS or not, and if it meets our
569 : // expectations to re-assemble it in a consistent way
570 0 : CPLDebug("MSGN", "raw col %d raster X %d raster Y %d", nRawVisIRColumns,
571 0 : poDS->nRasterXSize, poDS->nRasterYSize);
572 :
573 0 : if (idr.plannedCoverage_visir.southernLinePlanned == 1 &&
574 0 : idr.plannedCoverage_visir.northernLinePlanned ==
575 0 : poDS->nRasterYSize &&
576 0 : idr.plannedCoverage_visir.easternColumnPlanned >= 1 &&
577 0 : idr.plannedCoverage_visir.westernColumnPlanned ==
578 0 : idr.plannedCoverage_visir.easternColumnPlanned +
579 0 : nRawVisIRColumns - 1 &&
580 0 : idr.plannedCoverage_visir.westernColumnPlanned <=
581 0 : poDS->nRasterXSize)
582 : {
583 0 : poDS->m_Shape = WHOLE_DISK;
584 : }
585 0 : else if (idr.plannedCoverage_visir.northernLinePlanned ==
586 0 : idr.referencegrid_visir.numberOfLines && // start at max N
587 : // full expected width
588 0 : idr.plannedCoverage_visir.westernColumnPlanned ==
589 0 : idr.plannedCoverage_visir.easternColumnPlanned +
590 0 : nRawVisIRColumns - 1 &&
591 0 : idr.plannedCoverage_visir.southernLinePlanned > 1 &&
592 0 : idr.plannedCoverage_visir.easternColumnPlanned >= 1 &&
593 0 : idr.plannedCoverage_visir.westernColumnPlanned <=
594 0 : poDS->nRasterXSize &&
595 : // full height
596 0 : idr.plannedCoverage_visir.northernLinePlanned ==
597 0 : idr.plannedCoverage_visir.southernLinePlanned +
598 0 : poDS->nRasterYSize - 1)
599 : {
600 0 : poDS->m_Shape = RSS;
601 : }
602 : else
603 : {
604 0 : CPLError(CE_Failure, CPLE_AppDefined,
605 : "Neither Whole Disk nor RSS - don't know how to handle");
606 0 : return nullptr;
607 : }
608 : }
609 :
610 0 : CPLDebug("MSGN", "Shape %s",
611 0 : poDS->m_Shape == RSS
612 : ? "RSS"
613 0 : : (poDS->m_Shape == WHOLE_DISK ? "whole" : "split HRV"));
614 :
615 : /* -------------------------------------------------------------------- */
616 : /* Create band information objects. */
617 : /* -------------------------------------------------------------------- */
618 : unsigned int i;
619 0 : unsigned int band_count = 1;
620 0 : unsigned int missing_band_count = 0;
621 0 : const unsigned char *bands = poDS->msg_reader_core->get_band_map();
622 0 : unsigned char band_map[MSG_NUM_CHANNELS + 1] = {
623 : 0}; // map GDAL band numbers to MSG channels
624 0 : for (i = 0; i < MSG_NUM_CHANNELS; i++)
625 : {
626 0 : if (bands[i])
627 : {
628 0 : bool ok_to_add = false;
629 0 : switch (open_mode)
630 : {
631 0 : case MODE_VISIR:
632 0 : ok_to_add = i < MSG_NUM_CHANNELS - 1;
633 0 : break;
634 0 : case MODE_RAD:
635 0 : ok_to_add = (i <= 2) ||
636 0 : (Msg_reader_core::Blackbody_LUT[i + 1].B != 0);
637 0 : break;
638 0 : case MODE_HRV:
639 0 : ok_to_add = i == MSG_NUM_CHANNELS - 1;
640 0 : break;
641 : }
642 0 : if (ok_to_add)
643 : {
644 0 : poDS->SetBand(band_count,
645 0 : new MSGNRasterBand(poDS.get(), band_count,
646 0 : open_mode, i + 1,
647 0 : i + 1 - missing_band_count));
648 0 : band_map[band_count] = (unsigned char)(i + 1);
649 0 : band_count++;
650 : }
651 : }
652 : else
653 : {
654 0 : missing_band_count++;
655 : }
656 : }
657 :
658 : double pixel_gsd_x;
659 : double pixel_gsd_y;
660 : double origin_x;
661 : double origin_y;
662 :
663 : {
664 0 : const auto &idr = poDS->msg_reader_core->get_image_description_record();
665 : /* there are a number of 'magic' constants below
666 : I trimmed them to get registration for MSG4, MSG3, MSG2 with country
667 : outlines from
668 : http://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units
669 :
670 : Adjust in two phases P1, P2.
671 : I describe direction as outline being NSEW of coast shape when number
672 : is changed
673 : */
674 :
675 0 : if (open_mode != MODE_HRV)
676 : {
677 0 : pixel_gsd_x =
678 0 : 1000.0 * poDS->msg_reader_core
679 0 : ->get_col_dir_step(); // convert from km to m
680 0 : pixel_gsd_y =
681 0 : 1000.0 * poDS->msg_reader_core
682 0 : ->get_line_dir_step(); // convert from km to m
683 0 : origin_x = -pixel_gsd_x * (-(Conversions::nlines / 2.0) +
684 0 : poDS->msg_reader_core->get_col_start() -
685 : 1); // all vis/NIR E-W -ve E
686 0 : origin_y = -pixel_gsd_y * ((Conversions::nlines / 2.0) -
687 0 : poDS->msg_reader_core->get_line_start() +
688 : 1.5); // set with 4 N-S +ve S
689 : }
690 : else
691 : {
692 0 : pixel_gsd_x =
693 0 : 1000.0 * poDS->msg_reader_core
694 0 : ->get_hrv_col_dir_step(); // convert from km to m
695 0 : pixel_gsd_y =
696 0 : 1000.0 * poDS->msg_reader_core
697 0 : ->get_hrv_line_dir_step(); // convert from km to m
698 0 : if (poDS->m_Shape == RSS)
699 : {
700 0 : origin_x = -pixel_gsd_x *
701 0 : (-(3 * Conversions::nlines / 2.0) -
702 0 : idr.plannedCoverage_hrv.lowerEastColumnPlanned -
703 : 1); // MSG3 HRV E-W -ve E
704 0 : origin_y = -pixel_gsd_y *
705 0 : ((3 * Conversions::nlines / 2.0) -
706 0 : idr.plannedCoverage_hrv.lowerSouthLinePlanned +
707 : 2); // N-S -ve S
708 : }
709 : else
710 : {
711 0 : origin_x =
712 0 : -pixel_gsd_x * (-(3 * Conversions::nlines / 2.0) +
713 0 : 1 * poDS->msg_reader_core->get_col_start() -
714 : 3); // MSG4, MSG2 HRV E-W -ve E
715 0 : origin_y = -pixel_gsd_y *
716 0 : ((3 * Conversions::nlines / 2.0) -
717 0 : 1 * poDS->msg_reader_core->get_line_start() +
718 : 4); // N-S +ve S
719 : }
720 : }
721 :
722 : /* the conversion to lat/long is in two parts:
723 : pixels to m (around imaginary circle r=sta height) in the geo
724 : projection (affine transformation) geo to lat/long via the GEOS
725 : projection (in WKT) and the ellipsoid
726 :
727 : CGMS/DOC/12/0017 section 4.4.2
728 : */
729 :
730 0 : poDS->adfGeoTransform[0] = -origin_x;
731 0 : poDS->adfGeoTransform[1] = pixel_gsd_x;
732 0 : poDS->adfGeoTransform[2] = 0.0;
733 :
734 0 : poDS->adfGeoTransform[3] = -origin_y;
735 0 : poDS->adfGeoTransform[4] = 0.0;
736 0 : poDS->adfGeoTransform[5] = -pixel_gsd_y;
737 :
738 0 : poDS->m_oSRS.SetProjCS("Geostationary projection (MSG)");
739 :
740 0 : poDS->m_oSRS.SetGeogCS(
741 : "MSG Ellipsoid", "MSG_DATUM", "MSG_SPHEROID",
742 0 : Conversions::req *
743 : 1000.0, // SetGeogCS doesn't specify length units, so all must
744 : // be the same - here m.
745 0 : 1.0 / Conversions::oblate // 1 / ( 1 -
746 : // Conversions::rpol/Conversions::req)
747 : );
748 :
749 0 : poDS->m_oSRS.SetGEOS(
750 0 : idr.longitudeOfSSP,
751 0 : (Conversions::altitude - Conversions::req) *
752 : 1000.0, // we're using meters as length unit
753 : 0.0,
754 : // false northing to handle the fact RSS is only 1/3 disk
755 : pixel_gsd_y *
756 0 : ((poDS->m_Shape == RSS)
757 0 : ? ((open_mode != MODE_HRV)
758 0 : ? -(idr.plannedCoverage_visir.southernLinePlanned -
759 : 1)
760 : : // MSG-3 vis/NIR N-S P2
761 0 : -(idr.plannedCoverage_hrv.lowerSouthLinePlanned +
762 : 1))
763 : : // MSG-3 HRV N-S P2 -ve N
764 : 0.0));
765 : }
766 :
767 : const CALIBRATION *cal =
768 0 : poDS->msg_reader_core->get_calibration_parameters();
769 : char tagname[30];
770 : char field[300];
771 :
772 0 : poDS->SetMetadataItem("Radiometric parameters format", "offset slope");
773 0 : for (i = 1; i < band_count; i++)
774 : {
775 0 : snprintf(tagname, sizeof(tagname), "ch%02u_cal", band_map[i]);
776 0 : CPLsnprintf(field, sizeof(field), "%.12e %.12e",
777 0 : cal[band_map[i] - 1].cal_offset,
778 0 : cal[band_map[i] - 1].cal_slope);
779 0 : poDS->SetMetadataItem(tagname, field);
780 : }
781 :
782 0 : snprintf(
783 : field, sizeof(field), "%04u%02u%02u/%02u:%02u",
784 0 : poDS->msg_reader_core->get_year(), poDS->msg_reader_core->get_month(),
785 0 : poDS->msg_reader_core->get_day(), poDS->msg_reader_core->get_hour(),
786 0 : poDS->msg_reader_core->get_minute());
787 0 : poDS->SetMetadataItem("Date/Time", field);
788 :
789 0 : snprintf(field, sizeof(field), "%u %u",
790 0 : poDS->msg_reader_core->get_line_start(),
791 0 : poDS->msg_reader_core->get_col_start());
792 0 : poDS->SetMetadataItem("Origin", field);
793 :
794 0 : return poDS.release();
795 : }
796 :
797 : /************************************************************************/
798 : /* GDALRegister_MSGN() */
799 : /************************************************************************/
800 :
801 1686 : void GDALRegister_MSGN()
802 :
803 : {
804 1686 : if (GDALGetDriverByName("MSGN") != nullptr)
805 302 : return;
806 :
807 1384 : GDALDriver *poDriver = new GDALDriver();
808 :
809 1384 : poDriver->SetDescription("MSGN");
810 1384 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
811 1384 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
812 1384 : "EUMETSAT Archive native (.nat)");
813 1384 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/msgn.html");
814 1384 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "nat");
815 1384 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
816 :
817 1384 : poDriver->pfnOpen = MSGNDataset::Open;
818 :
819 1384 : GetGDALDriverManager()->RegisterDriver(poDriver);
820 : }
|
