Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: NOAA Polar Orbiter Level 1b Dataset Reader (AVHRR)
4 : * Purpose: Can read NOAA-9(F)-NOAA-17(M) AVHRR datasets
5 : * Author: Andrey Kiselev, dron@ak4719.spb.edu
6 : *
7 : * Some format info at: http://www.sat.dundee.ac.uk/noaa1b.html
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2002, Andrey Kiselev <dron@ak4719.spb.edu>
11 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
12 : *
13 : * ----------------------------------------------------------------------------
14 : * Lagrange interpolation suitable for NOAA level 1B file formats.
15 : * Submitted by Andrew Brooks <arb@sat.dundee.ac.uk>
16 : *
17 : * SPDX-License-Identifier: MIT
18 : ****************************************************************************/
19 :
20 : #include "cpl_string.h"
21 : #include "gdal_frmts.h"
22 : #include "gdal_pam.h"
23 : #include "ogr_srs_api.h"
24 :
25 : #include <algorithm>
26 :
27 : typedef enum
28 : { // File formats
29 : L1B_NONE, // Not a L1B format
30 : L1B_NOAA9, // NOAA-9/14
31 : L1B_NOAA15, // NOAA-15/METOP-2
32 : L1B_NOAA15_NOHDR // NOAA-15/METOP-2 without ARS header
33 : } L1BFileFormat;
34 :
35 : typedef enum
36 : { // Spacecrafts:
37 : TIROSN, // TIROS-N
38 : // NOAA are given a letter before launch and a number after launch
39 : NOAA6, // NOAA-6(A)
40 : NOAAB, // NOAA-B
41 : NOAA7, // NOAA-7(C)
42 : NOAA8, // NOAA-8(E)
43 : NOAA9_UNKNOWN, // Some NOAA-18 and NOAA-19 HRPT are recognized like that...
44 : NOAA9, // NOAA-9(F)
45 : NOAA10, // NOAA-10(G)
46 : NOAA11, // NOAA-11(H)
47 : NOAA12, // NOAA-12(D)
48 : NOAA13, // NOAA-13(I)
49 : NOAA14, // NOAA-14(J)
50 : NOAA15, // NOAA-15(K)
51 : NOAA16, // NOAA-16(L)
52 : NOAA17, // NOAA-17(M)
53 : NOAA18, // NOAA-18(N)
54 : NOAA19, // NOAA-19(N')
55 : // MetOp are given a number before launch and a letter after launch
56 : METOP2, // METOP-A(2)
57 : METOP1, // METOP-B(1)
58 : METOP3, // METOP-C(3)
59 : } L1BSpaceCraftdID;
60 :
61 : typedef enum
62 : { // Product types
63 : HRPT,
64 : LAC,
65 : GAC,
66 : FRAC
67 : } L1BProductType;
68 :
69 : typedef enum
70 : { // Data format
71 : PACKED10BIT,
72 : UNPACKED8BIT,
73 : UNPACKED16BIT
74 : } L1BDataFormat;
75 :
76 : typedef enum
77 : { // Receiving stations names:
78 : DU, // Dundee, Scotland, UK
79 : GC, // Fairbanks, Alaska, USA (formerly Gilmore Creek)
80 : HO, // Honolulu, Hawaii, USA
81 : MO, // Monterey, California, USA
82 : WE, // Western Europe CDA, Lannion, France
83 : SO, // SOCC (Satellite Operations Control Center), Suitland, Maryland, USA
84 : WI, // Wallops Island, Virginia, USA
85 : SV, // Svalbard, Norway
86 : UNKNOWN_STATION
87 : } L1BReceivingStation;
88 :
89 : typedef enum
90 : { // Data processing centers:
91 : CMS, // Centre de Meteorologie Spatiale - Lannion, France
92 : DSS, // Dundee Satellite Receiving Station - Dundee, Scotland, UK
93 : NSS, // NOAA/NESDIS - Suitland, Maryland, USA
94 : UKM, // United Kingdom Meteorological Office - Bracknell, England, UK
95 : UNKNOWN_CENTER
96 : } L1BProcessingCenter;
97 :
98 : typedef enum
99 : { // AVHRR Earth location indication
100 : ASCEND,
101 : DESCEND
102 : } L1BAscendOrDescend;
103 :
104 : /************************************************************************/
105 : /* AVHRR band widths */
106 : /************************************************************************/
107 :
108 : static const char *const apszBandDesc[] = {
109 : // NOAA-7 -- METOP-2 channels
110 : "AVHRR Channel 1: 0.58 micrometers -- 0.68 micrometers",
111 : "AVHRR Channel 2: 0.725 micrometers -- 1.10 micrometers",
112 : "AVHRR Channel 3: 3.55 micrometers -- 3.93 micrometers",
113 : "AVHRR Channel 4: 10.3 micrometers -- 11.3 micrometers",
114 : "AVHRR Channel 5: 11.5 micrometers -- 12.5 micrometers", // not in
115 : // NOAA-6,-8,-10
116 : // NOAA-13
117 : "AVHRR Channel 5: 11.4 micrometers -- 12.4 micrometers",
118 : // NOAA-15 -- METOP-2
119 : "AVHRR Channel 3A: 1.58 micrometers -- 1.64 micrometers",
120 : "AVHRR Channel 3B: 3.55 micrometers -- 3.93 micrometers"};
121 :
122 : /************************************************************************/
123 : /* L1B file format related constants */
124 : /************************************************************************/
125 :
126 : // Length of the string containing dataset name
127 : #define L1B_DATASET_NAME_SIZE 42
128 : #define L1B_NOAA9_HEADER_SIZE 122 // Terabit memory (TBM) header length
129 : #define L1B_NOAA9_HDR_NAME_OFF 30 // Dataset name offset
130 : #define L1B_NOAA9_HDR_SRC_OFF 70 // Receiving station name offset
131 : #define L1B_NOAA9_HDR_CHAN_OFF 97 // Selected channels map offset
132 : #define L1B_NOAA9_HDR_CHAN_SIZE 20 // Length of selected channels map
133 : #define L1B_NOAA9_HDR_WORD_OFF 117 // Sensor data word size offset
134 :
135 : // Archive Retrieval System (ARS) header
136 : #define L1B_NOAA15_HEADER_SIZE 512
137 : #define L1B_NOAA15_HDR_CHAN_OFF 97 // Selected channels map offset
138 : #define L1B_NOAA15_HDR_CHAN_SIZE 20 // Length of selected channels map
139 : #define L1B_NOAA15_HDR_WORD_OFF 117 // Sensor data word size offset
140 :
141 : // Length of header record filled with the data
142 : #define L1B_NOAA9_HDR_REC_SIZE 146
143 : #define L1B_NOAA9_HDR_REC_ID_OFF 0 // Spacecraft ID offset
144 : #define L1B_NOAA9_HDR_REC_PROD_OFF 1 // Data type offset
145 : #define L1B_NOAA9_HDR_REC_DSTAT_OFF 34 // DACS status offset
146 :
147 : /* See
148 : * http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c8/sec83132-2.htm */
149 : // Length of header record filled with the data
150 : #define L1B_NOAA15_HDR_REC_SIZE 992
151 : #define L1B_NOAA15_HDR_REC_SITE_OFF 0 // Dataset creation site ID offset
152 : #define L1B_NOAA15_HDR_REC_FORMAT_VERSION_OFF \
153 : 4 // NOAA Level 1b Format Version Number
154 : #define L1B_NOAA15_HDR_REC_FORMAT_VERSION_YEAR_OFF \
155 : 6 // Level 1b Format Version Year (e.g., 1999)
156 : #define L1B_NOAA15_HDR_REC_FORMAT_VERSION_DAY_OFF \
157 : 8 // Level 1b Format Version Day of Year (e.g., 365)
158 : #define L1B_NOAA15_HDR_REC_LOGICAL_REC_LENGTH_OFF \
159 : 10 // Logical Record Length of source Level 1b data set prior to processing
160 : #define L1B_NOAA15_HDR_REC_BLOCK_SIZE_OFF \
161 : 12 // Block Size of source Level 1b data set prior to processing
162 : #define L1B_NOAA15_HDR_REC_HDR_REC_COUNT_OFF \
163 : 14 // Count of Header Records in this Data Set
164 : #define L1B_NOAA15_HDR_REC_NAME_OFF 22 // Dataset name
165 : #define L1B_NOAA15_HDR_REC_ID_OFF 72 // Spacecraft ID offset
166 : #define L1B_NOAA15_HDR_REC_PROD_OFF 76 // Data type offset
167 : #define L1B_NOAA15_HDR_REC_STAT_OFF 116 // Instrument status offset
168 : #define L1B_NOAA15_HDR_REC_DATA_RECORD_COUNT_OFF 128
169 : #define L1B_NOAA15_HDR_REC_CALIBRATED_SCANLINE_COUNT_OFF 130
170 : #define L1B_NOAA15_HDR_REC_MISSING_SCANLINE_COUNT_OFF 132
171 : #define L1B_NOAA15_HDR_REC_SRC_OFF 154 // Receiving station name offset
172 : #define L1B_NOAA15_HDR_REC_ELLIPSOID_OFF 328
173 :
174 : /* This only apply if L1B_HIGH_GCP_DENSITY is explicitly set to NO */
175 : /* otherwise we will report more GCPs */
176 : #define DESIRED_GCPS_PER_LINE 11
177 : #define DESIRED_LINES_OF_GCPS 20
178 :
179 : // Fixed values used to scale GCPs coordinates in AVHRR records
180 : #define L1B_NOAA9_GCP_SCALE 128.0
181 : #define L1B_NOAA15_GCP_SCALE 10000.0
182 :
183 : /************************************************************************/
184 : /* ==================================================================== */
185 : /* TimeCode (helper class) */
186 : /* ==================================================================== */
187 : /************************************************************************/
188 :
189 : #define L1B_TIMECODE_LENGTH 100
190 :
191 : class TimeCode
192 : {
193 : long lYear;
194 : long lDay;
195 : long lMillisecond;
196 : char szString[L1B_TIMECODE_LENGTH];
197 :
198 : public:
199 6 : TimeCode() : lYear(0), lDay(0), lMillisecond(0)
200 : {
201 6 : memset(szString, 0, sizeof(szString));
202 6 : }
203 :
204 6 : void SetYear(long year)
205 : {
206 6 : lYear = year;
207 6 : }
208 :
209 6 : void SetDay(long day)
210 : {
211 6 : lDay = day;
212 6 : }
213 :
214 6 : void SetMillisecond(long millisecond)
215 : {
216 6 : lMillisecond = millisecond;
217 6 : }
218 :
219 0 : long GetYear() const
220 : {
221 0 : return lYear;
222 : }
223 :
224 0 : long GetDay() const
225 : {
226 0 : return lDay;
227 : }
228 :
229 0 : long GetMillisecond() const
230 : {
231 0 : return lMillisecond;
232 : }
233 :
234 6 : char *PrintTime()
235 : {
236 6 : snprintf(szString, L1B_TIMECODE_LENGTH,
237 : "year: %ld, day: %ld, millisecond: %ld", lYear, lDay,
238 : lMillisecond);
239 6 : return szString;
240 : }
241 : };
242 :
243 : #undef L1B_TIMECODE_LENGTH
244 :
245 : /************************************************************************/
246 : /* ==================================================================== */
247 : /* L1BDataset */
248 : /* ==================================================================== */
249 : /************************************************************************/
250 : class L1BGeolocDataset;
251 : class L1BGeolocRasterBand;
252 : class L1BSolarZenithAnglesDataset;
253 : class L1BSolarZenithAnglesRasterBand;
254 : class L1BNOAA15AnglesDataset;
255 : class L1BNOAA15AnglesRasterBand;
256 : class L1BCloudsDataset;
257 : class L1BCloudsRasterBand;
258 :
259 : class L1BDataset final : public GDALPamDataset
260 : {
261 : friend class L1BRasterBand;
262 : friend class L1BMaskBand;
263 : friend class L1BGeolocDataset;
264 : friend class L1BGeolocRasterBand;
265 : friend class L1BSolarZenithAnglesDataset;
266 : friend class L1BSolarZenithAnglesRasterBand;
267 : friend class L1BNOAA15AnglesDataset;
268 : friend class L1BNOAA15AnglesRasterBand;
269 : friend class L1BCloudsDataset;
270 : friend class L1BCloudsRasterBand;
271 :
272 : // char pszRevolution[6]; // Five-digit number identifying spacecraft
273 : // revolution
274 : L1BReceivingStation eSource; // Source of data (receiving station name)
275 : L1BProcessingCenter eProcCenter; // Data processing center
276 : TimeCode sStartTime;
277 : TimeCode sStopTime;
278 :
279 : int bHighGCPDensityStrategy;
280 : GDAL_GCP *pasGCPList;
281 : int nGCPCount;
282 : int iGCPOffset;
283 : int iGCPCodeOffset;
284 : int iCLAVRStart;
285 : int nGCPsPerLine;
286 : int eLocationIndicator;
287 : int iGCPStart;
288 : int iGCPStep;
289 :
290 : L1BFileFormat eL1BFormat;
291 : int nBufferSize;
292 : L1BSpaceCraftdID eSpacecraftID;
293 : L1BProductType eProductType; // LAC, GAC, HRPT, FRAC
294 : L1BDataFormat iDataFormat; // 10-bit packed or 16-bit unpacked
295 : int nRecordDataStart;
296 : int nRecordDataEnd;
297 : int nDataStartOffset;
298 : int nRecordSize;
299 : int nRecordSizeFromHeader;
300 : GUInt32 iInstrumentStatus;
301 : GUInt32 iChannelsMask;
302 :
303 : OGRSpatialReference m_oGCPSRS{};
304 :
305 : VSILFILE *fp;
306 :
307 : int bGuessDataFormat;
308 :
309 : int bByteSwap;
310 :
311 : int bExposeMaskBand;
312 : GDALRasterBand *poMaskBand;
313 :
314 : void ProcessRecordHeaders();
315 : int FetchGCPs(GDAL_GCP *, GByte *, int) const;
316 : static void FetchNOAA9TimeCode(TimeCode *, const GByte *, int *);
317 : void FetchNOAA15TimeCode(TimeCode *, const GByte *, int *) const;
318 : void FetchTimeCode(TimeCode *psTime, const void *pRecordHeader,
319 : int *peLocationIndicator) const;
320 : CPLErr ProcessDatasetHeader(const char *pszFilename);
321 : int ComputeFileOffsets();
322 :
323 : void FetchMetadata();
324 : void FetchMetadataNOAA15();
325 :
326 : vsi_l_offset GetLineOffset(int nBlockYOff) const;
327 :
328 : GUInt16 GetUInt16(const void *pabyData) const;
329 : GInt16 GetInt16(const void *pabyData) const;
330 : GUInt32 GetUInt32(const void *pabyData) const;
331 : GInt32 GetInt32(const void *pabyData) const;
332 :
333 : static L1BFileFormat DetectFormat(const char *pszFilename,
334 : const GByte *pabyHeader,
335 : int nHeaderBytes);
336 :
337 : public:
338 : explicit L1BDataset(L1BFileFormat);
339 : virtual ~L1BDataset();
340 :
341 : virtual int GetGCPCount() override;
342 : const OGRSpatialReference *GetGCPSpatialRef() const override;
343 : virtual const GDAL_GCP *GetGCPs() override;
344 :
345 : static int Identify(GDALOpenInfo *);
346 : static GDALDataset *Open(GDALOpenInfo *);
347 : };
348 :
349 : /************************************************************************/
350 : /* ==================================================================== */
351 : /* L1BRasterBand */
352 : /* ==================================================================== */
353 : /************************************************************************/
354 :
355 : class L1BRasterBand final : public GDALPamRasterBand
356 : {
357 : friend class L1BDataset;
358 :
359 : public:
360 : L1BRasterBand(L1BDataset *, int);
361 :
362 : // virtual double GetNoDataValue( int *pbSuccess = NULL );
363 : virtual CPLErr IReadBlock(int, int, void *) override;
364 : virtual GDALRasterBand *GetMaskBand() override;
365 : virtual int GetMaskFlags() override;
366 : };
367 :
368 : /************************************************************************/
369 : /* ==================================================================== */
370 : /* L1BMaskBand */
371 : /* ==================================================================== */
372 : /************************************************************************/
373 :
374 : class L1BMaskBand final : public GDALPamRasterBand
375 : {
376 : friend class L1BDataset;
377 :
378 : public:
379 : explicit L1BMaskBand(L1BDataset *);
380 :
381 : virtual CPLErr IReadBlock(int, int, void *) override;
382 : };
383 :
384 : /************************************************************************/
385 : /* L1BMaskBand() */
386 : /************************************************************************/
387 :
388 1 : L1BMaskBand::L1BMaskBand(L1BDataset *poDSIn)
389 : {
390 1 : CPLAssert(poDSIn->eL1BFormat == L1B_NOAA15 ||
391 : poDSIn->eL1BFormat == L1B_NOAA15_NOHDR);
392 :
393 1 : poDS = poDSIn;
394 1 : eDataType = GDT_Byte;
395 :
396 1 : nRasterXSize = poDS->GetRasterXSize();
397 1 : nRasterYSize = poDS->GetRasterYSize();
398 1 : nBlockXSize = poDS->GetRasterXSize();
399 1 : nBlockYSize = 1;
400 1 : }
401 :
402 : /************************************************************************/
403 : /* IReadBlock() */
404 : /************************************************************************/
405 :
406 2 : CPLErr L1BMaskBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
407 : void *pImage)
408 : {
409 2 : L1BDataset *poGDS = cpl::down_cast<L1BDataset *>(poDS);
410 :
411 2 : CPL_IGNORE_RET_VAL(
412 2 : VSIFSeekL(poGDS->fp, poGDS->GetLineOffset(nBlockYOff) + 24, SEEK_SET));
413 :
414 : GByte abyData[4];
415 2 : CPL_IGNORE_RET_VAL(VSIFReadL(abyData, 1, 4, poGDS->fp));
416 2 : GUInt32 n32 = poGDS->GetUInt32(abyData);
417 :
418 2 : if ((n32 >> 31) != 0) /* fatal flag */
419 1 : memset(pImage, 0, nBlockXSize);
420 : else
421 1 : memset(pImage, 255, nBlockXSize);
422 :
423 2 : return CE_None;
424 : }
425 :
426 : /************************************************************************/
427 : /* L1BRasterBand() */
428 : /************************************************************************/
429 :
430 7 : L1BRasterBand::L1BRasterBand(L1BDataset *poDSIn, int nBandIn)
431 :
432 : {
433 7 : poDS = poDSIn;
434 7 : nBand = nBandIn;
435 7 : eDataType = GDT_UInt16;
436 :
437 7 : nBlockXSize = poDS->GetRasterXSize();
438 7 : nBlockYSize = 1;
439 7 : }
440 :
441 : /************************************************************************/
442 : /* GetMaskBand() */
443 : /************************************************************************/
444 :
445 1 : GDALRasterBand *L1BRasterBand::GetMaskBand()
446 : {
447 1 : L1BDataset *poGDS = cpl::down_cast<L1BDataset *>(poDS);
448 1 : if (poGDS->poMaskBand)
449 1 : return poGDS->poMaskBand;
450 0 : return GDALPamRasterBand::GetMaskBand();
451 : }
452 :
453 : /************************************************************************/
454 : /* GetMaskFlags() */
455 : /************************************************************************/
456 :
457 1 : int L1BRasterBand::GetMaskFlags()
458 : {
459 1 : L1BDataset *poGDS = cpl::down_cast<L1BDataset *>(poDS);
460 1 : if (poGDS->poMaskBand)
461 1 : return GMF_PER_DATASET;
462 0 : return GDALPamRasterBand::GetMaskFlags();
463 : }
464 :
465 : /************************************************************************/
466 : /* IReadBlock() */
467 : /************************************************************************/
468 :
469 2 : CPLErr L1BRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
470 : void *pImage)
471 : {
472 2 : L1BDataset *poGDS = cpl::down_cast<L1BDataset *>(poDS);
473 :
474 : /* -------------------------------------------------------------------- */
475 : /* Seek to data. */
476 : /* -------------------------------------------------------------------- */
477 2 : CPL_IGNORE_RET_VAL(
478 2 : VSIFSeekL(poGDS->fp, poGDS->GetLineOffset(nBlockYOff), SEEK_SET));
479 :
480 : /* -------------------------------------------------------------------- */
481 : /* Read data into the buffer. */
482 : /* -------------------------------------------------------------------- */
483 2 : GUInt16 *iScan = nullptr; // Unpacked 16-bit scanline buffer
484 : int i, j;
485 :
486 2 : switch (poGDS->iDataFormat)
487 : {
488 0 : case PACKED10BIT:
489 : {
490 : // Read packed scanline
491 0 : GUInt32 *iRawScan = (GUInt32 *)CPLMalloc(poGDS->nRecordSize);
492 0 : CPL_IGNORE_RET_VAL(
493 0 : VSIFReadL(iRawScan, 1, poGDS->nRecordSize, poGDS->fp));
494 :
495 0 : iScan = (GUInt16 *)CPLMalloc(poGDS->nBufferSize);
496 0 : j = 0;
497 0 : for (i = poGDS->nRecordDataStart / (int)sizeof(iRawScan[0]);
498 0 : i < poGDS->nRecordDataEnd / (int)sizeof(iRawScan[0]); i++)
499 : {
500 0 : GUInt32 iWord1 = poGDS->GetUInt32(&iRawScan[i]);
501 0 : GUInt32 iWord2 = iWord1 & 0x3FF00000;
502 :
503 0 : iScan[j++] = (GUInt16)(iWord2 >> 20);
504 0 : iWord2 = iWord1 & 0x000FFC00;
505 0 : iScan[j++] = (GUInt16)(iWord2 >> 10);
506 0 : iScan[j++] = (GUInt16)(iWord1 & 0x000003FF);
507 : }
508 0 : CPLFree(iRawScan);
509 : }
510 0 : break;
511 2 : case UNPACKED16BIT:
512 : {
513 : // Read unpacked scanline
514 2 : GUInt16 *iRawScan = (GUInt16 *)CPLMalloc(poGDS->nRecordSize);
515 2 : CPL_IGNORE_RET_VAL(
516 2 : VSIFReadL(iRawScan, 1, poGDS->nRecordSize, poGDS->fp));
517 :
518 4 : iScan = (GUInt16 *)CPLMalloc(
519 2 : sizeof(GUInt16) * poGDS->GetRasterXSize() * poGDS->nBands);
520 20482 : for (i = 0; i < poGDS->GetRasterXSize() * poGDS->nBands; i++)
521 : {
522 20480 : iScan[i] =
523 20480 : poGDS->GetUInt16(&iRawScan[poGDS->nRecordDataStart /
524 20480 : (int)sizeof(iRawScan[0]) +
525 20480 : i]);
526 : }
527 2 : CPLFree(iRawScan);
528 : }
529 2 : break;
530 0 : case UNPACKED8BIT:
531 : {
532 : // Read 8-bit unpacked scanline
533 0 : GByte *byRawScan = (GByte *)CPLMalloc(poGDS->nRecordSize);
534 0 : CPL_IGNORE_RET_VAL(
535 0 : VSIFReadL(byRawScan, 1, poGDS->nRecordSize, poGDS->fp));
536 :
537 0 : iScan = (GUInt16 *)CPLMalloc(
538 0 : sizeof(GUInt16) * poGDS->GetRasterXSize() * poGDS->nBands);
539 0 : for (i = 0; i < poGDS->GetRasterXSize() * poGDS->nBands; i++)
540 0 : iScan[i] = byRawScan[poGDS->nRecordDataStart /
541 0 : (int)sizeof(byRawScan[0]) +
542 0 : i];
543 0 : CPLFree(byRawScan);
544 : }
545 0 : break;
546 0 : default: // NOTREACHED
547 0 : break;
548 : }
549 :
550 2 : int nBlockSize = nBlockXSize * nBlockYSize;
551 2 : if (poGDS->eLocationIndicator == DESCEND)
552 : {
553 0 : for (i = 0, j = 0; i < nBlockSize; i++)
554 : {
555 0 : ((GUInt16 *)pImage)[i] = iScan[j + nBand - 1];
556 0 : j += poGDS->nBands;
557 : }
558 : }
559 : else
560 : {
561 4098 : for (i = nBlockSize - 1, j = 0; i >= 0; i--)
562 : {
563 4096 : ((GUInt16 *)pImage)[i] = iScan[j + nBand - 1];
564 4096 : j += poGDS->nBands;
565 : }
566 : }
567 :
568 2 : CPLFree(iScan);
569 2 : return CE_None;
570 : }
571 :
572 : /************************************************************************/
573 : /* L1BDataset() */
574 : /************************************************************************/
575 :
576 3 : L1BDataset::L1BDataset(L1BFileFormat eL1BFormatIn)
577 : : eSource(UNKNOWN_STATION), eProcCenter(UNKNOWN_CENTER),
578 : // sStartTime
579 : // sStopTime
580 : bHighGCPDensityStrategy(
581 3 : CPLTestBool(CPLGetConfigOption("L1B_HIGH_GCP_DENSITY", "TRUE"))),
582 : pasGCPList(nullptr), nGCPCount(0), iGCPOffset(0), iGCPCodeOffset(0),
583 : iCLAVRStart(0), nGCPsPerLine(0),
584 : eLocationIndicator(DESCEND), // XXX: should be initialized
585 : iGCPStart(0), iGCPStep(0), eL1BFormat(eL1BFormatIn), nBufferSize(0),
586 : eSpacecraftID(TIROSN), eProductType(HRPT), iDataFormat(PACKED10BIT),
587 : nRecordDataStart(0), nRecordDataEnd(0), nDataStartOffset(0),
588 : nRecordSize(0), nRecordSizeFromHeader(0), iInstrumentStatus(0),
589 : iChannelsMask(0), fp(nullptr), bGuessDataFormat(FALSE),
590 : // L1B is normally big-endian ordered, so byte-swap on little-endian CPU.
591 6 : bByteSwap(CPL_IS_LSB), bExposeMaskBand(FALSE), poMaskBand(nullptr)
592 : {
593 3 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
594 3 : m_oGCPSRS.importFromWkt(
595 : "GEOGCS[\"WGS 72\",DATUM[\"WGS_1972\","
596 : "SPHEROID[\"WGS 72\",6378135,298.26,AUTHORITY[\"EPSG\",7043]],"
597 : "TOWGS84[0,0,4.5,0,0,0.554,0.2263],AUTHORITY[\"EPSG\",6322]],"
598 : "PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",8901]],"
599 : "UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",9108]],"
600 : "AUTHORITY[\"EPSG\",4322]]");
601 3 : }
602 :
603 : /************************************************************************/
604 : /* ~L1BDataset() */
605 : /************************************************************************/
606 :
607 6 : L1BDataset::~L1BDataset()
608 :
609 : {
610 3 : FlushCache(true);
611 :
612 3 : if (nGCPCount > 0)
613 : {
614 1 : GDALDeinitGCPs(nGCPCount, pasGCPList);
615 1 : CPLFree(pasGCPList);
616 : }
617 3 : if (fp != nullptr)
618 3 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
619 3 : delete poMaskBand;
620 6 : }
621 :
622 : /************************************************************************/
623 : /* GetLineOffset() */
624 : /************************************************************************/
625 :
626 4 : vsi_l_offset L1BDataset::GetLineOffset(int nBlockYOff) const
627 : {
628 4 : return (eLocationIndicator == DESCEND)
629 4 : ? nDataStartOffset + (vsi_l_offset)nBlockYOff * nRecordSize
630 4 : : nDataStartOffset +
631 4 : (vsi_l_offset)(nRasterYSize - nBlockYOff - 1) *
632 4 : nRecordSize;
633 : }
634 :
635 : /************************************************************************/
636 : /* GetGCPCount() */
637 : /************************************************************************/
638 :
639 0 : int L1BDataset::GetGCPCount()
640 :
641 : {
642 0 : return nGCPCount;
643 : }
644 :
645 : /************************************************************************/
646 : /* GetGCPSpatialRef() */
647 : /************************************************************************/
648 :
649 1 : const OGRSpatialReference *L1BDataset::GetGCPSpatialRef() const
650 :
651 : {
652 1 : if (nGCPCount > 0 && !m_oGCPSRS.IsEmpty())
653 1 : return &m_oGCPSRS;
654 : else
655 0 : return nullptr;
656 : }
657 :
658 : /************************************************************************/
659 : /* GetGCPs() */
660 : /************************************************************************/
661 :
662 0 : const GDAL_GCP *L1BDataset::GetGCPs()
663 : {
664 0 : return pasGCPList;
665 : }
666 :
667 : /************************************************************************/
668 : /* Byte swapping helpers */
669 : /************************************************************************/
670 :
671 20500 : GUInt16 L1BDataset::GetUInt16(const void *pabyData) const
672 : {
673 : GUInt16 iTemp;
674 20500 : memcpy(&iTemp, pabyData, 2);
675 20500 : if (bByteSwap)
676 3 : return CPL_SWAP16(iTemp);
677 20497 : return iTemp;
678 : }
679 :
680 0 : GInt16 L1BDataset::GetInt16(const void *pabyData) const
681 : {
682 : GInt16 iTemp;
683 0 : memcpy(&iTemp, pabyData, 2);
684 0 : if (bByteSwap)
685 0 : return CPL_SWAP16(iTemp);
686 0 : return iTemp;
687 : }
688 :
689 5 : GUInt32 L1BDataset::GetUInt32(const void *pabyData) const
690 : {
691 : GUInt32 lTemp;
692 5 : memcpy(&lTemp, pabyData, 4);
693 5 : if (bByteSwap)
694 0 : return CPL_SWAP32(lTemp);
695 5 : return lTemp;
696 : }
697 :
698 204 : GInt32 L1BDataset::GetInt32(const void *pabyData) const
699 : {
700 : GInt32 lTemp;
701 204 : memcpy(&lTemp, pabyData, 4);
702 204 : if (bByteSwap)
703 0 : return CPL_SWAP32(lTemp);
704 204 : return lTemp;
705 : }
706 :
707 : /************************************************************************/
708 : /* Fetch timecode from the record header (NOAA9-NOAA14 version) */
709 : /************************************************************************/
710 :
711 4 : void L1BDataset::FetchNOAA9TimeCode(TimeCode *psTime,
712 : const GByte *piRecordHeader,
713 : int *peLocationIndicator)
714 : {
715 : GUInt32 lTemp;
716 :
717 4 : lTemp = ((piRecordHeader[2] >> 1) & 0x7F);
718 8 : psTime->SetYear((lTemp > 77)
719 2 : ? (lTemp + 1900)
720 2 : : (lTemp + 2000)); // Avoid `Year 2000' problem
721 4 : psTime->SetDay((GUInt32)(piRecordHeader[2] & 0x01) << 8 |
722 4 : (GUInt32)piRecordHeader[3]);
723 4 : psTime->SetMillisecond(((GUInt32)(piRecordHeader[4] & 0x07) << 24) |
724 4 : ((GUInt32)piRecordHeader[5] << 16) |
725 4 : ((GUInt32)piRecordHeader[6] << 8) |
726 4 : (GUInt32)piRecordHeader[7]);
727 4 : if (peLocationIndicator)
728 : {
729 2 : *peLocationIndicator =
730 2 : ((piRecordHeader[8] & 0x02) == 0) ? ASCEND : DESCEND;
731 : }
732 4 : }
733 :
734 : /************************************************************************/
735 : /* Fetch timecode from the record header (NOAA15-METOP2 version) */
736 : /************************************************************************/
737 :
738 2 : void L1BDataset::FetchNOAA15TimeCode(TimeCode *psTime,
739 : const GByte *pabyRecordHeader,
740 : int *peLocationIndicator) const
741 : {
742 2 : psTime->SetYear(GetUInt16(pabyRecordHeader + 2));
743 2 : psTime->SetDay(GetUInt16(pabyRecordHeader + 4));
744 2 : psTime->SetMillisecond(GetUInt32(pabyRecordHeader + 8));
745 2 : if (peLocationIndicator)
746 : {
747 : // FIXME: hemisphere
748 1 : *peLocationIndicator =
749 1 : ((GetUInt16(pabyRecordHeader + 12) & 0x8000) == 0) ? ASCEND
750 : : DESCEND;
751 : }
752 2 : }
753 :
754 : /************************************************************************/
755 : /* FetchTimeCode() */
756 : /************************************************************************/
757 :
758 6 : void L1BDataset::FetchTimeCode(TimeCode *psTime, const void *pRecordHeader,
759 : int *peLocationIndicator) const
760 : {
761 6 : if (eSpacecraftID <= NOAA14)
762 : {
763 4 : FetchNOAA9TimeCode(psTime, (const GByte *)pRecordHeader,
764 : peLocationIndicator);
765 : }
766 : else
767 : {
768 2 : FetchNOAA15TimeCode(psTime, (const GByte *)pRecordHeader,
769 : peLocationIndicator);
770 : }
771 6 : }
772 :
773 : /************************************************************************/
774 : /* Fetch GCPs from the individual scanlines */
775 : /************************************************************************/
776 :
777 2 : int L1BDataset::FetchGCPs(GDAL_GCP *pasGCPListRow, GByte *pabyRecordHeader,
778 : int iLine) const
779 : {
780 : // LAC and HRPT GCPs are tied to the center of pixel,
781 : // GAC ones are slightly displaced.
782 2 : double dfDelta = (eProductType == GAC) ? 0.9 : 0.5;
783 4 : double dfPixel = (eLocationIndicator == DESCEND)
784 2 : ? iGCPStart + dfDelta
785 2 : : (nRasterXSize - (iGCPStart + dfDelta));
786 :
787 : int nGCPs;
788 2 : if (eSpacecraftID <= NOAA14)
789 : {
790 : // NOAA9-NOAA14 records have an indicator of number of working GCPs.
791 : // Number of good GCPs may be smaller than the total amount of points.
792 0 : nGCPs = (*(pabyRecordHeader + iGCPCodeOffset) < nGCPsPerLine)
793 : ? *(pabyRecordHeader + iGCPCodeOffset)
794 : : nGCPsPerLine;
795 : #ifdef DEBUG_VERBOSE
796 : CPLDebug("L1B", "iGCPCodeOffset=%d, nGCPsPerLine=%d, nGoodGCPs=%d",
797 : iGCPCodeOffset, nGCPsPerLine, nGCPs);
798 : #endif
799 : }
800 : else
801 2 : nGCPs = nGCPsPerLine;
802 :
803 2 : pabyRecordHeader += iGCPOffset;
804 :
805 2 : int nGCPCountRow = 0;
806 104 : while (nGCPs--)
807 : {
808 102 : if (eSpacecraftID <= NOAA14)
809 : {
810 0 : GInt16 nRawY = GetInt16(pabyRecordHeader);
811 0 : pabyRecordHeader += sizeof(GInt16);
812 0 : GInt16 nRawX = GetInt16(pabyRecordHeader);
813 0 : pabyRecordHeader += sizeof(GInt16);
814 :
815 0 : pasGCPListRow[nGCPCountRow].dfGCPY = nRawY / L1B_NOAA9_GCP_SCALE;
816 0 : pasGCPListRow[nGCPCountRow].dfGCPX = nRawX / L1B_NOAA9_GCP_SCALE;
817 : }
818 : else
819 : {
820 102 : GInt32 nRawY = GetInt32(pabyRecordHeader);
821 102 : pabyRecordHeader += sizeof(GInt32);
822 102 : GInt32 nRawX = GetInt32(pabyRecordHeader);
823 102 : pabyRecordHeader += sizeof(GInt32);
824 :
825 102 : pasGCPListRow[nGCPCountRow].dfGCPY = nRawY / L1B_NOAA15_GCP_SCALE;
826 102 : pasGCPListRow[nGCPCountRow].dfGCPX = nRawX / L1B_NOAA15_GCP_SCALE;
827 : }
828 :
829 102 : if (pasGCPListRow[nGCPCountRow].dfGCPX < -180 ||
830 102 : pasGCPListRow[nGCPCountRow].dfGCPX > 180 ||
831 102 : pasGCPListRow[nGCPCountRow].dfGCPY < -90 ||
832 102 : pasGCPListRow[nGCPCountRow].dfGCPY > 90)
833 0 : continue;
834 :
835 102 : pasGCPListRow[nGCPCountRow].dfGCPZ = 0.0;
836 102 : pasGCPListRow[nGCPCountRow].dfGCPPixel = dfPixel;
837 102 : dfPixel += (eLocationIndicator == DESCEND) ? iGCPStep : -iGCPStep;
838 102 : pasGCPListRow[nGCPCountRow].dfGCPLine =
839 102 : (double)((eLocationIndicator == DESCEND)
840 : ? iLine
841 102 : : nRasterYSize - iLine - 1) +
842 : 0.5;
843 102 : nGCPCountRow++;
844 : }
845 2 : return nGCPCountRow;
846 : }
847 :
848 : /************************************************************************/
849 : /* ProcessRecordHeaders() */
850 : /************************************************************************/
851 :
852 3 : void L1BDataset::ProcessRecordHeaders()
853 : {
854 3 : void *pRecordHeader = CPLCalloc(1, nRecordDataStart);
855 :
856 3 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, nDataStartOffset, SEEK_SET));
857 3 : CPL_IGNORE_RET_VAL(VSIFReadL(pRecordHeader, 1, nRecordDataStart, fp));
858 :
859 3 : FetchTimeCode(&sStartTime, pRecordHeader, &eLocationIndicator);
860 :
861 3 : CPL_IGNORE_RET_VAL(VSIFSeekL(
862 3 : fp, nDataStartOffset + (nRasterYSize - 1) * nRecordSize, SEEK_SET));
863 3 : CPL_IGNORE_RET_VAL(VSIFReadL(pRecordHeader, 1, nRecordDataStart, fp));
864 :
865 3 : FetchTimeCode(&sStopTime, pRecordHeader, nullptr);
866 :
867 : /* -------------------------------------------------------------------- */
868 : /* Pick a skip factor so that we will get roughly 20 lines */
869 : /* worth of GCPs. That should give respectable coverage on all */
870 : /* but the longest swaths. */
871 : /* -------------------------------------------------------------------- */
872 : int nTargetLines;
873 3 : double dfLineStep = 0.0;
874 :
875 3 : if (bHighGCPDensityStrategy)
876 : {
877 3 : if (nRasterYSize < nGCPsPerLine)
878 : {
879 3 : nTargetLines = nRasterYSize;
880 : }
881 : else
882 : {
883 : int nColStep;
884 0 : nColStep = nRasterXSize / nGCPsPerLine;
885 0 : if (nRasterYSize >= nRasterXSize)
886 : {
887 0 : dfLineStep = nColStep;
888 : }
889 : else
890 : {
891 0 : dfLineStep = nRasterYSize / nGCPsPerLine;
892 : }
893 0 : nTargetLines = static_cast<int>(nRasterYSize / dfLineStep);
894 : }
895 : }
896 : else
897 : {
898 0 : nTargetLines = std::min(DESIRED_LINES_OF_GCPS, nRasterYSize);
899 : }
900 3 : if (nTargetLines > 1)
901 1 : dfLineStep = 1.0 * (nRasterYSize - 1) / (nTargetLines - 1);
902 :
903 : /* -------------------------------------------------------------------- */
904 : /* Initialize the GCP list. */
905 : /* -------------------------------------------------------------------- */
906 3 : const int nExpectedGCPs = nTargetLines * nGCPsPerLine;
907 3 : if (nExpectedGCPs > 0)
908 : {
909 1 : pasGCPList =
910 1 : (GDAL_GCP *)VSI_CALLOC_VERBOSE(nExpectedGCPs, sizeof(GDAL_GCP));
911 1 : if (pasGCPList == nullptr)
912 : {
913 0 : CPLFree(pRecordHeader);
914 0 : return;
915 : }
916 1 : GDALInitGCPs(nExpectedGCPs, pasGCPList);
917 : }
918 :
919 : /* -------------------------------------------------------------------- */
920 : /* Fetch the GCPs for each selected line. We force the last */
921 : /* line sampled to be the last line in the dataset even if that */
922 : /* leaves a bigger than expected gap. */
923 : /* -------------------------------------------------------------------- */
924 : int iStep;
925 3 : int iPrevLine = -1;
926 :
927 5 : for (iStep = 0; iStep < nTargetLines; iStep++)
928 : {
929 : int iLine;
930 :
931 2 : if (iStep == nTargetLines - 1)
932 1 : iLine = nRasterYSize - 1;
933 : else
934 1 : iLine = (int)(dfLineStep * iStep);
935 2 : if (iLine == iPrevLine)
936 0 : continue;
937 2 : iPrevLine = iLine;
938 :
939 2 : CPL_IGNORE_RET_VAL(
940 2 : VSIFSeekL(fp, nDataStartOffset + iLine * nRecordSize, SEEK_SET));
941 2 : CPL_IGNORE_RET_VAL(VSIFReadL(pRecordHeader, 1, nRecordDataStart, fp));
942 :
943 : int nGCPsOnThisLine =
944 2 : FetchGCPs(pasGCPList + nGCPCount, (GByte *)pRecordHeader, iLine);
945 :
946 2 : if (!bHighGCPDensityStrategy)
947 : {
948 : /* --------------------------------------------------------------------
949 : */
950 : /* We don't really want too many GCPs per line. Downsample to
951 : */
952 : /* 11 per line. */
953 : /* --------------------------------------------------------------------
954 : */
955 :
956 : const int nDesiredGCPsPerLine =
957 0 : std::min(DESIRED_GCPS_PER_LINE, nGCPsOnThisLine);
958 0 : int nGCPStep =
959 : (nDesiredGCPsPerLine > 1)
960 0 : ? (nGCPsOnThisLine - 1) / (nDesiredGCPsPerLine - 1)
961 : : 1;
962 0 : int iSrcGCP = nGCPCount;
963 0 : int iDstGCP = nGCPCount;
964 :
965 0 : if (nGCPStep == 0)
966 0 : nGCPStep = 1;
967 :
968 0 : for (int iGCP = 0; iGCP < nDesiredGCPsPerLine; iGCP++)
969 : {
970 0 : if (iGCP == nDesiredGCPsPerLine - 1)
971 0 : iSrcGCP = nGCPCount + nGCPsOnThisLine - 1;
972 : else
973 0 : iSrcGCP += nGCPStep;
974 0 : iDstGCP++;
975 :
976 0 : pasGCPList[iDstGCP].dfGCPX = pasGCPList[iSrcGCP].dfGCPX;
977 0 : pasGCPList[iDstGCP].dfGCPY = pasGCPList[iSrcGCP].dfGCPY;
978 0 : pasGCPList[iDstGCP].dfGCPPixel = pasGCPList[iSrcGCP].dfGCPPixel;
979 0 : pasGCPList[iDstGCP].dfGCPLine = pasGCPList[iSrcGCP].dfGCPLine;
980 : }
981 :
982 0 : nGCPCount += nDesiredGCPsPerLine;
983 : }
984 : else
985 : {
986 2 : nGCPCount += nGCPsOnThisLine;
987 : }
988 : }
989 :
990 3 : if (nGCPCount < nExpectedGCPs)
991 : {
992 0 : GDALDeinitGCPs(nExpectedGCPs - nGCPCount, pasGCPList + nGCPCount);
993 0 : if (nGCPCount == 0)
994 : {
995 0 : CPLFree(pasGCPList);
996 0 : pasGCPList = nullptr;
997 : }
998 : }
999 :
1000 3 : CPLFree(pRecordHeader);
1001 :
1002 : /* -------------------------------------------------------------------- */
1003 : /* Set fetched information as metadata records */
1004 : /* -------------------------------------------------------------------- */
1005 : // Time of first scanline
1006 3 : SetMetadataItem("START", sStartTime.PrintTime());
1007 : // Time of last scanline
1008 3 : SetMetadataItem("STOP", sStopTime.PrintTime());
1009 : // AVHRR Earth location indication
1010 :
1011 3 : switch (eLocationIndicator)
1012 : {
1013 3 : case ASCEND:
1014 3 : SetMetadataItem("LOCATION", "Ascending");
1015 3 : break;
1016 0 : case DESCEND:
1017 : default:
1018 0 : SetMetadataItem("LOCATION", "Descending");
1019 0 : break;
1020 : }
1021 : }
1022 :
1023 : /************************************************************************/
1024 : /* FetchMetadata() */
1025 : /************************************************************************/
1026 :
1027 0 : void L1BDataset::FetchMetadata()
1028 : {
1029 0 : if (eL1BFormat != L1B_NOAA9)
1030 : {
1031 0 : FetchMetadataNOAA15();
1032 0 : return;
1033 : }
1034 :
1035 0 : std::string osDir = CPLGetConfigOption("L1B_METADATA_DIRECTORY", "");
1036 0 : if (osDir.empty())
1037 : {
1038 0 : osDir = CPLGetPathSafe(GetDescription());
1039 0 : if (osDir.empty())
1040 0 : osDir = ".";
1041 : }
1042 : const std::string osMetadataFile =
1043 0 : std::string(osDir)
1044 0 : .append("/")
1045 0 : .append(CPLGetFilename(GetDescription()))
1046 0 : .append("_metadata.csv");
1047 0 : CPL_IGNORE_RET_VAL(osDir);
1048 0 : VSILFILE *fpCSV = VSIFOpenL(osMetadataFile.c_str(), "wb");
1049 0 : if (fpCSV == nullptr)
1050 : {
1051 0 : CPLError(CE_Failure, CPLE_AppDefined,
1052 : "Cannot create metadata file : %s", osMetadataFile.c_str());
1053 0 : return;
1054 : }
1055 :
1056 0 : CPL_IGNORE_RET_VAL(
1057 0 : VSIFPrintfL(fpCSV, "SCANLINE,NBLOCKYOFF,YEAR,DAY,MS_IN_DAY,"));
1058 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1059 : fpCSV, "FATAL_FLAG,TIME_ERROR,DATA_GAP,DATA_JITTER,INSUFFICIENT_DATA_"
1060 : "FOR_CAL,NO_EARTH_LOCATION,DESCEND,P_N_STATUS,"));
1061 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1062 : fpCSV, "BIT_SYNC_STATUS,SYNC_ERROR,FRAME_SYNC_ERROR,FLYWHEELING,BIT_"
1063 : "SLIPPAGE,C3_SBBC,C4_SBBC,C5_SBBC,"));
1064 0 : CPL_IGNORE_RET_VAL(
1065 0 : VSIFPrintfL(fpCSV, "TIP_PARITY_FRAME_1,TIP_PARITY_FRAME_2,TIP_PARITY_"
1066 : "FRAME_3,TIP_PARITY_FRAME_4,TIP_PARITY_FRAME_5,"));
1067 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "SYNC_ERRORS,"));
1068 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1069 : fpCSV, "CAL_SLOPE_C1,CAL_INTERCEPT_C1,CAL_SLOPE_C2,CAL_INTERCEPT_C2,"
1070 : "CAL_SLOPE_C3,CAL_INTERCEPT_C3,CAL_SLOPE_C4,CAL_INTERCEPT_C4,"
1071 : "CAL_SLOPE_C5,CAL_INTERCEPT_C5,"));
1072 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "NUM_SOLZENANGLES_EARTHLOCPNTS"));
1073 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));
1074 :
1075 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(nRecordDataStart);
1076 :
1077 0 : for (int nBlockYOff = 0; nBlockYOff < nRasterYSize; nBlockYOff++)
1078 : {
1079 : /* --------------------------------------------------------------------
1080 : */
1081 : /* Seek to data. */
1082 : /* --------------------------------------------------------------------
1083 : */
1084 0 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, GetLineOffset(nBlockYOff), SEEK_SET));
1085 :
1086 0 : CPL_IGNORE_RET_VAL(
1087 0 : VSIFReadL(pabyRecordHeader, 1, nRecordDataStart, fp));
1088 :
1089 0 : GUInt16 nScanlineNumber = GetUInt16(pabyRecordHeader);
1090 :
1091 0 : TimeCode timeCode;
1092 0 : FetchTimeCode(&timeCode, pabyRecordHeader, nullptr);
1093 :
1094 0 : CPL_IGNORE_RET_VAL(
1095 0 : VSIFPrintfL(fpCSV, "%d,%d,%d,%d,%d,", nScanlineNumber, nBlockYOff,
1096 0 : (int)timeCode.GetYear(), (int)timeCode.GetDay(),
1097 0 : (int)timeCode.GetMillisecond()));
1098 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1099 0 : fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,", (pabyRecordHeader[8] >> 7) & 1,
1100 0 : (pabyRecordHeader[8] >> 6) & 1, (pabyRecordHeader[8] >> 5) & 1,
1101 0 : (pabyRecordHeader[8] >> 4) & 1, (pabyRecordHeader[8] >> 3) & 1,
1102 0 : (pabyRecordHeader[8] >> 2) & 1, (pabyRecordHeader[8] >> 1) & 1,
1103 0 : (pabyRecordHeader[8] >> 0) & 1));
1104 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1105 0 : fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,", (pabyRecordHeader[9] >> 7) & 1,
1106 0 : (pabyRecordHeader[9] >> 6) & 1, (pabyRecordHeader[9] >> 5) & 1,
1107 0 : (pabyRecordHeader[9] >> 4) & 1, (pabyRecordHeader[9] >> 3) & 1,
1108 0 : (pabyRecordHeader[9] >> 2) & 1, (pabyRecordHeader[9] >> 1) & 1,
1109 0 : (pabyRecordHeader[9] >> 0) & 1));
1110 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1111 0 : fpCSV, "%d,%d,%d,%d,%d,", (pabyRecordHeader[10] >> 7) & 1,
1112 0 : (pabyRecordHeader[10] >> 6) & 1, (pabyRecordHeader[10] >> 5) & 1,
1113 0 : (pabyRecordHeader[10] >> 4) & 1, (pabyRecordHeader[10] >> 3) & 1));
1114 0 : CPL_IGNORE_RET_VAL(
1115 0 : VSIFPrintfL(fpCSV, "%d,", pabyRecordHeader[11] >> 2));
1116 : GInt32 i32;
1117 0 : for (int i = 0; i < 10; i++)
1118 : {
1119 0 : i32 = GetInt32(pabyRecordHeader + 12 + 4 * i);
1120 : /* Scales :
1121 : * http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c3/sec3-3.htm
1122 : */
1123 0 : if ((i % 2) == 0)
1124 0 : CPL_IGNORE_RET_VAL(
1125 0 : VSIFPrintfL(fpCSV, "%f,", i32 / pow(2.0, 30.0)));
1126 : else
1127 0 : CPL_IGNORE_RET_VAL(
1128 0 : VSIFPrintfL(fpCSV, "%f,", i32 / pow(2.0, 22.0)));
1129 : }
1130 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d", pabyRecordHeader[52]));
1131 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));
1132 : }
1133 :
1134 0 : CPLFree(pabyRecordHeader);
1135 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpCSV));
1136 : }
1137 :
1138 : /************************************************************************/
1139 : /* FetchMetadataNOAA15() */
1140 : /************************************************************************/
1141 :
1142 0 : void L1BDataset::FetchMetadataNOAA15()
1143 : {
1144 : int i, j;
1145 0 : std::string osDir = CPLGetConfigOption("L1B_METADATA_DIRECTORY", "");
1146 0 : if (osDir.empty())
1147 : {
1148 0 : osDir = CPLGetPathSafe(GetDescription());
1149 0 : if (osDir.empty())
1150 0 : osDir = ".";
1151 : }
1152 : const std::string osMetadataFile =
1153 0 : std::string(osDir)
1154 0 : .append("/")
1155 0 : .append(CPLGetFilename(GetDescription()))
1156 0 : .append("_metadata.csv");
1157 0 : CPL_IGNORE_RET_VAL(osDir);
1158 0 : VSILFILE *fpCSV = VSIFOpenL(osMetadataFile.c_str(), "wb");
1159 0 : if (fpCSV == nullptr)
1160 : {
1161 0 : CPLError(CE_Failure, CPLE_AppDefined,
1162 : "Cannot create metadata file : %s", osMetadataFile.c_str());
1163 0 : return;
1164 : }
1165 :
1166 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1167 : fpCSV, "SCANLINE,NBLOCKYOFF,YEAR,DAY,MS_IN_DAY,SAT_CLOCK_DRIF_DELTA,"
1168 : "SOUTHBOUND,SCANTIME_CORRECTED,C3_SELECT,"));
1169 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1170 : fpCSV,
1171 : "FATAL_FLAG,TIME_ERROR,DATA_GAP,INSUFFICIENT_DATA_FOR_CAL,"
1172 : "NO_EARTH_LOCATION,FIRST_GOOD_TIME_AFTER_CLOCK_UPDATE,"
1173 : "INSTRUMENT_STATUS_CHANGED,SYNC_LOCK_DROPPED,"
1174 : "FRAME_SYNC_ERROR,FRAME_SYNC_DROPPED_LOCK,FLYWHEELING,"
1175 : "BIT_SLIPPAGE,TIP_PARITY_ERROR,REFLECTED_SUNLIGHT_C3B,"
1176 : "REFLECTED_SUNLIGHT_C4,REFLECTED_SUNLIGHT_C5,RESYNC,P_N_STATUS,"));
1177 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1178 : fpCSV, "BAD_TIME_CAN_BE_INFERRED,BAD_TIME_CANNOT_BE_INFERRED,"
1179 : "TIME_DISCONTINUITY,REPEAT_SCAN_TIME,"));
1180 0 : CPL_IGNORE_RET_VAL(
1181 0 : VSIFPrintfL(fpCSV, "UNCALIBRATED_BAD_TIME,CALIBRATED_FEWER_SCANLINES,"
1182 : "UNCALIBRATED_BAD_PRT,CALIBRATED_MARGINAL_PRT,"
1183 : "UNCALIBRATED_CHANNELS,"));
1184 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1185 : fpCSV, "NO_EARTH_LOC_BAD_TIME,EARTH_LOC_QUESTIONABLE_TIME,"
1186 : "EARTH_LOC_QUESTIONABLE,EARTH_LOC_VERY_QUESTIONABLE,"));
1187 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1188 : fpCSV,
1189 : "C3B_UNCALIBRATED,C3B_QUESTIONABLE,C3B_ALL_BLACKBODY,"
1190 : "C3B_ALL_SPACEVIEW,C3B_MARGINAL_BLACKBODY,C3B_MARGINAL_SPACEVIEW,"));
1191 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1192 : fpCSV,
1193 : "C4_UNCALIBRATED,C4_QUESTIONABLE,C4_ALL_BLACKBODY,"
1194 : "C4_ALL_SPACEVIEW,C4_MARGINAL_BLACKBODY,C4_MARGINAL_SPACEVIEW,"));
1195 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1196 : fpCSV,
1197 : "C5_UNCALIBRATED,C5_QUESTIONABLE,C5_ALL_BLACKBODY,"
1198 : "C5_ALL_SPACEVIEW,C5_MARGINAL_BLACKBODY,C5_MARGINAL_SPACEVIEW,"));
1199 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "BIT_ERRORS,"));
1200 0 : for (i = 0; i < 3; i++)
1201 : {
1202 0 : const char *pszChannel = (i == 0) ? "C1" : (i == 1) ? "C2" : "C3A";
1203 0 : for (j = 0; j < 3; j++)
1204 : {
1205 0 : const char *pszType = (j == 0) ? "OP"
1206 0 : : (j == 1) ? "TEST"
1207 : : "PRELAUNCH";
1208 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_SLOPE_1,",
1209 : pszType, pszChannel));
1210 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_INTERCEPT_1,",
1211 : pszType, pszChannel));
1212 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_SLOPE_2,",
1213 : pszType, pszChannel));
1214 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_INTERCEPT_2,",
1215 : pszType, pszChannel));
1216 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_INTERSECTION,",
1217 : pszType, pszChannel));
1218 : }
1219 : }
1220 0 : for (i = 0; i < 3; i++)
1221 : {
1222 0 : const char *pszChannel = (i == 0) ? "C3B" : (i == 1) ? "C4" : "C5";
1223 0 : for (j = 0; j < 2; j++)
1224 : {
1225 0 : const char *pszType = (j == 0) ? "OP" : "TEST";
1226 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "IR_%s_CAL_%s_COEFF_1,",
1227 : pszType, pszChannel));
1228 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "IR_%s_CAL_%s_COEFF_2,",
1229 : pszType, pszChannel));
1230 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "IR_%s_CAL_%s_COEFF_3,",
1231 : pszType, pszChannel));
1232 : }
1233 : }
1234 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1235 : fpCSV, "EARTH_LOC_CORR_TIP_EULER,EARTH_LOC_IND,"
1236 : "SPACECRAFT_ATT_CTRL,ATT_SMODE,ATT_PASSIVE_WHEEL_TEST,"
1237 : "TIME_TIP_EULER,TIP_EULER_ROLL,TIP_EULER_PITCH,TIP_EULER_YAW,"
1238 : "SPACECRAFT_ALT"));
1239 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));
1240 :
1241 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(nRecordDataStart);
1242 : GInt16 i16;
1243 : GUInt16 n16;
1244 : GInt32 i32;
1245 : GUInt32 n32;
1246 :
1247 0 : for (int nBlockYOff = 0; nBlockYOff < nRasterYSize; nBlockYOff++)
1248 : {
1249 : /* --------------------------------------------------------------------
1250 : */
1251 : /* Seek to data. */
1252 : /* --------------------------------------------------------------------
1253 : */
1254 0 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, GetLineOffset(nBlockYOff), SEEK_SET));
1255 :
1256 0 : CPL_IGNORE_RET_VAL(
1257 0 : VSIFReadL(pabyRecordHeader, 1, nRecordDataStart, fp));
1258 :
1259 0 : GUInt16 nScanlineNumber = GetUInt16(pabyRecordHeader);
1260 :
1261 0 : TimeCode timeCode;
1262 0 : FetchTimeCode(&timeCode, pabyRecordHeader, nullptr);
1263 :
1264 : /* Clock drift delta */
1265 0 : i16 = GetInt16(pabyRecordHeader + 6);
1266 : /* Scanline bit field */
1267 0 : n16 = GetInt16(pabyRecordHeader + 12);
1268 :
1269 0 : CPL_IGNORE_RET_VAL(
1270 0 : VSIFPrintfL(fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,%d,", nScanlineNumber,
1271 0 : nBlockYOff, (int)timeCode.GetYear(),
1272 0 : (int)timeCode.GetDay(), (int)timeCode.GetMillisecond(),
1273 0 : i16, (n16 >> 15) & 1, (n16 >> 14) & 1, (n16)&3));
1274 :
1275 0 : n32 = GetUInt32(pabyRecordHeader + 24);
1276 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1277 : fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,",
1278 0 : (n32 >> 31) & 1, (n32 >> 30) & 1, (n32 >> 29) & 1, (n32 >> 28) & 1,
1279 0 : (n32 >> 27) & 1, (n32 >> 26) & 1, (n32 >> 25) & 1, (n32 >> 24) & 1,
1280 0 : (n32 >> 23) & 1, (n32 >> 22) & 1, (n32 >> 21) & 1, (n32 >> 20) & 1,
1281 0 : (n32 >> 8) & 1, (n32 >> 6) & 3, (n32 >> 4) & 3, (n32 >> 2) & 3,
1282 0 : (n32 >> 1) & 1, (n32 >> 0) & 1));
1283 :
1284 0 : n32 = GetUInt32(pabyRecordHeader + 28);
1285 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1286 0 : fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,", (n32 >> 23) & 1,
1287 0 : (n32 >> 22) & 1, (n32 >> 21) & 1, (n32 >> 20) & 1, (n32 >> 15) & 1,
1288 0 : (n32 >> 14) & 1, (n32 >> 13) & 1, (n32 >> 12) & 1, (n32 >> 11) & 1,
1289 0 : (n32 >> 7) & 1, (n32 >> 6) & 1, (n32 >> 5) & 1, (n32 >> 4) & 1));
1290 :
1291 0 : for (i = 0; i < 3; i++)
1292 : {
1293 0 : n16 = GetUInt16(pabyRecordHeader + 32 + 2 * i);
1294 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d,%d,%d,%d,%d,%d,",
1295 0 : (n16 >> 7) & 1, (n16 >> 6) & 1,
1296 0 : (n16 >> 5) & 1, (n16 >> 4) & 1,
1297 0 : (n16 >> 2) & 1, (n16 >> 1) & 1));
1298 : }
1299 :
1300 : /* Bit errors */
1301 0 : n16 = GetUInt16(pabyRecordHeader + 38);
1302 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d,", n16));
1303 :
1304 0 : int nOffset = 48;
1305 0 : for (i = 0; i < 3; i++)
1306 : {
1307 0 : for (j = 0; j < 3; j++)
1308 : {
1309 0 : i32 = GetInt32(pabyRecordHeader + nOffset);
1310 0 : nOffset += 4;
1311 0 : CPL_IGNORE_RET_VAL(
1312 0 : VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 7.0)));
1313 0 : i32 = GetInt32(pabyRecordHeader + nOffset);
1314 0 : nOffset += 4;
1315 0 : CPL_IGNORE_RET_VAL(
1316 0 : VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 6.0)));
1317 0 : i32 = GetInt32(pabyRecordHeader + nOffset);
1318 0 : nOffset += 4;
1319 0 : CPL_IGNORE_RET_VAL(
1320 0 : VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 7.0)));
1321 0 : i32 = GetInt32(pabyRecordHeader + nOffset);
1322 0 : nOffset += 4;
1323 0 : CPL_IGNORE_RET_VAL(
1324 0 : VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 6.0)));
1325 0 : i32 = GetInt32(pabyRecordHeader + nOffset);
1326 0 : nOffset += 4;
1327 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d,", i32));
1328 : }
1329 : }
1330 0 : for (i = 0; i < 18; i++)
1331 : {
1332 0 : i32 = GetInt32(pabyRecordHeader + nOffset);
1333 0 : nOffset += 4;
1334 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 6.0)));
1335 : }
1336 :
1337 0 : n32 = GetUInt32(pabyRecordHeader + 312);
1338 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1339 0 : fpCSV, "%d,%d,%d,%d,%d,", (n32 >> 16) & 1, (n32 >> 12) & 15,
1340 0 : (n32 >> 8) & 15, (n32 >> 4) & 15, (n32 >> 0) & 15));
1341 :
1342 0 : n32 = GetUInt32(pabyRecordHeader + 316);
1343 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d,", n32));
1344 :
1345 0 : for (i = 0; i < 3; i++)
1346 : {
1347 0 : i16 = GetUInt16(pabyRecordHeader + 320 + 2 * i);
1348 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%f,", i16 / pow(10.0, 3.0)));
1349 : }
1350 :
1351 0 : n16 = GetUInt16(pabyRecordHeader + 326);
1352 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%f", n16 / pow(10.0, 1.0)));
1353 :
1354 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));
1355 : }
1356 :
1357 0 : CPLFree(pabyRecordHeader);
1358 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpCSV));
1359 : }
1360 :
1361 : /************************************************************************/
1362 : /* EBCDICToASCII */
1363 : /************************************************************************/
1364 :
1365 : constexpr GByte EBCDICToASCII[] = {
1366 : 0x00, 0x01, 0x02, 0x03, 0x9C, 0x09, 0x86, 0x7F, 0x97, 0x8D, 0x8E, 0x0B,
1367 : 0x0C, 0x0D, 0x0E, 0x0F, 0x10, 0x11, 0x12, 0x13, 0x9D, 0x85, 0x08, 0x87,
1368 : 0x18, 0x19, 0x92, 0x8F, 0x1C, 0x1D, 0x1E, 0x1F, 0x80, 0x81, 0x82, 0x83,
1369 : 0x84, 0x0A, 0x17, 0x1B, 0x88, 0x89, 0x8A, 0x8B, 0x8C, 0x05, 0x06, 0x07,
1370 : 0x90, 0x91, 0x16, 0x93, 0x94, 0x95, 0x96, 0x04, 0x98, 0x99, 0x9A, 0x9B,
1371 : 0x14, 0x15, 0x9E, 0x1A, 0x20, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1372 : 0x00, 0x00, 0xA2, 0x2E, 0x3C, 0x28, 0x2B, 0x7C, 0x26, 0x00, 0x00, 0x00,
1373 : 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x21, 0x24, 0x2A, 0x29, 0x3B, 0xAC,
1374 : 0x2D, 0x2F, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xA6, 0x2C,
1375 : 0x25, 0x5F, 0x3E, 0x3F, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1376 : 0x00, 0x60, 0x3A, 0x23, 0x40, 0x27, 0x3D, 0x22, 0x00, 0x61, 0x62, 0x63,
1377 : 0x64, 0x65, 0x66, 0x67, 0x68, 0x69, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1378 : 0x00, 0x6A, 0x6B, 0x6C, 0x6D, 0x6E, 0x6F, 0x70, 0x71, 0x72, 0x00, 0x00,
1379 : 0x00, 0x00, 0x00, 0x00, 0x00, 0x7E, 0x73, 0x74, 0x75, 0x76, 0x77, 0x78,
1380 : 0x79, 0x7A, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1381 : 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1382 : 0x7B, 0x41, 0x42, 0x43, 0x44, 0x45, 0x46, 0x47, 0x48, 0x49, 0x00, 0x00,
1383 : 0x00, 0x00, 0x00, 0x00, 0x7D, 0x4A, 0x4B, 0x4C, 0x4D, 0x4E, 0x4F, 0x50,
1384 : 0x51, 0x52, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x5C, 0x00, 0x53, 0x54,
1385 : 0x55, 0x56, 0x57, 0x58, 0x59, 0x5A, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1386 : 0x30, 0x31, 0x32, 0x33, 0x34, 0x35, 0x36, 0x37, 0x38, 0x39, 0x00, 0x00,
1387 : 0x00, 0x00, 0x00, 0x9F,
1388 : };
1389 :
1390 : /************************************************************************/
1391 : /* ProcessDatasetHeader() */
1392 : /************************************************************************/
1393 :
1394 3 : CPLErr L1BDataset::ProcessDatasetHeader(const char *pszFilename)
1395 : {
1396 : char szDatasetName[L1B_DATASET_NAME_SIZE + 1];
1397 :
1398 3 : if (eL1BFormat == L1B_NOAA9)
1399 : {
1400 : GByte abyTBMHeader[L1B_NOAA9_HEADER_SIZE];
1401 :
1402 4 : if (VSIFSeekL(fp, 0, SEEK_SET) < 0 ||
1403 2 : VSIFReadL(abyTBMHeader, 1, L1B_NOAA9_HEADER_SIZE, fp) <
1404 : L1B_NOAA9_HEADER_SIZE)
1405 : {
1406 0 : CPLDebug("L1B", "Can't read NOAA-9/14 TBM header.");
1407 0 : return CE_Failure;
1408 : }
1409 :
1410 : // If dataset name in EBCDIC, decode it in ASCII
1411 2 : if (*(abyTBMHeader + 8 + 25) == 'K' &&
1412 0 : *(abyTBMHeader + 8 + 30) == 'K' &&
1413 0 : *(abyTBMHeader + 8 + 33) == 'K' &&
1414 0 : *(abyTBMHeader + 8 + 40) == 'K' &&
1415 0 : *(abyTBMHeader + 8 + 46) == 'K' &&
1416 0 : *(abyTBMHeader + 8 + 52) == 'K' && *(abyTBMHeader + 8 + 61) == 'K')
1417 : {
1418 0 : for (int i = 0; i < L1B_DATASET_NAME_SIZE; i++)
1419 0 : abyTBMHeader[L1B_NOAA9_HDR_NAME_OFF + i] =
1420 0 : EBCDICToASCII[abyTBMHeader[L1B_NOAA9_HDR_NAME_OFF + i]];
1421 : }
1422 :
1423 : // Fetch dataset name. NOAA-9/14 datasets contain the names in TBM
1424 : // header only, so read it there.
1425 2 : memcpy(szDatasetName, abyTBMHeader + L1B_NOAA9_HDR_NAME_OFF,
1426 : L1B_DATASET_NAME_SIZE);
1427 2 : szDatasetName[L1B_DATASET_NAME_SIZE] = '\0';
1428 :
1429 : // Deal with a few NOAA <= 9 datasets with no dataset name in TBM header
1430 2 : if (memcmp(szDatasetName,
1431 : "\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0"
1432 : "\0\0\0\0\0\0\0\0\0\0\0\0\0",
1433 0 : L1B_DATASET_NAME_SIZE) == 0 &&
1434 0 : strlen(pszFilename) == L1B_DATASET_NAME_SIZE)
1435 : {
1436 0 : strcpy(szDatasetName, pszFilename);
1437 : }
1438 :
1439 : // Determine processing center where the dataset was created
1440 2 : if (STARTS_WITH_CI(szDatasetName, "CMS"))
1441 0 : eProcCenter = CMS;
1442 2 : else if (STARTS_WITH_CI(szDatasetName, "DSS"))
1443 0 : eProcCenter = DSS;
1444 2 : else if (STARTS_WITH_CI(szDatasetName, "NSS"))
1445 2 : eProcCenter = NSS;
1446 0 : else if (STARTS_WITH_CI(szDatasetName, "UKM"))
1447 0 : eProcCenter = UKM;
1448 : else
1449 0 : eProcCenter = UNKNOWN_CENTER;
1450 :
1451 : // Determine number of bands
1452 : int i;
1453 42 : for (i = 0; i < L1B_NOAA9_HDR_CHAN_SIZE; i++)
1454 : {
1455 40 : if (abyTBMHeader[L1B_NOAA9_HDR_CHAN_OFF + i] == 1 ||
1456 38 : abyTBMHeader[L1B_NOAA9_HDR_CHAN_OFF + i] == 'Y')
1457 : {
1458 2 : nBands++;
1459 2 : iChannelsMask |= (1 << i);
1460 : }
1461 : }
1462 2 : if (nBands == 0 || nBands > 5)
1463 : {
1464 0 : nBands = 5;
1465 0 : iChannelsMask = 0x1F;
1466 : }
1467 :
1468 : // Determine data format (10-bit packed or 8/16-bit unpacked)
1469 2 : if (STARTS_WITH_CI((const char *)abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF,
1470 : "10"))
1471 0 : iDataFormat = PACKED10BIT;
1472 2 : else if (STARTS_WITH_CI(
1473 : (const char *)abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF, "16"))
1474 0 : iDataFormat = UNPACKED16BIT;
1475 2 : else if (STARTS_WITH_CI(
1476 : (const char *)abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF, "08"))
1477 2 : iDataFormat = UNPACKED8BIT;
1478 0 : else if (STARTS_WITH_CI((const char *)abyTBMHeader +
1479 : L1B_NOAA9_HDR_WORD_OFF,
1480 0 : " ") ||
1481 0 : abyTBMHeader[L1B_NOAA9_HDR_WORD_OFF] == '\0')
1482 : /* Empty string can be found in the following samples :
1483 : http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/franh.1b (10
1484 : bit) http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/frang.1b (10
1485 : bit) http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/calfilel.1b
1486 : (16 bit)
1487 : http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/rapnzg.1b (16
1488 : bit)
1489 : ftp://ftp.sat.dundee.ac.uk/misc/testdata/noaa12/hrptnoaa1b.dat
1490 : (10 bit)
1491 : */
1492 0 : bGuessDataFormat = TRUE;
1493 : else
1494 : {
1495 : #ifdef DEBUG
1496 0 : CPLDebug("L1B", "Unknown data format \"%.2s\".",
1497 : abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF);
1498 : #endif
1499 0 : return CE_Failure;
1500 : }
1501 :
1502 : // Now read the dataset header record
1503 : GByte abyRecHeader[L1B_NOAA9_HDR_REC_SIZE];
1504 4 : if (VSIFSeekL(fp, L1B_NOAA9_HEADER_SIZE, SEEK_SET) < 0 ||
1505 2 : VSIFReadL(abyRecHeader, 1, L1B_NOAA9_HDR_REC_SIZE, fp) <
1506 : L1B_NOAA9_HDR_REC_SIZE)
1507 : {
1508 0 : CPLDebug("L1B", "Can't read NOAA-9/14 record header.");
1509 0 : return CE_Failure;
1510 : }
1511 :
1512 : // Determine the spacecraft name
1513 2 : switch (abyRecHeader[L1B_NOAA9_HDR_REC_ID_OFF])
1514 : {
1515 0 : case 4:
1516 0 : eSpacecraftID = NOAA7;
1517 0 : break;
1518 0 : case 6:
1519 0 : eSpacecraftID = NOAA8;
1520 0 : break;
1521 0 : case 7:
1522 0 : eSpacecraftID = NOAA9;
1523 0 : break;
1524 0 : case 8:
1525 0 : eSpacecraftID = NOAA10;
1526 0 : break;
1527 0 : case 1:
1528 : {
1529 : /* We could also use the time code to determine TIROS-N */
1530 0 : if (strlen(pszFilename) == L1B_DATASET_NAME_SIZE &&
1531 0 : STARTS_WITH(pszFilename + 8, ".TN."))
1532 0 : eSpacecraftID = TIROSN;
1533 : else
1534 0 : eSpacecraftID = NOAA11;
1535 0 : break;
1536 : }
1537 2 : case 5:
1538 2 : eSpacecraftID = NOAA12;
1539 2 : break;
1540 0 : case 2:
1541 : {
1542 : /* We could also use the time code to determine NOAA6 */
1543 0 : if (strlen(pszFilename) == L1B_DATASET_NAME_SIZE &&
1544 0 : STARTS_WITH(pszFilename + 8, ".NA."))
1545 0 : eSpacecraftID = NOAA6;
1546 : else
1547 0 : eSpacecraftID = NOAA13;
1548 0 : break;
1549 : }
1550 0 : case 3:
1551 0 : eSpacecraftID = NOAA14;
1552 0 : break;
1553 0 : default:
1554 0 : CPLError(CE_Warning, CPLE_AppDefined,
1555 : "Unknown spacecraft ID \"%d\".",
1556 0 : abyRecHeader[L1B_NOAA9_HDR_REC_ID_OFF]);
1557 :
1558 0 : eSpacecraftID = NOAA9_UNKNOWN;
1559 0 : break;
1560 : }
1561 :
1562 : // Determine the product data type
1563 2 : int iWord = abyRecHeader[L1B_NOAA9_HDR_REC_PROD_OFF] >> 4;
1564 2 : switch (iWord)
1565 : {
1566 0 : case 1:
1567 0 : eProductType = LAC;
1568 0 : break;
1569 2 : case 2:
1570 2 : eProductType = GAC;
1571 2 : break;
1572 0 : case 3:
1573 0 : eProductType = HRPT;
1574 0 : break;
1575 0 : default:
1576 : #ifdef DEBUG
1577 0 : CPLDebug("L1B", "Unknown product type \"%d\".", iWord);
1578 : #endif
1579 0 : return CE_Failure;
1580 : }
1581 :
1582 : // Determine receiving station name
1583 2 : iWord = (abyRecHeader[L1B_NOAA9_HDR_REC_DSTAT_OFF] & 0x60) >> 5;
1584 2 : switch (iWord)
1585 : {
1586 0 : case 1:
1587 0 : eSource = GC;
1588 0 : break;
1589 2 : case 2:
1590 2 : eSource = WI;
1591 2 : break;
1592 0 : case 3:
1593 0 : eSource = SO;
1594 0 : break;
1595 0 : default:
1596 0 : eSource = UNKNOWN_STATION;
1597 0 : break;
1598 : }
1599 : }
1600 :
1601 1 : else if (eL1BFormat == L1B_NOAA15 || eL1BFormat == L1B_NOAA15_NOHDR)
1602 : {
1603 1 : if (eL1BFormat == L1B_NOAA15)
1604 : {
1605 : GByte abyARSHeader[L1B_NOAA15_HEADER_SIZE];
1606 :
1607 0 : if (VSIFSeekL(fp, 0, SEEK_SET) < 0 ||
1608 0 : VSIFReadL(abyARSHeader, 1, L1B_NOAA15_HEADER_SIZE, fp) <
1609 : L1B_NOAA15_HEADER_SIZE)
1610 : {
1611 0 : CPLDebug("L1B", "Can't read NOAA-15 ARS header.");
1612 0 : return CE_Failure;
1613 : }
1614 :
1615 : // Determine number of bands
1616 : int i;
1617 0 : for (i = 0; i < L1B_NOAA15_HDR_CHAN_SIZE; i++)
1618 : {
1619 0 : if (abyARSHeader[L1B_NOAA15_HDR_CHAN_OFF + i] == 1 ||
1620 0 : abyARSHeader[L1B_NOAA15_HDR_CHAN_OFF + i] == 'Y')
1621 : {
1622 0 : nBands++;
1623 0 : iChannelsMask |= (1 << i);
1624 : }
1625 : }
1626 0 : if (nBands == 0 || nBands > 5)
1627 : {
1628 0 : nBands = 5;
1629 0 : iChannelsMask = 0x1F;
1630 : }
1631 :
1632 : // Determine data format (10-bit packed or 8/16-bit unpacked)
1633 0 : if (STARTS_WITH_CI(
1634 : (const char *)abyARSHeader + L1B_NOAA15_HDR_WORD_OFF, "10"))
1635 0 : iDataFormat = PACKED10BIT;
1636 0 : else if (STARTS_WITH_CI((const char *)abyARSHeader +
1637 : L1B_NOAA15_HDR_WORD_OFF,
1638 : "16"))
1639 0 : iDataFormat = UNPACKED16BIT;
1640 0 : else if (STARTS_WITH_CI((const char *)abyARSHeader +
1641 : L1B_NOAA15_HDR_WORD_OFF,
1642 : "08"))
1643 0 : iDataFormat = UNPACKED8BIT;
1644 : else
1645 : {
1646 : #ifdef DEBUG
1647 0 : CPLDebug("L1B", "Unknown data format \"%.2s\".",
1648 : abyARSHeader + L1B_NOAA9_HDR_WORD_OFF);
1649 : #endif
1650 0 : return CE_Failure;
1651 : }
1652 : }
1653 : else
1654 : {
1655 1 : nBands = 5;
1656 1 : iChannelsMask = 0x1F;
1657 1 : iDataFormat = PACKED10BIT;
1658 : }
1659 :
1660 : // Now read the dataset header record
1661 : GByte abyRecHeader[L1B_NOAA15_HDR_REC_SIZE];
1662 1 : if (VSIFSeekL(fp,
1663 1 : (eL1BFormat == L1B_NOAA15) ? L1B_NOAA15_HEADER_SIZE : 0,
1664 2 : SEEK_SET) < 0 ||
1665 1 : VSIFReadL(abyRecHeader, 1, L1B_NOAA15_HDR_REC_SIZE, fp) <
1666 : L1B_NOAA15_HDR_REC_SIZE)
1667 : {
1668 0 : CPLDebug("L1B", "Can't read NOAA-9/14 record header.");
1669 0 : return CE_Failure;
1670 : }
1671 :
1672 : // Fetch dataset name
1673 1 : memcpy(szDatasetName, abyRecHeader + L1B_NOAA15_HDR_REC_NAME_OFF,
1674 : L1B_DATASET_NAME_SIZE);
1675 1 : szDatasetName[L1B_DATASET_NAME_SIZE] = '\0';
1676 :
1677 : // Determine processing center where the dataset was created
1678 1 : if (STARTS_WITH_CI((const char *)abyRecHeader +
1679 : L1B_NOAA15_HDR_REC_SITE_OFF,
1680 : "CMS"))
1681 0 : eProcCenter = CMS;
1682 1 : else if (STARTS_WITH_CI((const char *)abyRecHeader +
1683 : L1B_NOAA15_HDR_REC_SITE_OFF,
1684 : "DSS"))
1685 0 : eProcCenter = DSS;
1686 1 : else if (STARTS_WITH_CI((const char *)abyRecHeader +
1687 : L1B_NOAA15_HDR_REC_SITE_OFF,
1688 : "NSS"))
1689 0 : eProcCenter = NSS;
1690 1 : else if (STARTS_WITH_CI((const char *)abyRecHeader +
1691 : L1B_NOAA15_HDR_REC_SITE_OFF,
1692 : "UKM"))
1693 0 : eProcCenter = UKM;
1694 : else
1695 1 : eProcCenter = UNKNOWN_CENTER;
1696 :
1697 : int nFormatVersionYear, nFormatVersionDayOfYear, nHeaderRecCount;
1698 :
1699 : /* Some products from NOAA-18 and NOAA-19 coming from 'ess' processing
1700 : * station */
1701 : /* have little-endian ordering. Try to detect it with some consistency
1702 : * checks */
1703 1 : int i = 0;
1704 1 : do
1705 : {
1706 2 : nFormatVersionYear = GetUInt16(
1707 : abyRecHeader + L1B_NOAA15_HDR_REC_FORMAT_VERSION_YEAR_OFF);
1708 2 : nFormatVersionDayOfYear = GetUInt16(
1709 : abyRecHeader + L1B_NOAA15_HDR_REC_FORMAT_VERSION_DAY_OFF);
1710 2 : nHeaderRecCount =
1711 2 : GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_HDR_REC_COUNT_OFF);
1712 2 : if (i == 2)
1713 0 : break;
1714 2 : if (!(nFormatVersionYear >= 1980 && nFormatVersionYear <= 2100) &&
1715 1 : !(nFormatVersionDayOfYear <= 366) && !(nHeaderRecCount == 1))
1716 : {
1717 1 : if (i == 0)
1718 1 : CPLDebug("L1B", "Trying little-endian ordering");
1719 : else
1720 0 : CPLDebug("L1B", "Not completely convincing... Returning to "
1721 : "big-endian order");
1722 1 : bByteSwap = !bByteSwap;
1723 : }
1724 : else
1725 : break;
1726 1 : i++;
1727 1 : } while (i <= 2);
1728 1 : nRecordSizeFromHeader =
1729 1 : GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_LOGICAL_REC_LENGTH_OFF);
1730 : int nFormatVersion =
1731 1 : GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_FORMAT_VERSION_OFF);
1732 1 : CPLDebug("L1B", "NOAA Level 1b Format Version Number = %d",
1733 : nFormatVersion);
1734 1 : CPLDebug("L1B", "Level 1b Format Version Year = %d",
1735 : nFormatVersionYear);
1736 1 : CPLDebug("L1B", "Level 1b Format Version Day of Year = %d",
1737 : nFormatVersionDayOfYear);
1738 1 : CPLDebug("L1B",
1739 : "Logical Record Length of source Level 1b data set prior to "
1740 : "processing = %d",
1741 : nRecordSizeFromHeader);
1742 : int nBlockSize =
1743 1 : GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_BLOCK_SIZE_OFF);
1744 1 : CPLDebug(
1745 : "L1B",
1746 : "Block Size of source Level 1b data set prior to processing = %d",
1747 : nBlockSize);
1748 1 : CPLDebug("L1B", "Count of Header Records in this Data Set = %d",
1749 : nHeaderRecCount);
1750 :
1751 : int nDataRecordCount =
1752 1 : GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_DATA_RECORD_COUNT_OFF);
1753 1 : CPLDebug("L1B", "Count of Data Records = %d", nDataRecordCount);
1754 :
1755 1 : int nCalibratedScanlineCount = GetUInt16(
1756 1 : abyRecHeader + L1B_NOAA15_HDR_REC_CALIBRATED_SCANLINE_COUNT_OFF);
1757 1 : CPLDebug("L1B", "Count of Calibrated, Earth Located Scan Lines = %d",
1758 : nCalibratedScanlineCount);
1759 :
1760 1 : int nMissingScanlineCount = GetUInt16(
1761 1 : abyRecHeader + L1B_NOAA15_HDR_REC_MISSING_SCANLINE_COUNT_OFF);
1762 1 : CPLDebug("L1B", "Count of Missing Scan Lines = %d",
1763 : nMissingScanlineCount);
1764 1 : if (nMissingScanlineCount != 0)
1765 1 : bExposeMaskBand = TRUE;
1766 :
1767 : char szEllipsoid[8 + 1];
1768 1 : memcpy(szEllipsoid, abyRecHeader + L1B_NOAA15_HDR_REC_ELLIPSOID_OFF, 8);
1769 1 : szEllipsoid[8] = '\0';
1770 1 : CPLDebug("L1B", "Reference Ellipsoid Model ID = '%s'", szEllipsoid);
1771 1 : if (EQUAL(szEllipsoid, "WGS-84 "))
1772 : {
1773 0 : m_oGCPSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
1774 : }
1775 1 : else if (EQUAL(szEllipsoid, " GRS 80"))
1776 : {
1777 1 : m_oGCPSRS.importFromWkt(
1778 : "GEOGCS[\"GRS 1980(IUGG, "
1779 : "1980)\",DATUM[\"unknown\",SPHEROID[\"GRS80\",6378137,298."
1780 : "257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],"
1781 : "UNIT[\"degree\",0.0174532925199433]]");
1782 : }
1783 :
1784 : // Determine the spacecraft name
1785 : // See
1786 : // http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c8/sec83132-2.htm
1787 1 : int iWord = GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_ID_OFF);
1788 1 : switch (iWord)
1789 : {
1790 0 : case 2:
1791 0 : eSpacecraftID = NOAA16;
1792 0 : break;
1793 0 : case 4:
1794 0 : eSpacecraftID = NOAA15;
1795 0 : break;
1796 0 : case 6:
1797 0 : eSpacecraftID = NOAA17;
1798 0 : break;
1799 0 : case 7:
1800 0 : eSpacecraftID = NOAA18;
1801 0 : break;
1802 1 : case 8:
1803 1 : eSpacecraftID = NOAA19;
1804 1 : break;
1805 0 : case 11:
1806 0 : eSpacecraftID = METOP1;
1807 0 : break;
1808 0 : case 12:
1809 0 : eSpacecraftID = METOP2;
1810 0 : break;
1811 : // METOP3 is not documented yet
1812 0 : case 13:
1813 0 : eSpacecraftID = METOP3;
1814 0 : break;
1815 0 : case 14:
1816 0 : eSpacecraftID = METOP3;
1817 0 : break;
1818 0 : default:
1819 : #ifdef DEBUG
1820 0 : CPLDebug("L1B", "Unknown spacecraft ID \"%d\".", iWord);
1821 : #endif
1822 0 : return CE_Failure;
1823 : }
1824 :
1825 : // Determine the product data type
1826 1 : iWord = GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_PROD_OFF);
1827 1 : switch (iWord)
1828 : {
1829 0 : case 1:
1830 0 : eProductType = LAC;
1831 0 : break;
1832 0 : case 2:
1833 0 : eProductType = GAC;
1834 0 : break;
1835 1 : case 3:
1836 1 : eProductType = HRPT;
1837 1 : break;
1838 0 : case 4: // XXX: documentation specifies the code '4'
1839 : case 13: // for FRAC but real datasets contain '13 here.'
1840 0 : eProductType = FRAC;
1841 0 : break;
1842 0 : default:
1843 : #ifdef DEBUG
1844 0 : CPLDebug("L1B", "Unknown product type \"%d\".", iWord);
1845 : #endif
1846 0 : return CE_Failure;
1847 : }
1848 :
1849 : // Fetch instrument status. Helps to determine whether we have
1850 : // 3A or 3B channel in the dataset.
1851 1 : iInstrumentStatus =
1852 1 : GetUInt32(abyRecHeader + L1B_NOAA15_HDR_REC_STAT_OFF);
1853 :
1854 : // Determine receiving station name
1855 1 : iWord = GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_SRC_OFF);
1856 1 : switch (iWord)
1857 : {
1858 0 : case 1:
1859 0 : eSource = GC;
1860 0 : break;
1861 0 : case 2:
1862 0 : eSource = WI;
1863 0 : break;
1864 0 : case 3:
1865 0 : eSource = SO;
1866 0 : break;
1867 0 : case 4:
1868 0 : eSource = SV;
1869 0 : break;
1870 0 : case 5:
1871 0 : eSource = MO;
1872 0 : break;
1873 1 : default:
1874 1 : eSource = UNKNOWN_STATION;
1875 1 : break;
1876 1 : }
1877 : }
1878 : else
1879 0 : return CE_Failure;
1880 :
1881 : /* -------------------------------------------------------------------- */
1882 : /* Set fetched information as metadata records */
1883 : /* -------------------------------------------------------------------- */
1884 :
1885 3 : SetMetadataItem("DATASET_NAME", szDatasetName);
1886 :
1887 3 : const char *pszText = nullptr;
1888 3 : switch (eSpacecraftID)
1889 : {
1890 0 : case TIROSN:
1891 0 : pszText = "TIROS-N";
1892 0 : break;
1893 0 : case NOAA6:
1894 0 : pszText = "NOAA-6(A)";
1895 0 : break;
1896 0 : case NOAAB:
1897 0 : pszText = "NOAA-B";
1898 0 : break;
1899 0 : case NOAA7:
1900 0 : pszText = "NOAA-7(C)";
1901 0 : break;
1902 0 : case NOAA8:
1903 0 : pszText = "NOAA-8(E)";
1904 0 : break;
1905 0 : case NOAA9_UNKNOWN:
1906 0 : pszText = "UNKNOWN";
1907 0 : break;
1908 0 : case NOAA9:
1909 0 : pszText = "NOAA-9(F)";
1910 0 : break;
1911 0 : case NOAA10:
1912 0 : pszText = "NOAA-10(G)";
1913 0 : break;
1914 0 : case NOAA11:
1915 0 : pszText = "NOAA-11(H)";
1916 0 : break;
1917 2 : case NOAA12:
1918 2 : pszText = "NOAA-12(D)";
1919 2 : break;
1920 0 : case NOAA13:
1921 0 : pszText = "NOAA-13(I)";
1922 0 : break;
1923 0 : case NOAA14:
1924 0 : pszText = "NOAA-14(J)";
1925 0 : break;
1926 0 : case NOAA15:
1927 0 : pszText = "NOAA-15(K)";
1928 0 : break;
1929 0 : case NOAA16:
1930 0 : pszText = "NOAA-16(L)";
1931 0 : break;
1932 0 : case NOAA17:
1933 0 : pszText = "NOAA-17(M)";
1934 0 : break;
1935 0 : case NOAA18:
1936 0 : pszText = "NOAA-18(N)";
1937 0 : break;
1938 1 : case NOAA19:
1939 1 : pszText = "NOAA-19(N')";
1940 1 : break;
1941 0 : case METOP2:
1942 0 : pszText = "METOP-A(2)";
1943 0 : break;
1944 0 : case METOP1:
1945 0 : pszText = "METOP-B(1)";
1946 0 : break;
1947 0 : case METOP3:
1948 0 : pszText = "METOP-C(3)";
1949 0 : break;
1950 0 : default:
1951 0 : pszText = "Unknown";
1952 0 : break;
1953 : }
1954 3 : SetMetadataItem("SATELLITE", pszText);
1955 :
1956 3 : switch (eProductType)
1957 : {
1958 0 : case LAC:
1959 0 : pszText = "AVHRR LAC";
1960 0 : break;
1961 1 : case HRPT:
1962 1 : pszText = "AVHRR HRPT";
1963 1 : break;
1964 2 : case GAC:
1965 2 : pszText = "AVHRR GAC";
1966 2 : break;
1967 0 : case FRAC:
1968 0 : pszText = "AVHRR FRAC";
1969 0 : break;
1970 0 : default:
1971 0 : pszText = "Unknown";
1972 0 : break;
1973 : }
1974 3 : SetMetadataItem("DATA_TYPE", pszText);
1975 :
1976 : // Get revolution number as string, we don't need this value for processing
1977 : char szRevolution[6];
1978 3 : memcpy(szRevolution, szDatasetName + 32, 5);
1979 3 : szRevolution[5] = '\0';
1980 3 : SetMetadataItem("REVOLUTION", szRevolution);
1981 :
1982 3 : switch (eSource)
1983 : {
1984 0 : case DU:
1985 0 : pszText = "Dundee, Scotland, UK";
1986 0 : break;
1987 0 : case GC:
1988 0 : pszText = "Fairbanks, Alaska, USA (formerly Gilmore Creek)";
1989 0 : break;
1990 0 : case HO:
1991 0 : pszText = "Honolulu, Hawaii, USA";
1992 0 : break;
1993 0 : case MO:
1994 0 : pszText = "Monterey, California, USA";
1995 0 : break;
1996 0 : case WE:
1997 0 : pszText = "Western Europe CDA, Lannion, France";
1998 0 : break;
1999 0 : case SO:
2000 0 : pszText = "SOCC (Satellite Operations Control Center), Suitland, "
2001 : "Maryland, USA";
2002 0 : break;
2003 2 : case WI:
2004 2 : pszText = "Wallops Island, Virginia, USA";
2005 2 : break;
2006 1 : default:
2007 1 : pszText = "Unknown receiving station";
2008 1 : break;
2009 : }
2010 3 : SetMetadataItem("SOURCE", pszText);
2011 :
2012 3 : switch (eProcCenter)
2013 : {
2014 0 : case CMS:
2015 0 : pszText = "Centre de Meteorologie Spatiale - Lannion, France";
2016 0 : break;
2017 0 : case DSS:
2018 0 : pszText =
2019 : "Dundee Satellite Receiving Station - Dundee, Scotland, UK";
2020 0 : break;
2021 2 : case NSS:
2022 2 : pszText = "NOAA/NESDIS - Suitland, Maryland, USA";
2023 2 : break;
2024 0 : case UKM:
2025 0 : pszText =
2026 : "United Kingdom Meteorological Office - Bracknell, England, UK";
2027 0 : break;
2028 1 : default:
2029 1 : pszText = "Unknown processing center";
2030 1 : break;
2031 : }
2032 3 : SetMetadataItem("PROCESSING_CENTER", pszText);
2033 :
2034 3 : return CE_None;
2035 : }
2036 :
2037 : /************************************************************************/
2038 : /* ComputeFileOffsets() */
2039 : /************************************************************************/
2040 :
2041 3 : int L1BDataset::ComputeFileOffsets()
2042 : {
2043 3 : CPLDebug("L1B", "Data format = %s",
2044 3 : (iDataFormat == PACKED10BIT) ? "Packed 10 bit"
2045 3 : : (iDataFormat == UNPACKED16BIT) ? "Unpacked 16 bit"
2046 : : "Unpacked 8 bit");
2047 :
2048 3 : switch (eProductType)
2049 : {
2050 1 : case HRPT:
2051 : case LAC:
2052 : case FRAC:
2053 1 : nRasterXSize = 2048;
2054 1 : nBufferSize = 20484;
2055 1 : iGCPStart =
2056 : 25 -
2057 : 1; /* http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c2/sec2-4.htm
2058 : */
2059 1 : iGCPStep = 40;
2060 1 : nGCPsPerLine = 51;
2061 1 : if (eL1BFormat == L1B_NOAA9)
2062 : {
2063 0 : if (iDataFormat == PACKED10BIT)
2064 : {
2065 0 : nRecordSize = 14800;
2066 0 : nRecordDataEnd = 14104;
2067 : }
2068 0 : else if (iDataFormat == UNPACKED16BIT)
2069 : {
2070 0 : switch (nBands)
2071 : {
2072 0 : case 1:
2073 0 : nRecordSize = 4544;
2074 0 : nRecordDataEnd = 4544;
2075 0 : break;
2076 0 : case 2:
2077 0 : nRecordSize = 8640;
2078 0 : nRecordDataEnd = 8640;
2079 0 : break;
2080 0 : case 3:
2081 0 : nRecordSize = 12736;
2082 0 : nRecordDataEnd = 12736;
2083 0 : break;
2084 0 : case 4:
2085 0 : nRecordSize = 16832;
2086 0 : nRecordDataEnd = 16832;
2087 0 : break;
2088 0 : case 5:
2089 0 : nRecordSize = 20928;
2090 0 : nRecordDataEnd = 20928;
2091 0 : break;
2092 : }
2093 : }
2094 : else // UNPACKED8BIT
2095 : {
2096 0 : switch (nBands)
2097 : {
2098 0 : case 1:
2099 0 : nRecordSize = 2496;
2100 0 : nRecordDataEnd = 2496;
2101 0 : break;
2102 0 : case 2:
2103 0 : nRecordSize = 4544;
2104 0 : nRecordDataEnd = 4544;
2105 0 : break;
2106 0 : case 3:
2107 0 : nRecordSize = 6592;
2108 0 : nRecordDataEnd = 6592;
2109 0 : break;
2110 0 : case 4:
2111 0 : nRecordSize = 8640;
2112 0 : nRecordDataEnd = 8640;
2113 0 : break;
2114 0 : case 5:
2115 0 : nRecordSize = 10688;
2116 0 : nRecordDataEnd = 10688;
2117 0 : break;
2118 : }
2119 : }
2120 0 : nDataStartOffset = nRecordSize + L1B_NOAA9_HEADER_SIZE;
2121 0 : nRecordDataStart = 448;
2122 0 : iGCPCodeOffset = 52;
2123 0 : iGCPOffset = 104;
2124 : }
2125 :
2126 1 : else if (eL1BFormat == L1B_NOAA15 || eL1BFormat == L1B_NOAA15_NOHDR)
2127 : {
2128 1 : if (iDataFormat == PACKED10BIT)
2129 : {
2130 0 : nRecordSize = 15872;
2131 0 : nRecordDataEnd = 14920;
2132 0 : iCLAVRStart = 14984;
2133 : }
2134 1 : else if (iDataFormat == UNPACKED16BIT)
2135 : { /* Table 8.3.1.3.3.1-3 */
2136 1 : switch (nBands)
2137 : {
2138 0 : case 1:
2139 0 : nRecordSize = 6144;
2140 0 : nRecordDataEnd = 5360;
2141 0 : iCLAVRStart =
2142 : 5368 + 56; /* guessed but not verified */
2143 0 : break;
2144 0 : case 2:
2145 0 : nRecordSize = 10240;
2146 0 : nRecordDataEnd = 9456;
2147 0 : iCLAVRStart =
2148 : 9464 + 56; /* guessed but not verified */
2149 0 : break;
2150 0 : case 3:
2151 0 : nRecordSize = 14336;
2152 0 : nRecordDataEnd = 13552;
2153 0 : iCLAVRStart =
2154 : 13560 + 56; /* guessed but not verified */
2155 0 : break;
2156 0 : case 4:
2157 0 : nRecordSize = 18432;
2158 0 : nRecordDataEnd = 17648;
2159 0 : iCLAVRStart =
2160 : 17656 + 56; /* guessed but not verified */
2161 0 : break;
2162 1 : case 5:
2163 1 : nRecordSize = 22528;
2164 1 : nRecordDataEnd = 21744;
2165 1 : iCLAVRStart = 21752 + 56;
2166 1 : break;
2167 : }
2168 : }
2169 : else // UNPACKED8BIT
2170 : { /* Table 8.3.1.3.3.1-2 */
2171 0 : switch (nBands)
2172 : {
2173 0 : case 1:
2174 0 : nRecordSize = 4096;
2175 0 : nRecordDataEnd = 3312;
2176 0 : iCLAVRStart =
2177 : 3320 + 56; /* guessed but not verified */
2178 0 : break;
2179 0 : case 2:
2180 0 : nRecordSize = 6144;
2181 0 : nRecordDataEnd = 5360;
2182 0 : iCLAVRStart =
2183 : 5368 + 56; /* guessed but not verified */
2184 0 : break;
2185 0 : case 3:
2186 0 : nRecordSize = 8192;
2187 0 : nRecordDataEnd = 7408;
2188 0 : iCLAVRStart =
2189 : 7416 + 56; /* guessed but not verified */
2190 0 : break;
2191 0 : case 4:
2192 0 : nRecordSize = 10240;
2193 0 : nRecordDataEnd = 9456;
2194 0 : iCLAVRStart =
2195 : 9464 + 56; /* guessed but not verified */
2196 0 : break;
2197 0 : case 5:
2198 0 : nRecordSize = 12288;
2199 0 : nRecordDataEnd = 11504;
2200 0 : iCLAVRStart =
2201 : 11512 + 56; /* guessed but not verified */
2202 0 : break;
2203 : }
2204 : }
2205 2 : nDataStartOffset = (eL1BFormat == L1B_NOAA15_NOHDR)
2206 1 : ? nRecordDataEnd
2207 0 : : nRecordSize + L1B_NOAA15_HEADER_SIZE;
2208 1 : nRecordDataStart = 1264;
2209 1 : iGCPCodeOffset = 0; // XXX: not exist for NOAA15?
2210 1 : iGCPOffset = 640;
2211 : }
2212 : else
2213 0 : return 0;
2214 1 : break;
2215 :
2216 2 : case GAC:
2217 2 : nRasterXSize = 409;
2218 2 : nBufferSize = 4092;
2219 2 : iGCPStart = 5 - 1; // FIXME: depends of scan direction
2220 2 : iGCPStep = 8;
2221 2 : nGCPsPerLine = 51;
2222 2 : if (eL1BFormat == L1B_NOAA9)
2223 : {
2224 2 : if (iDataFormat == PACKED10BIT)
2225 : {
2226 0 : nRecordSize = 3220;
2227 0 : nRecordDataEnd = 3176;
2228 : }
2229 2 : else if (iDataFormat == UNPACKED16BIT)
2230 0 : switch (nBands)
2231 : {
2232 0 : case 1:
2233 0 : nRecordSize = 1268;
2234 0 : nRecordDataEnd = 1266;
2235 0 : break;
2236 0 : case 2:
2237 0 : nRecordSize = 2084;
2238 0 : nRecordDataEnd = 2084;
2239 0 : break;
2240 0 : case 3:
2241 0 : nRecordSize = 2904;
2242 0 : nRecordDataEnd = 2902;
2243 0 : break;
2244 0 : case 4:
2245 0 : nRecordSize = 3720;
2246 0 : nRecordDataEnd = 3720;
2247 0 : break;
2248 0 : case 5:
2249 0 : nRecordSize = 4540;
2250 0 : nRecordDataEnd = 4538;
2251 0 : break;
2252 : }
2253 : else // UNPACKED8BIT
2254 : {
2255 2 : switch (nBands)
2256 : {
2257 2 : case 1:
2258 2 : nRecordSize = 860;
2259 2 : nRecordDataEnd = 858;
2260 2 : break;
2261 0 : case 2:
2262 0 : nRecordSize = 1268;
2263 0 : nRecordDataEnd = 1266;
2264 0 : break;
2265 0 : case 3:
2266 0 : nRecordSize = 1676;
2267 0 : nRecordDataEnd = 1676;
2268 0 : break;
2269 0 : case 4:
2270 0 : nRecordSize = 2084;
2271 0 : nRecordDataEnd = 2084;
2272 0 : break;
2273 0 : case 5:
2274 0 : nRecordSize = 2496;
2275 0 : nRecordDataEnd = 2494;
2276 0 : break;
2277 : }
2278 : }
2279 2 : nDataStartOffset = nRecordSize * 2 + L1B_NOAA9_HEADER_SIZE;
2280 2 : nRecordDataStart = 448;
2281 2 : iGCPCodeOffset = 52;
2282 2 : iGCPOffset = 104;
2283 : }
2284 :
2285 0 : else if (eL1BFormat == L1B_NOAA15 || eL1BFormat == L1B_NOAA15_NOHDR)
2286 : {
2287 0 : if (iDataFormat == PACKED10BIT)
2288 : {
2289 0 : nRecordSize = 4608;
2290 0 : nRecordDataEnd = 3992;
2291 0 : iCLAVRStart = 4056;
2292 : }
2293 0 : else if (iDataFormat == UNPACKED16BIT)
2294 : { /* Table 8.3.1.4.3.1-3 */
2295 0 : switch (nBands)
2296 : {
2297 0 : case 1:
2298 0 : nRecordSize = 2360;
2299 0 : nRecordDataEnd = 2082;
2300 0 : iCLAVRStart =
2301 : 2088 + 56; /* guessed but not verified */
2302 0 : break;
2303 0 : case 2:
2304 0 : nRecordSize = 3176;
2305 0 : nRecordDataEnd = 2900;
2306 0 : iCLAVRStart =
2307 : 2904 + 56; /* guessed but not verified */
2308 0 : break;
2309 0 : case 3:
2310 0 : nRecordSize = 3992;
2311 0 : nRecordDataEnd = 3718;
2312 0 : iCLAVRStart =
2313 : 3720 + 56; /* guessed but not verified */
2314 0 : break;
2315 0 : case 4:
2316 0 : nRecordSize = 4816;
2317 0 : nRecordDataEnd = 4536;
2318 0 : iCLAVRStart =
2319 : 4544 + 56; /* guessed but not verified */
2320 0 : break;
2321 0 : case 5:
2322 0 : nRecordSize = 5632;
2323 0 : nRecordDataEnd = 5354;
2324 0 : iCLAVRStart = 5360 + 56;
2325 0 : break;
2326 : }
2327 : }
2328 : else // UNPACKED8BIT
2329 : { /* Table 8.3.1.4.3.1-2 but record length is wrong in the table
2330 : ! */
2331 0 : switch (nBands)
2332 : {
2333 0 : case 1:
2334 0 : nRecordSize = 1952;
2335 0 : nRecordDataEnd = 1673;
2336 0 : iCLAVRStart =
2337 : 1680 + 56; /* guessed but not verified */
2338 0 : break;
2339 0 : case 2:
2340 0 : nRecordSize = 2360;
2341 0 : nRecordDataEnd = 2082;
2342 0 : iCLAVRStart =
2343 : 2088 + 56; /* guessed but not verified */
2344 0 : break;
2345 0 : case 3:
2346 0 : nRecordSize = 2768;
2347 0 : nRecordDataEnd = 2491;
2348 0 : iCLAVRStart =
2349 : 2496 + 56; /* guessed but not verified */
2350 0 : break;
2351 0 : case 4:
2352 0 : nRecordSize = 3176;
2353 0 : nRecordDataEnd = 2900;
2354 0 : iCLAVRStart =
2355 : 2904 + 56; /* guessed but not verified */
2356 0 : break;
2357 0 : case 5:
2358 0 : nRecordSize = 3584;
2359 0 : nRecordDataEnd = 3309;
2360 0 : iCLAVRStart =
2361 : 3312 + 56; /* guessed but not verified */
2362 0 : break;
2363 : }
2364 : }
2365 0 : nDataStartOffset = (eL1BFormat == L1B_NOAA15_NOHDR)
2366 0 : ? nRecordDataEnd
2367 0 : : nRecordSize + L1B_NOAA15_HEADER_SIZE;
2368 0 : nRecordDataStart = 1264;
2369 0 : iGCPCodeOffset = 0; // XXX: not exist for NOAA15?
2370 0 : iGCPOffset = 640;
2371 : }
2372 : else
2373 0 : return 0;
2374 2 : break;
2375 0 : default:
2376 0 : return 0;
2377 : }
2378 :
2379 3 : return 1;
2380 : }
2381 :
2382 : /************************************************************************/
2383 : /* L1BGeolocDataset */
2384 : /************************************************************************/
2385 :
2386 : class L1BGeolocDataset final : public GDALDataset
2387 : {
2388 : friend class L1BGeolocRasterBand;
2389 :
2390 : L1BDataset *poL1BDS;
2391 : int bInterpolGeolocationDS;
2392 :
2393 : public:
2394 : L1BGeolocDataset(L1BDataset *poMainDS, int bInterpolGeolocationDS);
2395 : virtual ~L1BGeolocDataset();
2396 :
2397 : static GDALDataset *CreateGeolocationDS(L1BDataset *poL1BDS,
2398 : int bInterpolGeolocationDS);
2399 : };
2400 :
2401 : /************************************************************************/
2402 : /* L1BGeolocRasterBand */
2403 : /************************************************************************/
2404 :
2405 : class L1BGeolocRasterBand final : public GDALRasterBand
2406 : {
2407 : public:
2408 : L1BGeolocRasterBand(L1BGeolocDataset *poDS, int nBand);
2409 :
2410 : virtual CPLErr IReadBlock(int, int, void *) override;
2411 : virtual double GetNoDataValue(int *pbSuccess = nullptr) override;
2412 : };
2413 :
2414 : /************************************************************************/
2415 : /* L1BGeolocDataset() */
2416 : /************************************************************************/
2417 :
2418 0 : L1BGeolocDataset::L1BGeolocDataset(L1BDataset *poL1BDSIn,
2419 0 : int bInterpolGeolocationDSIn)
2420 0 : : poL1BDS(poL1BDSIn), bInterpolGeolocationDS(bInterpolGeolocationDSIn)
2421 : {
2422 0 : if (bInterpolGeolocationDS)
2423 0 : nRasterXSize = poL1BDS->nRasterXSize;
2424 : else
2425 0 : nRasterXSize = poL1BDS->nGCPsPerLine;
2426 0 : nRasterYSize = poL1BDS->nRasterYSize;
2427 0 : }
2428 :
2429 : /************************************************************************/
2430 : /* ~L1BGeolocDataset() */
2431 : /************************************************************************/
2432 :
2433 0 : L1BGeolocDataset::~L1BGeolocDataset()
2434 : {
2435 0 : delete poL1BDS;
2436 0 : }
2437 :
2438 : /************************************************************************/
2439 : /* L1BGeolocRasterBand() */
2440 : /************************************************************************/
2441 :
2442 0 : L1BGeolocRasterBand::L1BGeolocRasterBand(L1BGeolocDataset *poDSIn, int nBandIn)
2443 : {
2444 0 : poDS = poDSIn;
2445 0 : nBand = nBandIn;
2446 0 : nRasterXSize = poDSIn->nRasterXSize;
2447 0 : nRasterYSize = poDSIn->nRasterYSize;
2448 0 : eDataType = GDT_Float64;
2449 0 : nBlockXSize = nRasterXSize;
2450 0 : nBlockYSize = 1;
2451 0 : if (nBand == 1)
2452 0 : SetDescription("GEOLOC X");
2453 : else
2454 0 : SetDescription("GEOLOC Y");
2455 0 : }
2456 :
2457 : /************************************************************************/
2458 : /* LagrangeInterpol() */
2459 : /************************************************************************/
2460 :
2461 : /* ----------------------------------------------------------------------------
2462 : * Perform a Lagrangian interpolation through the given x,y coordinates
2463 : * and return the interpolated y value for the given x value.
2464 : * The array size and thus the polynomial order is defined by numpt.
2465 : * Input: x[] and y[] are of size numpt,
2466 : * x0 is the x value for which we calculate the corresponding y
2467 : * Returns: y value calculated for given x0.
2468 : */
2469 0 : static double LagrangeInterpol(const double x[], const double y[], double x0,
2470 : int numpt)
2471 : {
2472 : int i, j;
2473 : double L;
2474 0 : double y0 = 0;
2475 :
2476 0 : for (i = 0; i < numpt; i++)
2477 : {
2478 0 : L = 1.0;
2479 0 : for (j = 0; j < numpt; j++)
2480 : {
2481 0 : if (i == j)
2482 0 : continue;
2483 0 : L = L * (x0 - x[j]) / (x[i] - x[j]);
2484 : }
2485 0 : y0 = y0 + L * y[i];
2486 : }
2487 0 : return y0;
2488 : }
2489 :
2490 : /************************************************************************/
2491 : /* L1BInterpol() */
2492 : /************************************************************************/
2493 :
2494 : /* ----------------------------------------------------------------------------
2495 : * Interpolate an array of size numPoints where the only values set on input are
2496 : * at knownFirst, and intervals of knownStep thereafter.
2497 : * On return all the rest from 0..numPoints-1 will be filled in.
2498 : * Uses the LagrangeInterpol() function to do the interpolation; 5-point for the
2499 : * beginning and end of the array and 4-point for the rest.
2500 : * To use this function for NOAA level 1B data extract the 51 latitude values
2501 : * into their appropriate places in the vals array then call L1BInterpol to
2502 : * calculate the rest of the values. Do similarly for longitudes, solar zenith
2503 : * angles, and any others which are present in the file.
2504 : * Reference:
2505 : * http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c2/sec2-4.htm
2506 : */
2507 :
2508 : #define MIDDLE_INTERP_ORDER 4
2509 : #define END_INTERP_ORDER 5 /* Ensure this is an odd number, 5 is suitable.*/
2510 :
2511 : /* Convert number of known point to its index in full array */
2512 : #define IDX(N) ((N)*knownStep + knownFirst)
2513 :
2514 0 : static void L1BInterpol(
2515 : double vals[], int numKnown, /* Number of known points (typically 51) */
2516 : int knownFirst, /* Index in full array of first known point (24) */
2517 : int knownStep, /* Interval to next and subsequent known points (40) */
2518 : int numPoints /* Number of points in whole array (2048) */)
2519 : {
2520 : int i, j;
2521 : double x[END_INTERP_ORDER];
2522 : double y[END_INTERP_ORDER];
2523 :
2524 : /* First extrapolate first 24 points */
2525 0 : for (i = 0; i < END_INTERP_ORDER; i++)
2526 : {
2527 0 : x[i] = IDX(i);
2528 0 : y[i] = vals[IDX(i)];
2529 : }
2530 0 : for (i = 0; i < knownFirst; i++)
2531 : {
2532 0 : vals[i] = LagrangeInterpol(x, y, i, END_INTERP_ORDER);
2533 : }
2534 :
2535 : /* Next extrapolate last 23 points */
2536 0 : for (i = 0; i < END_INTERP_ORDER; i++)
2537 : {
2538 0 : x[i] = IDX(numKnown - END_INTERP_ORDER + i);
2539 0 : y[i] = vals[IDX(numKnown - END_INTERP_ORDER + i)];
2540 : }
2541 0 : for (i = IDX(numKnown - 1); i < numPoints; i++)
2542 : {
2543 0 : vals[i] = LagrangeInterpol(x, y, i, END_INTERP_ORDER);
2544 : }
2545 :
2546 : /* Interpolate all intermediate points using two before and two after */
2547 0 : for (i = knownFirst; i < IDX(numKnown - 1); i++)
2548 : {
2549 : double x2[MIDDLE_INTERP_ORDER];
2550 : double y2[MIDDLE_INTERP_ORDER];
2551 : int startpt;
2552 :
2553 : /* Find a suitable set of two known points before and two after */
2554 0 : startpt = (i / knownStep) - MIDDLE_INTERP_ORDER / 2;
2555 0 : if (startpt < 0)
2556 0 : startpt = 0;
2557 0 : if (startpt + MIDDLE_INTERP_ORDER - 1 >= numKnown)
2558 0 : startpt = numKnown - MIDDLE_INTERP_ORDER;
2559 0 : for (j = 0; j < MIDDLE_INTERP_ORDER; j++)
2560 : {
2561 0 : x2[j] = IDX(startpt + j);
2562 0 : y2[j] = vals[IDX(startpt + j)];
2563 : }
2564 0 : vals[i] = LagrangeInterpol(x2, y2, i, MIDDLE_INTERP_ORDER);
2565 : }
2566 0 : }
2567 :
2568 : /************************************************************************/
2569 : /* IReadBlock() */
2570 : /************************************************************************/
2571 :
2572 0 : CPLErr L1BGeolocRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
2573 : int nBlockYOff, void *pData)
2574 : {
2575 0 : L1BGeolocDataset *poGDS = cpl::down_cast<L1BGeolocDataset *>(poDS);
2576 0 : L1BDataset *poL1BDS = poGDS->poL1BDS;
2577 : GDAL_GCP *pasGCPList =
2578 0 : (GDAL_GCP *)CPLCalloc(poL1BDS->nGCPsPerLine, sizeof(GDAL_GCP));
2579 0 : GDALInitGCPs(poL1BDS->nGCPsPerLine, pasGCPList);
2580 :
2581 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);
2582 :
2583 : /* -------------------------------------------------------------------- */
2584 : /* Seek to data. */
2585 : /* -------------------------------------------------------------------- */
2586 0 : CPL_IGNORE_RET_VAL(
2587 0 : VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));
2588 :
2589 0 : CPL_IGNORE_RET_VAL(
2590 0 : VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordDataStart, poL1BDS->fp));
2591 :
2592 : /* Fetch the GCPs for the row */
2593 0 : int nGotGCPs = poL1BDS->FetchGCPs(pasGCPList, pabyRecordHeader, nBlockYOff);
2594 0 : double *padfData = (double *)pData;
2595 : int i;
2596 0 : if (poGDS->bInterpolGeolocationDS)
2597 : {
2598 : /* Fill the known position */
2599 0 : for (i = 0; i < nGotGCPs; i++)
2600 : {
2601 0 : double dfVal =
2602 0 : (nBand == 1) ? pasGCPList[i].dfGCPX : pasGCPList[i].dfGCPY;
2603 0 : padfData[poL1BDS->iGCPStart + i * poL1BDS->iGCPStep] = dfVal;
2604 : }
2605 :
2606 0 : if (nGotGCPs == poL1BDS->nGCPsPerLine)
2607 : {
2608 : /* And do Lagangian interpolation to fill the holes */
2609 0 : L1BInterpol(padfData, poL1BDS->nGCPsPerLine, poL1BDS->iGCPStart,
2610 : poL1BDS->iGCPStep, nRasterXSize);
2611 : }
2612 : else
2613 : {
2614 0 : int iFirstNonValid = 0;
2615 0 : if (nGotGCPs > 5)
2616 0 : iFirstNonValid = poL1BDS->iGCPStart +
2617 0 : nGotGCPs * poL1BDS->iGCPStep +
2618 0 : poL1BDS->iGCPStep / 2;
2619 0 : for (i = iFirstNonValid; i < nRasterXSize; i++)
2620 : {
2621 0 : padfData[i] = GetNoDataValue(nullptr);
2622 : }
2623 0 : if (iFirstNonValid > 0)
2624 : {
2625 0 : L1BInterpol(padfData, poL1BDS->nGCPsPerLine, poL1BDS->iGCPStart,
2626 : poL1BDS->iGCPStep, iFirstNonValid);
2627 : }
2628 : }
2629 : }
2630 : else
2631 : {
2632 0 : for (i = 0; i < nGotGCPs; i++)
2633 : {
2634 0 : padfData[i] =
2635 0 : (nBand == 1) ? pasGCPList[i].dfGCPX : pasGCPList[i].dfGCPY;
2636 : }
2637 0 : for (i = nGotGCPs; i < nRasterXSize; i++)
2638 0 : padfData[i] = GetNoDataValue(nullptr);
2639 : }
2640 :
2641 0 : if (poL1BDS->eLocationIndicator == ASCEND)
2642 : {
2643 0 : for (i = 0; i < nRasterXSize / 2; i++)
2644 : {
2645 0 : double dfTmp = padfData[i];
2646 0 : padfData[i] = padfData[nRasterXSize - 1 - i];
2647 0 : padfData[nRasterXSize - 1 - i] = dfTmp;
2648 : }
2649 : }
2650 :
2651 0 : CPLFree(pabyRecordHeader);
2652 0 : GDALDeinitGCPs(poL1BDS->nGCPsPerLine, pasGCPList);
2653 0 : CPLFree(pasGCPList);
2654 :
2655 0 : return CE_None;
2656 : }
2657 :
2658 : /************************************************************************/
2659 : /* GetNoDataValue() */
2660 : /************************************************************************/
2661 :
2662 0 : double L1BGeolocRasterBand::GetNoDataValue(int *pbSuccess)
2663 : {
2664 0 : if (pbSuccess)
2665 0 : *pbSuccess = TRUE;
2666 0 : return -200.0;
2667 : }
2668 :
2669 : /************************************************************************/
2670 : /* CreateGeolocationDS() */
2671 : /************************************************************************/
2672 :
2673 0 : GDALDataset *L1BGeolocDataset::CreateGeolocationDS(L1BDataset *poL1BDS,
2674 : int bInterpolGeolocationDS)
2675 : {
2676 : L1BGeolocDataset *poGeolocDS =
2677 0 : new L1BGeolocDataset(poL1BDS, bInterpolGeolocationDS);
2678 0 : for (int i = 1; i <= 2; i++)
2679 : {
2680 0 : poGeolocDS->SetBand(i, new L1BGeolocRasterBand(poGeolocDS, i));
2681 : }
2682 0 : return poGeolocDS;
2683 : }
2684 :
2685 : /************************************************************************/
2686 : /* L1BSolarZenithAnglesDataset */
2687 : /************************************************************************/
2688 :
2689 : class L1BSolarZenithAnglesDataset final : public GDALDataset
2690 : {
2691 : friend class L1BSolarZenithAnglesRasterBand;
2692 :
2693 : L1BDataset *poL1BDS;
2694 :
2695 : public:
2696 : explicit L1BSolarZenithAnglesDataset(L1BDataset *poMainDS);
2697 : virtual ~L1BSolarZenithAnglesDataset();
2698 :
2699 : static GDALDataset *CreateSolarZenithAnglesDS(L1BDataset *poL1BDS);
2700 : };
2701 :
2702 : /************************************************************************/
2703 : /* L1BSolarZenithAnglesRasterBand */
2704 : /************************************************************************/
2705 :
2706 : class L1BSolarZenithAnglesRasterBand final : public GDALRasterBand
2707 : {
2708 : public:
2709 : L1BSolarZenithAnglesRasterBand(L1BSolarZenithAnglesDataset *poDS,
2710 : int nBand);
2711 :
2712 : virtual CPLErr IReadBlock(int, int, void *) override;
2713 : virtual double GetNoDataValue(int *pbSuccess = nullptr) override;
2714 : };
2715 :
2716 : /************************************************************************/
2717 : /* L1BSolarZenithAnglesDataset() */
2718 : /************************************************************************/
2719 :
2720 0 : L1BSolarZenithAnglesDataset::L1BSolarZenithAnglesDataset(L1BDataset *poL1BDSIn)
2721 : {
2722 0 : poL1BDS = poL1BDSIn;
2723 0 : nRasterXSize = 51;
2724 0 : nRasterYSize = poL1BDSIn->nRasterYSize;
2725 0 : }
2726 :
2727 : /************************************************************************/
2728 : /* ~L1BSolarZenithAnglesDataset() */
2729 : /************************************************************************/
2730 :
2731 0 : L1BSolarZenithAnglesDataset::~L1BSolarZenithAnglesDataset()
2732 : {
2733 0 : delete poL1BDS;
2734 0 : }
2735 :
2736 : /************************************************************************/
2737 : /* L1BSolarZenithAnglesRasterBand() */
2738 : /************************************************************************/
2739 :
2740 0 : L1BSolarZenithAnglesRasterBand::L1BSolarZenithAnglesRasterBand(
2741 0 : L1BSolarZenithAnglesDataset *poDSIn, int nBandIn)
2742 : {
2743 0 : poDS = poDSIn;
2744 0 : nBand = nBandIn;
2745 0 : nRasterXSize = poDSIn->nRasterXSize;
2746 0 : nRasterYSize = poDSIn->nRasterYSize;
2747 0 : eDataType = GDT_Float32;
2748 0 : nBlockXSize = nRasterXSize;
2749 0 : nBlockYSize = 1;
2750 0 : }
2751 :
2752 : /************************************************************************/
2753 : /* IReadBlock() */
2754 : /************************************************************************/
2755 :
2756 0 : CPLErr L1BSolarZenithAnglesRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
2757 : int nBlockYOff, void *pData)
2758 : {
2759 : L1BSolarZenithAnglesDataset *poGDS =
2760 0 : cpl::down_cast<L1BSolarZenithAnglesDataset *>(poDS);
2761 0 : L1BDataset *poL1BDS = poGDS->poL1BDS;
2762 : int i;
2763 :
2764 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);
2765 :
2766 : /* -------------------------------------------------------------------- */
2767 : /* Seek to data. */
2768 : /* -------------------------------------------------------------------- */
2769 0 : CPL_IGNORE_RET_VAL(
2770 0 : VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));
2771 :
2772 0 : CPL_IGNORE_RET_VAL(
2773 0 : VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordSize, poL1BDS->fp));
2774 :
2775 : const int nValidValues =
2776 0 : std::min(nRasterXSize,
2777 0 : static_cast<int>(pabyRecordHeader[poL1BDS->iGCPCodeOffset]));
2778 0 : float *pafData = (float *)pData;
2779 :
2780 0 : int bHasFractional = (poL1BDS->nRecordDataEnd + 20 <= poL1BDS->nRecordSize);
2781 :
2782 : #ifdef notdef
2783 : if (bHasFractional)
2784 : {
2785 : for (i = 0; i < 20; i++)
2786 : {
2787 : GByte val = pabyRecordHeader[poL1BDS->nRecordDataEnd + i];
2788 : for (int j = 0; j < 8; j++)
2789 : fprintf(stderr, "%c", /*ok*/
2790 : ((val >> (7 - j)) & 1) ? '1' : '0');
2791 : fprintf(stderr, " "); /*ok*/
2792 : }
2793 : fprintf(stderr, "\n"); /*ok*/
2794 : }
2795 : #endif
2796 :
2797 0 : for (i = 0; i < nValidValues; i++)
2798 : {
2799 0 : pafData[i] = pabyRecordHeader[poL1BDS->iGCPCodeOffset + 1 + i] / 2.0f;
2800 :
2801 0 : if (bHasFractional)
2802 : {
2803 : /* Cf
2804 : * http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/l/app-l.htm#notl-2
2805 : */
2806 : /* This is not very clear on if bits must be counted from MSB or LSB
2807 : */
2808 : /* but when testing on n12gac10bit.l1b, it appears that the 3 bits
2809 : * for i=0 are the 3 MSB bits */
2810 0 : int nAddBitStart = i * 3;
2811 : int nFractional;
2812 :
2813 : #if 1
2814 0 : if ((nAddBitStart % 8) + 3 <= 8)
2815 : {
2816 0 : nFractional = (pabyRecordHeader[poL1BDS->nRecordDataEnd +
2817 0 : (nAddBitStart / 8)] >>
2818 0 : (8 - ((nAddBitStart % 8) + 3))) &
2819 : 0x7;
2820 : }
2821 : else
2822 : {
2823 0 : nFractional = (((pabyRecordHeader[poL1BDS->nRecordDataEnd +
2824 0 : (nAddBitStart / 8)]
2825 0 : << 8) |
2826 0 : pabyRecordHeader[poL1BDS->nRecordDataEnd +
2827 0 : (nAddBitStart / 8) + 1]) >>
2828 0 : (16 - ((nAddBitStart % 8) + 3))) &
2829 : 0x7;
2830 : }
2831 : #else
2832 : nFractional = (pabyRecordHeader[poL1BDS->nRecordDataEnd +
2833 : (nAddBitStart / 8)] >>
2834 : (nAddBitStart % 8)) &
2835 : 0x7;
2836 : if ((nAddBitStart % 8) + 3 > 8)
2837 : nFractional |= (pabyRecordHeader[poL1BDS->nRecordDataEnd +
2838 : (nAddBitStart / 8) + 1] &
2839 : ((1 << (((nAddBitStart % 8) + 3 - 8))) - 1))
2840 : << (3 - ((((nAddBitStart % 8) + 3 - 8))));
2841 : * /
2842 : #endif
2843 0 : if (nFractional > 4)
2844 : {
2845 0 : CPLDebug("L1B",
2846 : "For nBlockYOff=%d, i=%d, wrong fractional value : %d",
2847 : nBlockYOff, i, nFractional);
2848 : }
2849 :
2850 0 : pafData[i] += nFractional / 10.0f;
2851 : }
2852 : }
2853 :
2854 0 : for (; i < nRasterXSize; i++)
2855 0 : pafData[i] = static_cast<float>(GetNoDataValue(nullptr));
2856 :
2857 0 : if (poL1BDS->eLocationIndicator == ASCEND)
2858 : {
2859 0 : for (i = 0; i < nRasterXSize / 2; i++)
2860 : {
2861 0 : float fTmp = pafData[i];
2862 0 : pafData[i] = pafData[nRasterXSize - 1 - i];
2863 0 : pafData[nRasterXSize - 1 - i] = fTmp;
2864 : }
2865 : }
2866 :
2867 0 : CPLFree(pabyRecordHeader);
2868 :
2869 0 : return CE_None;
2870 : }
2871 :
2872 : /************************************************************************/
2873 : /* GetNoDataValue() */
2874 : /************************************************************************/
2875 :
2876 0 : double L1BSolarZenithAnglesRasterBand::GetNoDataValue(int *pbSuccess)
2877 : {
2878 0 : if (pbSuccess)
2879 0 : *pbSuccess = TRUE;
2880 0 : return -200.0;
2881 : }
2882 :
2883 : /************************************************************************/
2884 : /* CreateSolarZenithAnglesDS() */
2885 : /************************************************************************/
2886 :
2887 : GDALDataset *
2888 0 : L1BSolarZenithAnglesDataset::CreateSolarZenithAnglesDS(L1BDataset *poL1BDS)
2889 : {
2890 : L1BSolarZenithAnglesDataset *poGeolocDS =
2891 0 : new L1BSolarZenithAnglesDataset(poL1BDS);
2892 0 : for (int i = 1; i <= 1; i++)
2893 : {
2894 0 : poGeolocDS->SetBand(i,
2895 0 : new L1BSolarZenithAnglesRasterBand(poGeolocDS, i));
2896 : }
2897 0 : return poGeolocDS;
2898 : }
2899 :
2900 : /************************************************************************/
2901 : /* L1BNOAA15AnglesDataset */
2902 : /************************************************************************/
2903 :
2904 : class L1BNOAA15AnglesDataset final : public GDALDataset
2905 : {
2906 : friend class L1BNOAA15AnglesRasterBand;
2907 :
2908 : L1BDataset *poL1BDS;
2909 :
2910 : public:
2911 : explicit L1BNOAA15AnglesDataset(L1BDataset *poMainDS);
2912 : virtual ~L1BNOAA15AnglesDataset();
2913 :
2914 : static GDALDataset *CreateAnglesDS(L1BDataset *poL1BDS);
2915 : };
2916 :
2917 : /************************************************************************/
2918 : /* L1BNOAA15AnglesRasterBand */
2919 : /************************************************************************/
2920 :
2921 : class L1BNOAA15AnglesRasterBand final : public GDALRasterBand
2922 : {
2923 : public:
2924 : L1BNOAA15AnglesRasterBand(L1BNOAA15AnglesDataset *poDS, int nBand);
2925 :
2926 : virtual CPLErr IReadBlock(int, int, void *) override;
2927 : };
2928 :
2929 : /************************************************************************/
2930 : /* L1BNOAA15AnglesDataset() */
2931 : /************************************************************************/
2932 :
2933 0 : L1BNOAA15AnglesDataset::L1BNOAA15AnglesDataset(L1BDataset *poL1BDSIn)
2934 0 : : poL1BDS(poL1BDSIn)
2935 : {
2936 0 : nRasterXSize = 51;
2937 0 : nRasterYSize = poL1BDS->nRasterYSize;
2938 0 : }
2939 :
2940 : /************************************************************************/
2941 : /* ~L1BNOAA15AnglesDataset() */
2942 : /************************************************************************/
2943 :
2944 0 : L1BNOAA15AnglesDataset::~L1BNOAA15AnglesDataset()
2945 : {
2946 0 : delete poL1BDS;
2947 0 : }
2948 :
2949 : /************************************************************************/
2950 : /* L1BNOAA15AnglesRasterBand() */
2951 : /************************************************************************/
2952 :
2953 0 : L1BNOAA15AnglesRasterBand::L1BNOAA15AnglesRasterBand(
2954 0 : L1BNOAA15AnglesDataset *poDSIn, int nBandIn)
2955 : {
2956 0 : poDS = poDSIn;
2957 0 : nBand = nBandIn;
2958 0 : nRasterXSize = poDSIn->nRasterXSize;
2959 0 : nRasterYSize = poDSIn->nRasterYSize;
2960 0 : eDataType = GDT_Float32;
2961 0 : nBlockXSize = nRasterXSize;
2962 0 : nBlockYSize = 1;
2963 0 : if (nBand == 1)
2964 0 : SetDescription("Solar zenith angles");
2965 0 : else if (nBand == 2)
2966 0 : SetDescription("Satellite zenith angles");
2967 : else
2968 0 : SetDescription("Relative azimuth angles");
2969 0 : }
2970 :
2971 : /************************************************************************/
2972 : /* IReadBlock() */
2973 : /************************************************************************/
2974 :
2975 0 : CPLErr L1BNOAA15AnglesRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
2976 : int nBlockYOff, void *pData)
2977 : {
2978 : L1BNOAA15AnglesDataset *poGDS =
2979 0 : cpl::down_cast<L1BNOAA15AnglesDataset *>(poDS);
2980 0 : L1BDataset *poL1BDS = poGDS->poL1BDS;
2981 : int i;
2982 :
2983 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);
2984 :
2985 : /* -------------------------------------------------------------------- */
2986 : /* Seek to data. */
2987 : /* -------------------------------------------------------------------- */
2988 0 : CPL_IGNORE_RET_VAL(
2989 0 : VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));
2990 :
2991 0 : CPL_IGNORE_RET_VAL(
2992 0 : VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordSize, poL1BDS->fp));
2993 :
2994 0 : float *pafData = (float *)pData;
2995 :
2996 0 : for (i = 0; i < nRasterXSize; i++)
2997 : {
2998 : GInt16 i16 =
2999 0 : poL1BDS->GetInt16(pabyRecordHeader + 328 + 6 * i + 2 * (nBand - 1));
3000 0 : pafData[i] = i16 / 100.0f;
3001 : }
3002 :
3003 0 : if (poL1BDS->eLocationIndicator == ASCEND)
3004 : {
3005 0 : for (i = 0; i < nRasterXSize / 2; i++)
3006 : {
3007 0 : float fTmp = pafData[i];
3008 0 : pafData[i] = pafData[nRasterXSize - 1 - i];
3009 0 : pafData[nRasterXSize - 1 - i] = fTmp;
3010 : }
3011 : }
3012 :
3013 0 : CPLFree(pabyRecordHeader);
3014 :
3015 0 : return CE_None;
3016 : }
3017 :
3018 : /************************************************************************/
3019 : /* CreateAnglesDS() */
3020 : /************************************************************************/
3021 :
3022 0 : GDALDataset *L1BNOAA15AnglesDataset::CreateAnglesDS(L1BDataset *poL1BDS)
3023 : {
3024 0 : L1BNOAA15AnglesDataset *poGeolocDS = new L1BNOAA15AnglesDataset(poL1BDS);
3025 0 : for (int i = 1; i <= 3; i++)
3026 : {
3027 0 : poGeolocDS->SetBand(i, new L1BNOAA15AnglesRasterBand(poGeolocDS, i));
3028 : }
3029 0 : return poGeolocDS;
3030 : }
3031 :
3032 : /************************************************************************/
3033 : /* L1BCloudsDataset */
3034 : /************************************************************************/
3035 :
3036 : class L1BCloudsDataset final : public GDALDataset
3037 : {
3038 : friend class L1BCloudsRasterBand;
3039 :
3040 : L1BDataset *poL1BDS;
3041 :
3042 : public:
3043 : explicit L1BCloudsDataset(L1BDataset *poMainDS);
3044 : virtual ~L1BCloudsDataset();
3045 :
3046 : static GDALDataset *CreateCloudsDS(L1BDataset *poL1BDS);
3047 : };
3048 :
3049 : /************************************************************************/
3050 : /* L1BCloudsRasterBand */
3051 : /************************************************************************/
3052 :
3053 : class L1BCloudsRasterBand final : public GDALRasterBand
3054 : {
3055 : public:
3056 : L1BCloudsRasterBand(L1BCloudsDataset *poDS, int nBand);
3057 :
3058 : virtual CPLErr IReadBlock(int, int, void *) override;
3059 : };
3060 :
3061 : /************************************************************************/
3062 : /* L1BCloudsDataset() */
3063 : /************************************************************************/
3064 :
3065 0 : L1BCloudsDataset::L1BCloudsDataset(L1BDataset *poL1BDSIn) : poL1BDS(poL1BDSIn)
3066 : {
3067 0 : nRasterXSize = poL1BDSIn->nRasterXSize;
3068 0 : nRasterYSize = poL1BDSIn->nRasterYSize;
3069 0 : }
3070 :
3071 : /************************************************************************/
3072 : /* ~L1BCloudsDataset() */
3073 : /************************************************************************/
3074 :
3075 0 : L1BCloudsDataset::~L1BCloudsDataset()
3076 : {
3077 0 : delete poL1BDS;
3078 0 : }
3079 :
3080 : /************************************************************************/
3081 : /* L1BCloudsRasterBand() */
3082 : /************************************************************************/
3083 :
3084 0 : L1BCloudsRasterBand::L1BCloudsRasterBand(L1BCloudsDataset *poDSIn, int nBandIn)
3085 : {
3086 0 : poDS = poDSIn;
3087 0 : nBand = nBandIn;
3088 0 : nRasterXSize = poDSIn->nRasterXSize;
3089 0 : nRasterYSize = poDSIn->nRasterYSize;
3090 0 : eDataType = GDT_Byte;
3091 0 : nBlockXSize = nRasterXSize;
3092 0 : nBlockYSize = 1;
3093 0 : }
3094 :
3095 : /************************************************************************/
3096 : /* IReadBlock() */
3097 : /************************************************************************/
3098 :
3099 0 : CPLErr L1BCloudsRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
3100 : int nBlockYOff, void *pData)
3101 : {
3102 0 : L1BCloudsDataset *poGDS = cpl::down_cast<L1BCloudsDataset *>(poDS);
3103 0 : L1BDataset *poL1BDS = poGDS->poL1BDS;
3104 : int i;
3105 :
3106 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);
3107 :
3108 : /* -------------------------------------------------------------------- */
3109 : /* Seek to data. */
3110 : /* -------------------------------------------------------------------- */
3111 0 : CPL_IGNORE_RET_VAL(
3112 0 : VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));
3113 :
3114 0 : CPL_IGNORE_RET_VAL(
3115 0 : VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordSize, poL1BDS->fp));
3116 :
3117 0 : GByte *pabyData = (GByte *)pData;
3118 :
3119 0 : for (i = 0; i < nRasterXSize; i++)
3120 : {
3121 0 : pabyData[i] = ((pabyRecordHeader[poL1BDS->iCLAVRStart + (i / 4)] >>
3122 0 : (8 - ((i % 4) * 2 + 2))) &
3123 : 0x3);
3124 : }
3125 :
3126 0 : if (poL1BDS->eLocationIndicator == ASCEND)
3127 : {
3128 0 : for (i = 0; i < nRasterXSize / 2; i++)
3129 : {
3130 0 : GByte byTmp = pabyData[i];
3131 0 : pabyData[i] = pabyData[nRasterXSize - 1 - i];
3132 0 : pabyData[nRasterXSize - 1 - i] = byTmp;
3133 : }
3134 : }
3135 :
3136 0 : CPLFree(pabyRecordHeader);
3137 :
3138 0 : return CE_None;
3139 : }
3140 :
3141 : /************************************************************************/
3142 : /* CreateCloudsDS() */
3143 : /************************************************************************/
3144 :
3145 0 : GDALDataset *L1BCloudsDataset::CreateCloudsDS(L1BDataset *poL1BDS)
3146 : {
3147 0 : L1BCloudsDataset *poGeolocDS = new L1BCloudsDataset(poL1BDS);
3148 0 : for (int i = 1; i <= 1; i++)
3149 : {
3150 0 : poGeolocDS->SetBand(i, new L1BCloudsRasterBand(poGeolocDS, i));
3151 : }
3152 0 : return poGeolocDS;
3153 : }
3154 :
3155 : /************************************************************************/
3156 : /* DetectFormat() */
3157 : /************************************************************************/
3158 :
3159 58755 : L1BFileFormat L1BDataset::DetectFormat(const char *pszFilename,
3160 : const GByte *pabyHeader,
3161 : int nHeaderBytes)
3162 :
3163 : {
3164 58755 : if (pabyHeader == nullptr || nHeaderBytes < L1B_NOAA9_HEADER_SIZE)
3165 54179 : return L1B_NONE;
3166 :
3167 : // try NOAA-18 formats
3168 4576 : if (*(pabyHeader + 0) == '\0' && *(pabyHeader + 1) == '\0' &&
3169 816 : *(pabyHeader + 2) == '\0' && *(pabyHeader + 3) == '\0' &&
3170 140 : *(pabyHeader + 4) == '\0' && *(pabyHeader + 5) == '\0' &&
3171 130 : EQUALN((const char *)(pabyHeader + 22), "/N1BD/N18/", 10))
3172 0 : return L1B_NOAA15_NOHDR;
3173 :
3174 : // We will try the NOAA-15 and later formats first
3175 4576 : if (nHeaderBytes > L1B_NOAA15_HEADER_SIZE + 61 &&
3176 3536 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 25) == '.' &&
3177 12 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 30) == '.' &&
3178 0 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 33) == '.' &&
3179 0 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 40) == '.' &&
3180 0 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 46) == '.' &&
3181 0 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 52) == '.' &&
3182 0 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 61) == '.')
3183 0 : return L1B_NOAA15;
3184 :
3185 : // Next try the NOAA-9/14 formats
3186 4576 : if (*(pabyHeader + 8 + 25) == '.' && *(pabyHeader + 8 + 30) == '.' &&
3187 4 : *(pabyHeader + 8 + 33) == '.' && *(pabyHeader + 8 + 40) == '.' &&
3188 4 : *(pabyHeader + 8 + 46) == '.' && *(pabyHeader + 8 + 52) == '.' &&
3189 4 : *(pabyHeader + 8 + 61) == '.')
3190 4 : return L1B_NOAA9;
3191 :
3192 : // Next try the NOAA-9/14 formats with dataset name in EBCDIC
3193 4572 : if (*(pabyHeader + 8 + 25) == 'K' && *(pabyHeader + 8 + 30) == 'K' &&
3194 0 : *(pabyHeader + 8 + 33) == 'K' && *(pabyHeader + 8 + 40) == 'K' &&
3195 0 : *(pabyHeader + 8 + 46) == 'K' && *(pabyHeader + 8 + 52) == 'K' &&
3196 0 : *(pabyHeader + 8 + 61) == 'K')
3197 0 : return L1B_NOAA9;
3198 :
3199 : // Finally try the AAPP formats
3200 4572 : if (*(pabyHeader + 25) == '.' && *(pabyHeader + 30) == '.' &&
3201 2 : *(pabyHeader + 33) == '.' && *(pabyHeader + 40) == '.' &&
3202 2 : *(pabyHeader + 46) == '.' && *(pabyHeader + 52) == '.' &&
3203 2 : *(pabyHeader + 61) == '.')
3204 2 : return L1B_NOAA15_NOHDR;
3205 :
3206 : // A few NOAA <= 9 datasets with no dataset name in TBM header
3207 4570 : if (strlen(pszFilename) == L1B_DATASET_NAME_SIZE && pszFilename[3] == '.' &&
3208 1 : pszFilename[8] == '.' && pszFilename[11] == '.' &&
3209 0 : pszFilename[18] == '.' && pszFilename[24] == '.' &&
3210 0 : pszFilename[30] == '.' && pszFilename[39] == '.' &&
3211 0 : memcmp(pabyHeader + 30,
3212 : "\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0"
3213 : "\0\0\0\0\0\0\0\0\0\0\0",
3214 0 : L1B_DATASET_NAME_SIZE) == 0 &&
3215 0 : (pabyHeader[75] == '+' || pabyHeader[75] == '-') &&
3216 0 : (pabyHeader[78] == '+' || pabyHeader[78] == '-') &&
3217 0 : (pabyHeader[81] == '+' || pabyHeader[81] == '-') &&
3218 0 : (pabyHeader[85] == '+' || pabyHeader[85] == '-'))
3219 0 : return L1B_NOAA9;
3220 :
3221 4570 : return L1B_NONE;
3222 : }
3223 :
3224 : /************************************************************************/
3225 : /* Identify() */
3226 : /************************************************************************/
3227 :
3228 58752 : int L1BDataset::Identify(GDALOpenInfo *poOpenInfo)
3229 :
3230 : {
3231 58752 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS:"))
3232 0 : return TRUE;
3233 58752 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS_INTERPOL:"))
3234 0 : return TRUE;
3235 58752 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_SOLAR_ZENITH_ANGLES:"))
3236 0 : return TRUE;
3237 58752 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_ANGLES:"))
3238 0 : return TRUE;
3239 58752 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_CLOUDS:"))
3240 0 : return TRUE;
3241 :
3242 58753 : if (DetectFormat(CPLGetFilename(poOpenInfo->pszFilename),
3243 58752 : poOpenInfo->pabyHeader,
3244 58752 : poOpenInfo->nHeaderBytes) == L1B_NONE)
3245 58749 : return FALSE;
3246 :
3247 3 : return TRUE;
3248 : }
3249 :
3250 : /************************************************************************/
3251 : /* Open() */
3252 : /************************************************************************/
3253 :
3254 3 : GDALDataset *L1BDataset::Open(GDALOpenInfo *poOpenInfo)
3255 :
3256 : {
3257 3 : GDALDataset *poOutDS = nullptr;
3258 3 : VSILFILE *fp = nullptr;
3259 6 : CPLString osFilename = poOpenInfo->pszFilename;
3260 3 : int bAskGeolocationDS = FALSE;
3261 3 : int bInterpolGeolocationDS = FALSE;
3262 3 : int bAskSolarZenithAnglesDS = FALSE;
3263 3 : int bAskAnglesDS = FALSE;
3264 3 : int bAskCloudsDS = FALSE;
3265 : L1BFileFormat eL1BFormat;
3266 :
3267 3 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS:") ||
3268 3 : STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS_INTERPOL:") ||
3269 3 : STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_SOLAR_ZENITH_ANGLES:") ||
3270 3 : STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_ANGLES:") ||
3271 3 : STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_CLOUDS:"))
3272 : {
3273 : GByte abyHeader[1024];
3274 0 : const char *pszFilename = nullptr;
3275 0 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS_INTERPOL:"))
3276 : {
3277 0 : bAskGeolocationDS = TRUE;
3278 0 : bInterpolGeolocationDS = TRUE;
3279 0 : pszFilename = poOpenInfo->pszFilename + strlen("L1BGCPS_INTERPOL:");
3280 : }
3281 0 : else if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS:"))
3282 : {
3283 0 : bAskGeolocationDS = TRUE;
3284 0 : pszFilename = poOpenInfo->pszFilename + strlen("L1BGCPS:");
3285 : }
3286 0 : else if (STARTS_WITH_CI(poOpenInfo->pszFilename,
3287 : "L1B_SOLAR_ZENITH_ANGLES:"))
3288 : {
3289 0 : bAskSolarZenithAnglesDS = TRUE;
3290 0 : pszFilename =
3291 0 : poOpenInfo->pszFilename + strlen("L1B_SOLAR_ZENITH_ANGLES:");
3292 : }
3293 0 : else if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_ANGLES:"))
3294 : {
3295 0 : bAskAnglesDS = TRUE;
3296 0 : pszFilename = poOpenInfo->pszFilename + strlen("L1B_ANGLES:");
3297 : }
3298 : else
3299 : {
3300 0 : bAskCloudsDS = TRUE;
3301 0 : pszFilename = poOpenInfo->pszFilename + strlen("L1B_CLOUDS:");
3302 : }
3303 0 : if (pszFilename[0] == '"')
3304 0 : pszFilename++;
3305 0 : osFilename = pszFilename;
3306 0 : if (!osFilename.empty() && osFilename.back() == '"')
3307 0 : osFilename.pop_back();
3308 0 : fp = VSIFOpenL(osFilename, "rb");
3309 0 : if (!fp)
3310 : {
3311 0 : CPLError(CE_Failure, CPLE_AppDefined, "Can't open file \"%s\".",
3312 : osFilename.c_str());
3313 0 : return nullptr;
3314 : }
3315 0 : CPL_IGNORE_RET_VAL(VSIFReadL(abyHeader, 1, sizeof(abyHeader) - 1, fp));
3316 0 : abyHeader[sizeof(abyHeader) - 1] = '\0';
3317 0 : eL1BFormat = DetectFormat(CPLGetFilename(osFilename), abyHeader,
3318 0 : sizeof(abyHeader));
3319 : }
3320 : else
3321 : eL1BFormat =
3322 3 : DetectFormat(CPLGetFilename(osFilename), poOpenInfo->pabyHeader,
3323 : poOpenInfo->nHeaderBytes);
3324 :
3325 3 : if (eL1BFormat == L1B_NONE)
3326 : {
3327 0 : if (fp != nullptr)
3328 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
3329 0 : return nullptr;
3330 : }
3331 :
3332 : /* -------------------------------------------------------------------- */
3333 : /* Confirm the requested access is supported. */
3334 : /* -------------------------------------------------------------------- */
3335 3 : if (poOpenInfo->eAccess == GA_Update)
3336 : {
3337 0 : ReportUpdateNotSupportedByDriver("L1B");
3338 0 : if (fp != nullptr)
3339 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
3340 0 : return nullptr;
3341 : }
3342 :
3343 : /* -------------------------------------------------------------------- */
3344 : /* Create a corresponding GDALDataset. */
3345 : /* -------------------------------------------------------------------- */
3346 : L1BDataset *poDS;
3347 : VSIStatBufL sStat;
3348 :
3349 3 : poDS = new L1BDataset(eL1BFormat);
3350 :
3351 3 : if (fp == nullptr)
3352 3 : fp = VSIFOpenL(osFilename, "rb");
3353 3 : poDS->fp = fp;
3354 3 : if (!poDS->fp || VSIStatL(osFilename, &sStat) != 0)
3355 : {
3356 0 : CPLDebug("L1B", "Can't open file \"%s\".", osFilename.c_str());
3357 0 : goto bad;
3358 : }
3359 :
3360 : /* -------------------------------------------------------------------- */
3361 : /* Read the header. */
3362 : /* -------------------------------------------------------------------- */
3363 3 : if (poDS->ProcessDatasetHeader(CPLGetFilename(osFilename)) != CE_None)
3364 : {
3365 0 : CPLDebug("L1B", "Error reading L1B record header.");
3366 0 : goto bad;
3367 : }
3368 :
3369 3 : if (poDS->eL1BFormat == L1B_NOAA15_NOHDR &&
3370 1 : poDS->nRecordSizeFromHeader == 22016 &&
3371 1 : (sStat.st_size % 22016 /* poDS->nRecordSizeFromHeader*/) == 0)
3372 : {
3373 1 : poDS->iDataFormat = UNPACKED16BIT;
3374 1 : poDS->ComputeFileOffsets();
3375 1 : poDS->nDataStartOffset = poDS->nRecordSizeFromHeader;
3376 1 : poDS->nRecordSize = poDS->nRecordSizeFromHeader;
3377 1 : poDS->iCLAVRStart = 0;
3378 : }
3379 2 : else if (poDS->bGuessDataFormat)
3380 : {
3381 : int nTempYSize;
3382 : GUInt16 nScanlineNumber;
3383 : int j;
3384 :
3385 : /* If the data format is unspecified, try each one of the 3 known data
3386 : * formats */
3387 : /* It is considered valid when the spacing between the first 5 scanline
3388 : * numbers */
3389 : /* is a constant */
3390 :
3391 0 : for (j = 0; j < 3; j++)
3392 : {
3393 0 : poDS->iDataFormat = (L1BDataFormat)(PACKED10BIT + j);
3394 0 : if (!poDS->ComputeFileOffsets())
3395 0 : goto bad;
3396 :
3397 0 : nTempYSize = static_cast<int>(
3398 0 : (sStat.st_size - poDS->nDataStartOffset) / poDS->nRecordSize);
3399 0 : if (nTempYSize < 5)
3400 0 : continue;
3401 :
3402 0 : int nLastScanlineNumber = 0;
3403 0 : int nDiffLine = 0;
3404 : int i;
3405 0 : for (i = 0; i < 5; i++)
3406 : {
3407 0 : nScanlineNumber = 0;
3408 :
3409 0 : CPL_IGNORE_RET_VAL(VSIFSeekL(
3410 0 : poDS->fp, poDS->nDataStartOffset + i * poDS->nRecordSize,
3411 : SEEK_SET));
3412 0 : CPL_IGNORE_RET_VAL(VSIFReadL(&nScanlineNumber, 1, 2, poDS->fp));
3413 0 : nScanlineNumber = poDS->GetUInt16(&nScanlineNumber);
3414 :
3415 0 : if (i == 1)
3416 : {
3417 0 : nDiffLine = nScanlineNumber - nLastScanlineNumber;
3418 0 : if (nDiffLine == 0)
3419 0 : break;
3420 : }
3421 0 : else if (i > 1)
3422 : {
3423 0 : if (nDiffLine != nScanlineNumber - nLastScanlineNumber)
3424 0 : break;
3425 : }
3426 :
3427 0 : nLastScanlineNumber = nScanlineNumber;
3428 : }
3429 :
3430 0 : if (i == 5)
3431 : {
3432 0 : CPLDebug("L1B", "Guessed data format : %s",
3433 0 : (poDS->iDataFormat == PACKED10BIT) ? "10"
3434 0 : : (poDS->iDataFormat == UNPACKED8BIT) ? "08"
3435 : : "16");
3436 0 : break;
3437 : }
3438 : }
3439 :
3440 0 : if (j == 3)
3441 : {
3442 0 : CPLError(CE_Failure, CPLE_AppDefined,
3443 : "Could not guess data format of L1B product");
3444 0 : goto bad;
3445 : }
3446 : }
3447 : else
3448 : {
3449 2 : if (!poDS->ComputeFileOffsets())
3450 0 : goto bad;
3451 : }
3452 :
3453 3 : CPLDebug("L1B", "nRecordDataStart = %d", poDS->nRecordDataStart);
3454 3 : CPLDebug("L1B", "nRecordDataEnd = %d", poDS->nRecordDataEnd);
3455 3 : CPLDebug("L1B", "nDataStartOffset = %d", poDS->nDataStartOffset);
3456 3 : CPLDebug("L1B", "iCLAVRStart = %d", poDS->iCLAVRStart);
3457 3 : CPLDebug("L1B", "nRecordSize = %d", poDS->nRecordSize);
3458 :
3459 : // Compute number of lines dynamically, so we can read partially
3460 : // downloaded files.
3461 3 : if (poDS->nDataStartOffset > sStat.st_size)
3462 0 : goto bad;
3463 3 : poDS->nRasterYSize = static_cast<int>(
3464 3 : (sStat.st_size - poDS->nDataStartOffset) / poDS->nRecordSize);
3465 :
3466 : /* -------------------------------------------------------------------- */
3467 : /* Deal with GCPs */
3468 : /* -------------------------------------------------------------------- */
3469 3 : poDS->ProcessRecordHeaders();
3470 :
3471 3 : if (bAskGeolocationDS)
3472 : {
3473 0 : return L1BGeolocDataset::CreateGeolocationDS(poDS,
3474 0 : bInterpolGeolocationDS);
3475 : }
3476 3 : else if (bAskSolarZenithAnglesDS)
3477 : {
3478 0 : if (eL1BFormat == L1B_NOAA9)
3479 0 : return L1BSolarZenithAnglesDataset::CreateSolarZenithAnglesDS(poDS);
3480 : else
3481 : {
3482 0 : delete poDS;
3483 0 : return nullptr;
3484 : }
3485 : }
3486 3 : else if (bAskAnglesDS)
3487 : {
3488 0 : if (eL1BFormat != L1B_NOAA9)
3489 0 : return L1BNOAA15AnglesDataset::CreateAnglesDS(poDS);
3490 : else
3491 : {
3492 0 : delete poDS;
3493 0 : return nullptr;
3494 : }
3495 : }
3496 3 : else if (bAskCloudsDS)
3497 : {
3498 0 : if (poDS->iCLAVRStart > 0)
3499 0 : poOutDS = L1BCloudsDataset::CreateCloudsDS(poDS);
3500 : else
3501 : {
3502 0 : delete poDS;
3503 0 : return nullptr;
3504 : }
3505 : }
3506 : else
3507 : {
3508 3 : poOutDS = poDS;
3509 : }
3510 :
3511 : {
3512 6 : CPLString osTMP;
3513 : int bInterpol =
3514 3 : CPLTestBool(CPLGetConfigOption("L1B_INTERPOL_GCPS", "TRUE"));
3515 :
3516 3 : char *pszWKT = nullptr;
3517 3 : poDS->m_oGCPSRS.exportToWkt(&pszWKT);
3518 3 : poOutDS->SetMetadataItem("SRS", pszWKT,
3519 3 : "GEOLOCATION"); /* unused by gdalgeoloc.cpp */
3520 3 : CPLFree(pszWKT);
3521 :
3522 3 : if (bInterpol)
3523 3 : osTMP.Printf("L1BGCPS_INTERPOL:\"%s\"", osFilename.c_str());
3524 : else
3525 0 : osTMP.Printf("L1BGCPS:\"%s\"", osFilename.c_str());
3526 3 : poOutDS->SetMetadataItem("X_DATASET", osTMP, "GEOLOCATION");
3527 3 : poOutDS->SetMetadataItem("X_BAND", "1", "GEOLOCATION");
3528 3 : poOutDS->SetMetadataItem("Y_DATASET", osTMP, "GEOLOCATION");
3529 3 : poOutDS->SetMetadataItem("Y_BAND", "2", "GEOLOCATION");
3530 :
3531 3 : if (bInterpol)
3532 : {
3533 3 : poOutDS->SetMetadataItem("PIXEL_OFFSET", "0", "GEOLOCATION");
3534 3 : poOutDS->SetMetadataItem("PIXEL_STEP", "1", "GEOLOCATION");
3535 : }
3536 : else
3537 : {
3538 0 : osTMP.Printf("%d", poDS->iGCPStart);
3539 0 : poOutDS->SetMetadataItem("PIXEL_OFFSET", osTMP, "GEOLOCATION");
3540 0 : osTMP.Printf("%d", poDS->iGCPStep);
3541 0 : poOutDS->SetMetadataItem("PIXEL_STEP", osTMP, "GEOLOCATION");
3542 : }
3543 :
3544 3 : poOutDS->SetMetadataItem("LINE_OFFSET", "0", "GEOLOCATION");
3545 3 : poOutDS->SetMetadataItem("LINE_STEP", "1", "GEOLOCATION");
3546 : }
3547 :
3548 3 : if (poOutDS != poDS)
3549 0 : return poOutDS;
3550 :
3551 3 : if (eL1BFormat == L1B_NOAA9)
3552 : {
3553 2 : char **papszSubdatasets = nullptr;
3554 2 : papszSubdatasets = CSLSetNameValue(
3555 : papszSubdatasets, "SUBDATASET_1_NAME",
3556 : CPLSPrintf("L1B_SOLAR_ZENITH_ANGLES:\"%s\"", osFilename.c_str()));
3557 2 : papszSubdatasets = CSLSetNameValue(
3558 : papszSubdatasets, "SUBDATASET_1_DESC", "Solar zenith angles");
3559 2 : poDS->SetMetadata(papszSubdatasets, "SUBDATASETS");
3560 2 : CSLDestroy(papszSubdatasets);
3561 : }
3562 : else
3563 : {
3564 1 : char **papszSubdatasets = nullptr;
3565 1 : papszSubdatasets = CSLSetNameValue(
3566 : papszSubdatasets, "SUBDATASET_1_NAME",
3567 : CPLSPrintf("L1B_ANGLES:\"%s\"", osFilename.c_str()));
3568 : papszSubdatasets =
3569 1 : CSLSetNameValue(papszSubdatasets, "SUBDATASET_1_DESC",
3570 : "Solar zenith angles, satellite zenith angles and "
3571 : "relative azimuth angles");
3572 :
3573 1 : if (poDS->iCLAVRStart > 0)
3574 : {
3575 0 : papszSubdatasets = CSLSetNameValue(
3576 : papszSubdatasets, "SUBDATASET_2_NAME",
3577 : CPLSPrintf("L1B_CLOUDS:\"%s\"", osFilename.c_str()));
3578 : papszSubdatasets =
3579 0 : CSLSetNameValue(papszSubdatasets, "SUBDATASET_2_DESC",
3580 : "Clouds from AVHRR (CLAVR)");
3581 : }
3582 :
3583 1 : poDS->SetMetadata(papszSubdatasets, "SUBDATASETS");
3584 1 : CSLDestroy(papszSubdatasets);
3585 : }
3586 :
3587 : /* -------------------------------------------------------------------- */
3588 : /* Create band information objects. */
3589 : /* -------------------------------------------------------------------- */
3590 : int iBand, i;
3591 :
3592 10 : for (iBand = 1, i = 0; iBand <= poDS->nBands; iBand++)
3593 : {
3594 7 : poDS->SetBand(iBand, new L1BRasterBand(poDS, iBand));
3595 :
3596 : // Channels descriptions
3597 7 : if (poDS->eSpacecraftID >= NOAA6 && poDS->eSpacecraftID <= METOP3)
3598 : {
3599 7 : if (!(i & 0x01) && poDS->iChannelsMask & 0x01)
3600 : {
3601 3 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[0]);
3602 3 : i |= 0x01;
3603 3 : continue;
3604 : }
3605 4 : if (!(i & 0x02) && poDS->iChannelsMask & 0x02)
3606 : {
3607 1 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[1]);
3608 1 : i |= 0x02;
3609 1 : continue;
3610 : }
3611 3 : if (!(i & 0x04) && poDS->iChannelsMask & 0x04)
3612 : {
3613 1 : if (poDS->eSpacecraftID >= NOAA15 &&
3614 1 : poDS->eSpacecraftID <= METOP3)
3615 1 : if (poDS->iInstrumentStatus & 0x0400)
3616 1 : poDS->GetRasterBand(iBand)->SetDescription(
3617 1 : apszBandDesc[7]);
3618 : else
3619 0 : poDS->GetRasterBand(iBand)->SetDescription(
3620 0 : apszBandDesc[6]);
3621 : else
3622 0 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[2]);
3623 1 : i |= 0x04;
3624 1 : continue;
3625 : }
3626 2 : if (!(i & 0x08) && poDS->iChannelsMask & 0x08)
3627 : {
3628 1 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[3]);
3629 1 : i |= 0x08;
3630 1 : continue;
3631 : }
3632 1 : if (!(i & 0x10) && poDS->iChannelsMask & 0x10)
3633 : {
3634 1 : if (poDS->eSpacecraftID == NOAA13) // 5 NOAA-13
3635 0 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[5]);
3636 1 : else if (poDS->eSpacecraftID == NOAA6 ||
3637 1 : poDS->eSpacecraftID == NOAA8 ||
3638 1 : poDS->eSpacecraftID == NOAA10) // 4 NOAA-6,-8,-10
3639 0 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[3]);
3640 : else
3641 1 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[4]);
3642 1 : i |= 0x10;
3643 1 : continue;
3644 : }
3645 : }
3646 : }
3647 :
3648 3 : if (poDS->bExposeMaskBand)
3649 1 : poDS->poMaskBand = new L1BMaskBand(poDS);
3650 :
3651 : /* -------------------------------------------------------------------- */
3652 : /* Initialize any PAM information. */
3653 : /* -------------------------------------------------------------------- */
3654 3 : poDS->SetDescription(poOpenInfo->pszFilename);
3655 3 : poDS->TryLoadXML();
3656 :
3657 : /* -------------------------------------------------------------------- */
3658 : /* Check for external overviews. */
3659 : /* -------------------------------------------------------------------- */
3660 6 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
3661 3 : poOpenInfo->GetSiblingFiles());
3662 :
3663 : /* -------------------------------------------------------------------- */
3664 : /* Fetch metadata in CSV file */
3665 : /* -------------------------------------------------------------------- */
3666 3 : if (CPLTestBool(CPLGetConfigOption("L1B_FETCH_METADATA", "NO")))
3667 : {
3668 0 : poDS->FetchMetadata();
3669 : }
3670 :
3671 3 : return poDS;
3672 :
3673 0 : bad:
3674 0 : delete poDS;
3675 0 : return nullptr;
3676 : }
3677 :
3678 : /************************************************************************/
3679 : /* GDALRegister_L1B() */
3680 : /************************************************************************/
3681 :
3682 1889 : void GDALRegister_L1B()
3683 :
3684 : {
3685 1889 : if (GDALGetDriverByName("L1B") != nullptr)
3686 282 : return;
3687 :
3688 1607 : GDALDriver *poDriver = new GDALDriver();
3689 :
3690 1607 : poDriver->SetDescription("L1B");
3691 1607 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
3692 1607 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
3693 1607 : "NOAA Polar Orbiter Level 1b Data Set");
3694 1607 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/l1b.html");
3695 1607 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
3696 1607 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
3697 :
3698 1607 : poDriver->pfnOpen = L1BDataset::Open;
3699 1607 : poDriver->pfnIdentify = L1BDataset::Identify;
3700 :
3701 1607 : GetGDALDriverManager()->RegisterDriver(poDriver);
3702 : }
|
