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