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 2 : TimeCode() : lYear(0), lDay(0), lMillisecond(0)
200 : {
201 2 : memset(szString, 0, sizeof(szString));
202 2 : }
203 :
204 2 : void SetYear(long year)
205 : {
206 2 : lYear = year;
207 2 : }
208 :
209 2 : void SetDay(long day)
210 : {
211 2 : lDay = day;
212 2 : }
213 :
214 2 : void SetMillisecond(long millisecond)
215 : {
216 2 : lMillisecond = millisecond;
217 2 : }
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 2 : char *PrintTime()
235 : {
236 2 : snprintf(szString, L1B_TIMECODE_LENGTH,
237 : "year: %ld, day: %ld, millisecond: %ld", lYear, lDay,
238 : lMillisecond);
239 2 : 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 = (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 5 : L1BRasterBand::L1BRasterBand(L1BDataset *poDSIn, int nBandIn)
431 :
432 : {
433 5 : poDS = poDSIn;
434 5 : nBand = nBandIn;
435 5 : eDataType = GDT_UInt16;
436 :
437 5 : nBlockXSize = poDS->GetRasterXSize();
438 5 : nBlockYSize = 1;
439 5 : }
440 :
441 : /************************************************************************/
442 : /* GetMaskBand() */
443 : /************************************************************************/
444 :
445 1 : GDALRasterBand *L1BRasterBand::GetMaskBand()
446 : {
447 1 : L1BDataset *poGDS = (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 = (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 = (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 1 : L1BDataset::L1BDataset(L1BFileFormat eL1BFormatIn)
577 : : eSource(UNKNOWN_STATION), eProcCenter(UNKNOWN_CENTER),
578 : // sStartTime
579 : // sStopTime
580 : bHighGCPDensityStrategy(
581 1 : 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 2 : bByteSwap(CPL_IS_LSB), bExposeMaskBand(FALSE), poMaskBand(nullptr)
592 : {
593 1 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
594 1 : 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 1 : }
602 :
603 : /************************************************************************/
604 : /* ~L1BDataset() */
605 : /************************************************************************/
606 :
607 2 : L1BDataset::~L1BDataset()
608 :
609 : {
610 1 : FlushCache(true);
611 :
612 1 : if (nGCPCount > 0)
613 : {
614 1 : GDALDeinitGCPs(nGCPCount, pasGCPList);
615 1 : CPLFree(pasGCPList);
616 : }
617 1 : if (fp != nullptr)
618 1 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
619 1 : delete poMaskBand;
620 2 : }
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 0 : void L1BDataset::FetchNOAA9TimeCode(TimeCode *psTime,
712 : const GByte *piRecordHeader,
713 : int *peLocationIndicator)
714 : {
715 : GUInt32 lTemp;
716 :
717 0 : lTemp = ((piRecordHeader[2] >> 1) & 0x7F);
718 0 : psTime->SetYear((lTemp > 77)
719 0 : ? (lTemp + 1900)
720 0 : : (lTemp + 2000)); // Avoid `Year 2000' problem
721 0 : psTime->SetDay((GUInt32)(piRecordHeader[2] & 0x01) << 8 |
722 0 : (GUInt32)piRecordHeader[3]);
723 0 : psTime->SetMillisecond(((GUInt32)(piRecordHeader[4] & 0x07) << 24) |
724 0 : ((GUInt32)piRecordHeader[5] << 16) |
725 0 : ((GUInt32)piRecordHeader[6] << 8) |
726 0 : (GUInt32)piRecordHeader[7]);
727 0 : if (peLocationIndicator)
728 : {
729 0 : *peLocationIndicator =
730 0 : ((piRecordHeader[8] & 0x02) == 0) ? ASCEND : DESCEND;
731 : }
732 0 : }
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 2 : void L1BDataset::FetchTimeCode(TimeCode *psTime, const void *pRecordHeader,
759 : int *peLocationIndicator) const
760 : {
761 2 : if (eSpacecraftID <= NOAA14)
762 : {
763 0 : FetchNOAA9TimeCode(psTime, (const GByte *)pRecordHeader,
764 : peLocationIndicator);
765 : }
766 : else
767 : {
768 2 : FetchNOAA15TimeCode(psTime, (const GByte *)pRecordHeader,
769 : peLocationIndicator);
770 : }
771 2 : }
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 1 : void L1BDataset::ProcessRecordHeaders()
853 : {
854 1 : void *pRecordHeader = CPLCalloc(1, nRecordDataStart);
855 :
856 1 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, nDataStartOffset, SEEK_SET));
857 1 : CPL_IGNORE_RET_VAL(VSIFReadL(pRecordHeader, 1, nRecordDataStart, fp));
858 :
859 1 : FetchTimeCode(&sStartTime, pRecordHeader, &eLocationIndicator);
860 :
861 1 : CPL_IGNORE_RET_VAL(VSIFSeekL(
862 1 : fp, nDataStartOffset + (nRasterYSize - 1) * nRecordSize, SEEK_SET));
863 1 : CPL_IGNORE_RET_VAL(VSIFReadL(pRecordHeader, 1, nRecordDataStart, fp));
864 :
865 1 : 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 1 : double dfLineStep = 0.0;
874 :
875 1 : if (bHighGCPDensityStrategy)
876 : {
877 1 : if (nRasterYSize < nGCPsPerLine)
878 : {
879 1 : 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 1 : if (nTargetLines > 1)
901 1 : dfLineStep = 1.0 * (nRasterYSize - 1) / (nTargetLines - 1);
902 :
903 : /* -------------------------------------------------------------------- */
904 : /* Initialize the GCP list. */
905 : /* -------------------------------------------------------------------- */
906 1 : const int nExpectedGCPs = nTargetLines * nGCPsPerLine;
907 1 : 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 1 : int iPrevLine = -1;
926 :
927 3 : 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 1 : 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 1 : CPLFree(pRecordHeader);
1001 :
1002 : /* -------------------------------------------------------------------- */
1003 : /* Set fetched information as metadata records */
1004 : /* -------------------------------------------------------------------- */
1005 : // Time of first scanline
1006 1 : SetMetadataItem("START", sStartTime.PrintTime());
1007 : // Time of last scanline
1008 1 : SetMetadataItem("STOP", sStopTime.PrintTime());
1009 : // AVHRR Earth location indication
1010 :
1011 1 : switch (eLocationIndicator)
1012 : {
1013 1 : case ASCEND:
1014 1 : SetMetadataItem("LOCATION", "Ascending");
1015 1 : 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 : const char *pszDir = CPLGetConfigOption("L1B_METADATA_DIRECTORY", nullptr);
1036 0 : if (pszDir == nullptr)
1037 : {
1038 0 : pszDir = CPLGetPath(GetDescription());
1039 0 : if (pszDir[0] == '\0')
1040 0 : pszDir = ".";
1041 : }
1042 : CPLString osMetadataFile(CPLSPrintf("%s/%s_metadata.csv", pszDir,
1043 0 : CPLGetFilename(GetDescription())));
1044 0 : VSILFILE *fpCSV = VSIFOpenL(osMetadataFile, "wb");
1045 0 : if (fpCSV == nullptr)
1046 : {
1047 0 : CPLError(CE_Failure, CPLE_AppDefined,
1048 : "Cannot create metadata file : %s", osMetadataFile.c_str());
1049 0 : return;
1050 : }
1051 :
1052 0 : CPL_IGNORE_RET_VAL(
1053 0 : VSIFPrintfL(fpCSV, "SCANLINE,NBLOCKYOFF,YEAR,DAY,MS_IN_DAY,"));
1054 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1055 : fpCSV, "FATAL_FLAG,TIME_ERROR,DATA_GAP,DATA_JITTER,INSUFFICIENT_DATA_"
1056 : "FOR_CAL,NO_EARTH_LOCATION,DESCEND,P_N_STATUS,"));
1057 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1058 : fpCSV, "BIT_SYNC_STATUS,SYNC_ERROR,FRAME_SYNC_ERROR,FLYWHEELING,BIT_"
1059 : "SLIPPAGE,C3_SBBC,C4_SBBC,C5_SBBC,"));
1060 0 : CPL_IGNORE_RET_VAL(
1061 0 : VSIFPrintfL(fpCSV, "TIP_PARITY_FRAME_1,TIP_PARITY_FRAME_2,TIP_PARITY_"
1062 : "FRAME_3,TIP_PARITY_FRAME_4,TIP_PARITY_FRAME_5,"));
1063 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "SYNC_ERRORS,"));
1064 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1065 : fpCSV, "CAL_SLOPE_C1,CAL_INTERCEPT_C1,CAL_SLOPE_C2,CAL_INTERCEPT_C2,"
1066 : "CAL_SLOPE_C3,CAL_INTERCEPT_C3,CAL_SLOPE_C4,CAL_INTERCEPT_C4,"
1067 : "CAL_SLOPE_C5,CAL_INTERCEPT_C5,"));
1068 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "NUM_SOLZENANGLES_EARTHLOCPNTS"));
1069 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));
1070 :
1071 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(nRecordDataStart);
1072 :
1073 0 : for (int nBlockYOff = 0; nBlockYOff < nRasterYSize; nBlockYOff++)
1074 : {
1075 : /* --------------------------------------------------------------------
1076 : */
1077 : /* Seek to data. */
1078 : /* --------------------------------------------------------------------
1079 : */
1080 0 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, GetLineOffset(nBlockYOff), SEEK_SET));
1081 :
1082 0 : CPL_IGNORE_RET_VAL(
1083 0 : VSIFReadL(pabyRecordHeader, 1, nRecordDataStart, fp));
1084 :
1085 0 : GUInt16 nScanlineNumber = GetUInt16(pabyRecordHeader);
1086 :
1087 0 : TimeCode timeCode;
1088 0 : FetchTimeCode(&timeCode, pabyRecordHeader, nullptr);
1089 :
1090 0 : CPL_IGNORE_RET_VAL(
1091 0 : VSIFPrintfL(fpCSV, "%d,%d,%d,%d,%d,", nScanlineNumber, nBlockYOff,
1092 0 : (int)timeCode.GetYear(), (int)timeCode.GetDay(),
1093 0 : (int)timeCode.GetMillisecond()));
1094 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1095 0 : fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,", (pabyRecordHeader[8] >> 7) & 1,
1096 0 : (pabyRecordHeader[8] >> 6) & 1, (pabyRecordHeader[8] >> 5) & 1,
1097 0 : (pabyRecordHeader[8] >> 4) & 1, (pabyRecordHeader[8] >> 3) & 1,
1098 0 : (pabyRecordHeader[8] >> 2) & 1, (pabyRecordHeader[8] >> 1) & 1,
1099 0 : (pabyRecordHeader[8] >> 0) & 1));
1100 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1101 0 : fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,", (pabyRecordHeader[9] >> 7) & 1,
1102 0 : (pabyRecordHeader[9] >> 6) & 1, (pabyRecordHeader[9] >> 5) & 1,
1103 0 : (pabyRecordHeader[9] >> 4) & 1, (pabyRecordHeader[9] >> 3) & 1,
1104 0 : (pabyRecordHeader[9] >> 2) & 1, (pabyRecordHeader[9] >> 1) & 1,
1105 0 : (pabyRecordHeader[9] >> 0) & 1));
1106 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1107 0 : fpCSV, "%d,%d,%d,%d,%d,", (pabyRecordHeader[10] >> 7) & 1,
1108 0 : (pabyRecordHeader[10] >> 6) & 1, (pabyRecordHeader[10] >> 5) & 1,
1109 0 : (pabyRecordHeader[10] >> 4) & 1, (pabyRecordHeader[10] >> 3) & 1));
1110 0 : CPL_IGNORE_RET_VAL(
1111 0 : VSIFPrintfL(fpCSV, "%d,", pabyRecordHeader[11] >> 2));
1112 : GInt32 i32;
1113 0 : for (int i = 0; i < 10; i++)
1114 : {
1115 0 : i32 = GetInt32(pabyRecordHeader + 12 + 4 * i);
1116 : /* Scales :
1117 : * http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c3/sec3-3.htm
1118 : */
1119 0 : if ((i % 2) == 0)
1120 0 : CPL_IGNORE_RET_VAL(
1121 0 : VSIFPrintfL(fpCSV, "%f,", i32 / pow(2.0, 30.0)));
1122 : else
1123 0 : CPL_IGNORE_RET_VAL(
1124 0 : VSIFPrintfL(fpCSV, "%f,", i32 / pow(2.0, 22.0)));
1125 : }
1126 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d", pabyRecordHeader[52]));
1127 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));
1128 : }
1129 :
1130 0 : CPLFree(pabyRecordHeader);
1131 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpCSV));
1132 : }
1133 :
1134 : /************************************************************************/
1135 : /* FetchMetadataNOAA15() */
1136 : /************************************************************************/
1137 :
1138 0 : void L1BDataset::FetchMetadataNOAA15()
1139 : {
1140 : int i, j;
1141 0 : const char *pszDir = CPLGetConfigOption("L1B_METADATA_DIRECTORY", nullptr);
1142 0 : if (pszDir == nullptr)
1143 : {
1144 0 : pszDir = CPLGetPath(GetDescription());
1145 0 : if (pszDir[0] == '\0')
1146 0 : pszDir = ".";
1147 : }
1148 : CPLString osMetadataFile(CPLSPrintf("%s/%s_metadata.csv", pszDir,
1149 0 : CPLGetFilename(GetDescription())));
1150 0 : VSILFILE *fpCSV = VSIFOpenL(osMetadataFile, "wb");
1151 0 : if (fpCSV == nullptr)
1152 : {
1153 0 : CPLError(CE_Failure, CPLE_AppDefined,
1154 : "Cannot create metadata file : %s", osMetadataFile.c_str());
1155 0 : return;
1156 : }
1157 :
1158 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1159 : fpCSV, "SCANLINE,NBLOCKYOFF,YEAR,DAY,MS_IN_DAY,SAT_CLOCK_DRIF_DELTA,"
1160 : "SOUTHBOUND,SCANTIME_CORRECTED,C3_SELECT,"));
1161 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1162 : fpCSV,
1163 : "FATAL_FLAG,TIME_ERROR,DATA_GAP,INSUFFICIENT_DATA_FOR_CAL,"
1164 : "NO_EARTH_LOCATION,FIRST_GOOD_TIME_AFTER_CLOCK_UPDATE,"
1165 : "INSTRUMENT_STATUS_CHANGED,SYNC_LOCK_DROPPED,"
1166 : "FRAME_SYNC_ERROR,FRAME_SYNC_DROPPED_LOCK,FLYWHEELING,"
1167 : "BIT_SLIPPAGE,TIP_PARITY_ERROR,REFLECTED_SUNLIGHT_C3B,"
1168 : "REFLECTED_SUNLIGHT_C4,REFLECTED_SUNLIGHT_C5,RESYNC,P_N_STATUS,"));
1169 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1170 : fpCSV, "BAD_TIME_CAN_BE_INFERRED,BAD_TIME_CANNOT_BE_INFERRED,"
1171 : "TIME_DISCONTINUITY,REPEAT_SCAN_TIME,"));
1172 0 : CPL_IGNORE_RET_VAL(
1173 0 : VSIFPrintfL(fpCSV, "UNCALIBRATED_BAD_TIME,CALIBRATED_FEWER_SCANLINES,"
1174 : "UNCALIBRATED_BAD_PRT,CALIBRATED_MARGINAL_PRT,"
1175 : "UNCALIBRATED_CHANNELS,"));
1176 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1177 : fpCSV, "NO_EARTH_LOC_BAD_TIME,EARTH_LOC_QUESTIONABLE_TIME,"
1178 : "EARTH_LOC_QUESTIONABLE,EARTH_LOC_VERY_QUESTIONABLE,"));
1179 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1180 : fpCSV,
1181 : "C3B_UNCALIBRATED,C3B_QUESTIONABLE,C3B_ALL_BLACKBODY,"
1182 : "C3B_ALL_SPACEVIEW,C3B_MARGINAL_BLACKBODY,C3B_MARGINAL_SPACEVIEW,"));
1183 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1184 : fpCSV,
1185 : "C4_UNCALIBRATED,C4_QUESTIONABLE,C4_ALL_BLACKBODY,"
1186 : "C4_ALL_SPACEVIEW,C4_MARGINAL_BLACKBODY,C4_MARGINAL_SPACEVIEW,"));
1187 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1188 : fpCSV,
1189 : "C5_UNCALIBRATED,C5_QUESTIONABLE,C5_ALL_BLACKBODY,"
1190 : "C5_ALL_SPACEVIEW,C5_MARGINAL_BLACKBODY,C5_MARGINAL_SPACEVIEW,"));
1191 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "BIT_ERRORS,"));
1192 0 : for (i = 0; i < 3; i++)
1193 : {
1194 0 : const char *pszChannel = (i == 0) ? "C1" : (i == 1) ? "C2" : "C3A";
1195 0 : for (j = 0; j < 3; j++)
1196 : {
1197 0 : const char *pszType = (j == 0) ? "OP"
1198 0 : : (j == 1) ? "TEST"
1199 : : "PRELAUNCH";
1200 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_SLOPE_1,",
1201 : pszType, pszChannel));
1202 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_INTERCEPT_1,",
1203 : pszType, pszChannel));
1204 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_SLOPE_2,",
1205 : pszType, pszChannel));
1206 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_INTERCEPT_2,",
1207 : pszType, pszChannel));
1208 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "VIS_%s_CAL_%s_INTERSECTION,",
1209 : pszType, pszChannel));
1210 : }
1211 : }
1212 0 : for (i = 0; i < 3; i++)
1213 : {
1214 0 : const char *pszChannel = (i == 0) ? "C3B" : (i == 1) ? "C4" : "C5";
1215 0 : for (j = 0; j < 2; j++)
1216 : {
1217 0 : const char *pszType = (j == 0) ? "OP" : "TEST";
1218 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "IR_%s_CAL_%s_COEFF_1,",
1219 : pszType, pszChannel));
1220 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "IR_%s_CAL_%s_COEFF_2,",
1221 : pszType, pszChannel));
1222 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "IR_%s_CAL_%s_COEFF_3,",
1223 : pszType, pszChannel));
1224 : }
1225 : }
1226 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1227 : fpCSV, "EARTH_LOC_CORR_TIP_EULER,EARTH_LOC_IND,"
1228 : "SPACECRAFT_ATT_CTRL,ATT_SMODE,ATT_PASSIVE_WHEEL_TEST,"
1229 : "TIME_TIP_EULER,TIP_EULER_ROLL,TIP_EULER_PITCH,TIP_EULER_YAW,"
1230 : "SPACECRAFT_ALT"));
1231 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));
1232 :
1233 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(nRecordDataStart);
1234 : GInt16 i16;
1235 : GUInt16 n16;
1236 : GInt32 i32;
1237 : GUInt32 n32;
1238 :
1239 0 : for (int nBlockYOff = 0; nBlockYOff < nRasterYSize; nBlockYOff++)
1240 : {
1241 : /* --------------------------------------------------------------------
1242 : */
1243 : /* Seek to data. */
1244 : /* --------------------------------------------------------------------
1245 : */
1246 0 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, GetLineOffset(nBlockYOff), SEEK_SET));
1247 :
1248 0 : CPL_IGNORE_RET_VAL(
1249 0 : VSIFReadL(pabyRecordHeader, 1, nRecordDataStart, fp));
1250 :
1251 0 : GUInt16 nScanlineNumber = GetUInt16(pabyRecordHeader);
1252 :
1253 0 : TimeCode timeCode;
1254 0 : FetchTimeCode(&timeCode, pabyRecordHeader, nullptr);
1255 :
1256 : /* Clock drift delta */
1257 0 : i16 = GetInt16(pabyRecordHeader + 6);
1258 : /* Scanline bit field */
1259 0 : n16 = GetInt16(pabyRecordHeader + 12);
1260 :
1261 0 : CPL_IGNORE_RET_VAL(
1262 0 : VSIFPrintfL(fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,%d,", nScanlineNumber,
1263 0 : nBlockYOff, (int)timeCode.GetYear(),
1264 0 : (int)timeCode.GetDay(), (int)timeCode.GetMillisecond(),
1265 0 : i16, (n16 >> 15) & 1, (n16 >> 14) & 1, (n16)&3));
1266 :
1267 0 : n32 = GetUInt32(pabyRecordHeader + 24);
1268 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1269 : fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,",
1270 0 : (n32 >> 31) & 1, (n32 >> 30) & 1, (n32 >> 29) & 1, (n32 >> 28) & 1,
1271 0 : (n32 >> 27) & 1, (n32 >> 26) & 1, (n32 >> 25) & 1, (n32 >> 24) & 1,
1272 0 : (n32 >> 23) & 1, (n32 >> 22) & 1, (n32 >> 21) & 1, (n32 >> 20) & 1,
1273 0 : (n32 >> 8) & 1, (n32 >> 6) & 3, (n32 >> 4) & 3, (n32 >> 2) & 3,
1274 0 : (n32 >> 1) & 1, (n32 >> 0) & 1));
1275 :
1276 0 : n32 = GetUInt32(pabyRecordHeader + 28);
1277 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1278 0 : fpCSV, "%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,", (n32 >> 23) & 1,
1279 0 : (n32 >> 22) & 1, (n32 >> 21) & 1, (n32 >> 20) & 1, (n32 >> 15) & 1,
1280 0 : (n32 >> 14) & 1, (n32 >> 13) & 1, (n32 >> 12) & 1, (n32 >> 11) & 1,
1281 0 : (n32 >> 7) & 1, (n32 >> 6) & 1, (n32 >> 5) & 1, (n32 >> 4) & 1));
1282 :
1283 0 : for (i = 0; i < 3; i++)
1284 : {
1285 0 : n16 = GetUInt16(pabyRecordHeader + 32 + 2 * i);
1286 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d,%d,%d,%d,%d,%d,",
1287 0 : (n16 >> 7) & 1, (n16 >> 6) & 1,
1288 0 : (n16 >> 5) & 1, (n16 >> 4) & 1,
1289 0 : (n16 >> 2) & 1, (n16 >> 1) & 1));
1290 : }
1291 :
1292 : /* Bit errors */
1293 0 : n16 = GetUInt16(pabyRecordHeader + 38);
1294 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d,", n16));
1295 :
1296 0 : int nOffset = 48;
1297 0 : for (i = 0; i < 3; i++)
1298 : {
1299 0 : for (j = 0; j < 3; j++)
1300 : {
1301 0 : i32 = GetInt32(pabyRecordHeader + nOffset);
1302 0 : nOffset += 4;
1303 0 : CPL_IGNORE_RET_VAL(
1304 0 : VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 7.0)));
1305 0 : i32 = GetInt32(pabyRecordHeader + nOffset);
1306 0 : nOffset += 4;
1307 0 : CPL_IGNORE_RET_VAL(
1308 0 : VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 6.0)));
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(VSIFPrintfL(fpCSV, "%d,", i32));
1320 : }
1321 : }
1322 0 : for (i = 0; i < 18; i++)
1323 : {
1324 0 : i32 = GetInt32(pabyRecordHeader + nOffset);
1325 0 : nOffset += 4;
1326 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%f,", i32 / pow(10.0, 6.0)));
1327 : }
1328 :
1329 0 : n32 = GetUInt32(pabyRecordHeader + 312);
1330 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(
1331 0 : fpCSV, "%d,%d,%d,%d,%d,", (n32 >> 16) & 1, (n32 >> 12) & 15,
1332 0 : (n32 >> 8) & 15, (n32 >> 4) & 15, (n32 >> 0) & 15));
1333 :
1334 0 : n32 = GetUInt32(pabyRecordHeader + 316);
1335 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%d,", n32));
1336 :
1337 0 : for (i = 0; i < 3; i++)
1338 : {
1339 0 : i16 = GetUInt16(pabyRecordHeader + 320 + 2 * i);
1340 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%f,", i16 / pow(10.0, 3.0)));
1341 : }
1342 :
1343 0 : n16 = GetUInt16(pabyRecordHeader + 326);
1344 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "%f", n16 / pow(10.0, 1.0)));
1345 :
1346 0 : CPL_IGNORE_RET_VAL(VSIFPrintfL(fpCSV, "\n"));
1347 : }
1348 :
1349 0 : CPLFree(pabyRecordHeader);
1350 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpCSV));
1351 : }
1352 :
1353 : /************************************************************************/
1354 : /* EBCDICToASCII */
1355 : /************************************************************************/
1356 :
1357 : constexpr GByte EBCDICToASCII[] = {
1358 : 0x00, 0x01, 0x02, 0x03, 0x9C, 0x09, 0x86, 0x7F, 0x97, 0x8D, 0x8E, 0x0B,
1359 : 0x0C, 0x0D, 0x0E, 0x0F, 0x10, 0x11, 0x12, 0x13, 0x9D, 0x85, 0x08, 0x87,
1360 : 0x18, 0x19, 0x92, 0x8F, 0x1C, 0x1D, 0x1E, 0x1F, 0x80, 0x81, 0x82, 0x83,
1361 : 0x84, 0x0A, 0x17, 0x1B, 0x88, 0x89, 0x8A, 0x8B, 0x8C, 0x05, 0x06, 0x07,
1362 : 0x90, 0x91, 0x16, 0x93, 0x94, 0x95, 0x96, 0x04, 0x98, 0x99, 0x9A, 0x9B,
1363 : 0x14, 0x15, 0x9E, 0x1A, 0x20, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1364 : 0x00, 0x00, 0xA2, 0x2E, 0x3C, 0x28, 0x2B, 0x7C, 0x26, 0x00, 0x00, 0x00,
1365 : 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x21, 0x24, 0x2A, 0x29, 0x3B, 0xAC,
1366 : 0x2D, 0x2F, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xA6, 0x2C,
1367 : 0x25, 0x5F, 0x3E, 0x3F, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1368 : 0x00, 0x60, 0x3A, 0x23, 0x40, 0x27, 0x3D, 0x22, 0x00, 0x61, 0x62, 0x63,
1369 : 0x64, 0x65, 0x66, 0x67, 0x68, 0x69, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1370 : 0x00, 0x6A, 0x6B, 0x6C, 0x6D, 0x6E, 0x6F, 0x70, 0x71, 0x72, 0x00, 0x00,
1371 : 0x00, 0x00, 0x00, 0x00, 0x00, 0x7E, 0x73, 0x74, 0x75, 0x76, 0x77, 0x78,
1372 : 0x79, 0x7A, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1373 : 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1374 : 0x7B, 0x41, 0x42, 0x43, 0x44, 0x45, 0x46, 0x47, 0x48, 0x49, 0x00, 0x00,
1375 : 0x00, 0x00, 0x00, 0x00, 0x7D, 0x4A, 0x4B, 0x4C, 0x4D, 0x4E, 0x4F, 0x50,
1376 : 0x51, 0x52, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x5C, 0x00, 0x53, 0x54,
1377 : 0x55, 0x56, 0x57, 0x58, 0x59, 0x5A, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
1378 : 0x30, 0x31, 0x32, 0x33, 0x34, 0x35, 0x36, 0x37, 0x38, 0x39, 0x00, 0x00,
1379 : 0x00, 0x00, 0x00, 0x9F,
1380 : };
1381 :
1382 : /************************************************************************/
1383 : /* ProcessDatasetHeader() */
1384 : /************************************************************************/
1385 :
1386 1 : CPLErr L1BDataset::ProcessDatasetHeader(const char *pszFilename)
1387 : {
1388 : char szDatasetName[L1B_DATASET_NAME_SIZE + 1];
1389 :
1390 1 : if (eL1BFormat == L1B_NOAA9)
1391 : {
1392 : GByte abyTBMHeader[L1B_NOAA9_HEADER_SIZE];
1393 :
1394 0 : if (VSIFSeekL(fp, 0, SEEK_SET) < 0 ||
1395 0 : VSIFReadL(abyTBMHeader, 1, L1B_NOAA9_HEADER_SIZE, fp) <
1396 : L1B_NOAA9_HEADER_SIZE)
1397 : {
1398 0 : CPLDebug("L1B", "Can't read NOAA-9/14 TBM header.");
1399 0 : return CE_Failure;
1400 : }
1401 :
1402 : // If dataset name in EBCDIC, decode it in ASCII
1403 0 : if (*(abyTBMHeader + 8 + 25) == 'K' &&
1404 0 : *(abyTBMHeader + 8 + 30) == 'K' &&
1405 0 : *(abyTBMHeader + 8 + 33) == 'K' &&
1406 0 : *(abyTBMHeader + 8 + 40) == 'K' &&
1407 0 : *(abyTBMHeader + 8 + 46) == 'K' &&
1408 0 : *(abyTBMHeader + 8 + 52) == 'K' && *(abyTBMHeader + 8 + 61) == 'K')
1409 : {
1410 0 : for (int i = 0; i < L1B_DATASET_NAME_SIZE; i++)
1411 0 : abyTBMHeader[L1B_NOAA9_HDR_NAME_OFF + i] =
1412 0 : EBCDICToASCII[abyTBMHeader[L1B_NOAA9_HDR_NAME_OFF + i]];
1413 : }
1414 :
1415 : // Fetch dataset name. NOAA-9/14 datasets contain the names in TBM
1416 : // header only, so read it there.
1417 0 : memcpy(szDatasetName, abyTBMHeader + L1B_NOAA9_HDR_NAME_OFF,
1418 : L1B_DATASET_NAME_SIZE);
1419 0 : szDatasetName[L1B_DATASET_NAME_SIZE] = '\0';
1420 :
1421 : // Deal with a few NOAA <= 9 datasets with no dataset name in TBM header
1422 0 : if (memcmp(szDatasetName,
1423 : "\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"
1424 : "\0\0\0\0\0\0\0\0\0\0\0\0\0",
1425 0 : L1B_DATASET_NAME_SIZE) == 0 &&
1426 0 : strlen(pszFilename) == L1B_DATASET_NAME_SIZE)
1427 : {
1428 0 : strcpy(szDatasetName, pszFilename);
1429 : }
1430 :
1431 : // Determine processing center where the dataset was created
1432 0 : if (STARTS_WITH_CI(szDatasetName, "CMS"))
1433 0 : eProcCenter = CMS;
1434 0 : else if (STARTS_WITH_CI(szDatasetName, "DSS"))
1435 0 : eProcCenter = DSS;
1436 0 : else if (STARTS_WITH_CI(szDatasetName, "NSS"))
1437 0 : eProcCenter = NSS;
1438 0 : else if (STARTS_WITH_CI(szDatasetName, "UKM"))
1439 0 : eProcCenter = UKM;
1440 : else
1441 0 : eProcCenter = UNKNOWN_CENTER;
1442 :
1443 : // Determine number of bands
1444 : int i;
1445 0 : for (i = 0; i < L1B_NOAA9_HDR_CHAN_SIZE; i++)
1446 : {
1447 0 : if (abyTBMHeader[L1B_NOAA9_HDR_CHAN_OFF + i] == 1 ||
1448 0 : abyTBMHeader[L1B_NOAA9_HDR_CHAN_OFF + i] == 'Y')
1449 : {
1450 0 : nBands++;
1451 0 : iChannelsMask |= (1 << i);
1452 : }
1453 : }
1454 0 : if (nBands == 0 || nBands > 5)
1455 : {
1456 0 : nBands = 5;
1457 0 : iChannelsMask = 0x1F;
1458 : }
1459 :
1460 : // Determine data format (10-bit packed or 8/16-bit unpacked)
1461 0 : if (STARTS_WITH_CI((const char *)abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF,
1462 : "10"))
1463 0 : iDataFormat = PACKED10BIT;
1464 0 : else if (STARTS_WITH_CI(
1465 : (const char *)abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF, "16"))
1466 0 : iDataFormat = UNPACKED16BIT;
1467 0 : else if (STARTS_WITH_CI(
1468 : (const char *)abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF, "08"))
1469 0 : iDataFormat = UNPACKED8BIT;
1470 0 : else if (STARTS_WITH_CI((const char *)abyTBMHeader +
1471 : L1B_NOAA9_HDR_WORD_OFF,
1472 0 : " ") ||
1473 0 : abyTBMHeader[L1B_NOAA9_HDR_WORD_OFF] == '\0')
1474 : /* Empty string can be found in the following samples :
1475 : http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/franh.1b (10
1476 : bit) http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/frang.1b (10
1477 : bit) http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/calfilel.1b
1478 : (16 bit)
1479 : http://www2.ncdc.noaa.gov/docs/podug/data/avhrr/rapnzg.1b (16
1480 : bit)
1481 : ftp://ftp.sat.dundee.ac.uk/misc/testdata/noaa12/hrptnoaa1b.dat
1482 : (10 bit)
1483 : */
1484 0 : bGuessDataFormat = TRUE;
1485 : else
1486 : {
1487 : #ifdef DEBUG
1488 0 : CPLDebug("L1B", "Unknown data format \"%.2s\".",
1489 : abyTBMHeader + L1B_NOAA9_HDR_WORD_OFF);
1490 : #endif
1491 0 : return CE_Failure;
1492 : }
1493 :
1494 : // Now read the dataset header record
1495 : GByte abyRecHeader[L1B_NOAA9_HDR_REC_SIZE];
1496 0 : if (VSIFSeekL(fp, L1B_NOAA9_HEADER_SIZE, SEEK_SET) < 0 ||
1497 0 : VSIFReadL(abyRecHeader, 1, L1B_NOAA9_HDR_REC_SIZE, fp) <
1498 : L1B_NOAA9_HDR_REC_SIZE)
1499 : {
1500 0 : CPLDebug("L1B", "Can't read NOAA-9/14 record header.");
1501 0 : return CE_Failure;
1502 : }
1503 :
1504 : // Determine the spacecraft name
1505 0 : switch (abyRecHeader[L1B_NOAA9_HDR_REC_ID_OFF])
1506 : {
1507 0 : case 4:
1508 0 : eSpacecraftID = NOAA7;
1509 0 : break;
1510 0 : case 6:
1511 0 : eSpacecraftID = NOAA8;
1512 0 : break;
1513 0 : case 7:
1514 0 : eSpacecraftID = NOAA9;
1515 0 : break;
1516 0 : case 8:
1517 0 : eSpacecraftID = NOAA10;
1518 0 : break;
1519 0 : case 1:
1520 : {
1521 : /* We could also use the time code to determine TIROS-N */
1522 0 : if (strlen(pszFilename) == L1B_DATASET_NAME_SIZE &&
1523 0 : STARTS_WITH(pszFilename + 8, ".TN."))
1524 0 : eSpacecraftID = TIROSN;
1525 : else
1526 0 : eSpacecraftID = NOAA11;
1527 0 : break;
1528 : }
1529 0 : case 5:
1530 0 : eSpacecraftID = NOAA12;
1531 0 : break;
1532 0 : case 2:
1533 : {
1534 : /* We could also use the time code to determine NOAA6 */
1535 0 : if (strlen(pszFilename) == L1B_DATASET_NAME_SIZE &&
1536 0 : STARTS_WITH(pszFilename + 8, ".NA."))
1537 0 : eSpacecraftID = NOAA6;
1538 : else
1539 0 : eSpacecraftID = NOAA13;
1540 0 : break;
1541 : }
1542 0 : case 3:
1543 0 : eSpacecraftID = NOAA14;
1544 0 : break;
1545 0 : default:
1546 0 : CPLError(CE_Warning, CPLE_AppDefined,
1547 : "Unknown spacecraft ID \"%d\".",
1548 0 : abyRecHeader[L1B_NOAA9_HDR_REC_ID_OFF]);
1549 :
1550 0 : eSpacecraftID = NOAA9_UNKNOWN;
1551 0 : break;
1552 : }
1553 :
1554 : // Determine the product data type
1555 0 : int iWord = abyRecHeader[L1B_NOAA9_HDR_REC_PROD_OFF] >> 4;
1556 0 : switch (iWord)
1557 : {
1558 0 : case 1:
1559 0 : eProductType = LAC;
1560 0 : break;
1561 0 : case 2:
1562 0 : eProductType = GAC;
1563 0 : break;
1564 0 : case 3:
1565 0 : eProductType = HRPT;
1566 0 : break;
1567 0 : default:
1568 : #ifdef DEBUG
1569 0 : CPLDebug("L1B", "Unknown product type \"%d\".", iWord);
1570 : #endif
1571 0 : return CE_Failure;
1572 : }
1573 :
1574 : // Determine receiving station name
1575 0 : iWord = (abyRecHeader[L1B_NOAA9_HDR_REC_DSTAT_OFF] & 0x60) >> 5;
1576 0 : switch (iWord)
1577 : {
1578 0 : case 1:
1579 0 : eSource = GC;
1580 0 : break;
1581 0 : case 2:
1582 0 : eSource = WI;
1583 0 : break;
1584 0 : case 3:
1585 0 : eSource = SO;
1586 0 : break;
1587 0 : default:
1588 0 : eSource = UNKNOWN_STATION;
1589 0 : break;
1590 : }
1591 : }
1592 :
1593 1 : else if (eL1BFormat == L1B_NOAA15 || eL1BFormat == L1B_NOAA15_NOHDR)
1594 : {
1595 1 : if (eL1BFormat == L1B_NOAA15)
1596 : {
1597 : GByte abyARSHeader[L1B_NOAA15_HEADER_SIZE];
1598 :
1599 0 : if (VSIFSeekL(fp, 0, SEEK_SET) < 0 ||
1600 0 : VSIFReadL(abyARSHeader, 1, L1B_NOAA15_HEADER_SIZE, fp) <
1601 : L1B_NOAA15_HEADER_SIZE)
1602 : {
1603 0 : CPLDebug("L1B", "Can't read NOAA-15 ARS header.");
1604 0 : return CE_Failure;
1605 : }
1606 :
1607 : // Determine number of bands
1608 : int i;
1609 0 : for (i = 0; i < L1B_NOAA15_HDR_CHAN_SIZE; i++)
1610 : {
1611 0 : if (abyARSHeader[L1B_NOAA15_HDR_CHAN_OFF + i] == 1 ||
1612 0 : abyARSHeader[L1B_NOAA15_HDR_CHAN_OFF + i] == 'Y')
1613 : {
1614 0 : nBands++;
1615 0 : iChannelsMask |= (1 << i);
1616 : }
1617 : }
1618 0 : if (nBands == 0 || nBands > 5)
1619 : {
1620 0 : nBands = 5;
1621 0 : iChannelsMask = 0x1F;
1622 : }
1623 :
1624 : // Determine data format (10-bit packed or 8/16-bit unpacked)
1625 0 : if (STARTS_WITH_CI(
1626 : (const char *)abyARSHeader + L1B_NOAA15_HDR_WORD_OFF, "10"))
1627 0 : iDataFormat = PACKED10BIT;
1628 0 : else if (STARTS_WITH_CI((const char *)abyARSHeader +
1629 : L1B_NOAA15_HDR_WORD_OFF,
1630 : "16"))
1631 0 : iDataFormat = UNPACKED16BIT;
1632 0 : else if (STARTS_WITH_CI((const char *)abyARSHeader +
1633 : L1B_NOAA15_HDR_WORD_OFF,
1634 : "08"))
1635 0 : iDataFormat = UNPACKED8BIT;
1636 : else
1637 : {
1638 : #ifdef DEBUG
1639 0 : CPLDebug("L1B", "Unknown data format \"%.2s\".",
1640 : abyARSHeader + L1B_NOAA9_HDR_WORD_OFF);
1641 : #endif
1642 0 : return CE_Failure;
1643 : }
1644 : }
1645 : else
1646 : {
1647 1 : nBands = 5;
1648 1 : iChannelsMask = 0x1F;
1649 1 : iDataFormat = PACKED10BIT;
1650 : }
1651 :
1652 : // Now read the dataset header record
1653 : GByte abyRecHeader[L1B_NOAA15_HDR_REC_SIZE];
1654 1 : if (VSIFSeekL(fp,
1655 1 : (eL1BFormat == L1B_NOAA15) ? L1B_NOAA15_HEADER_SIZE : 0,
1656 2 : SEEK_SET) < 0 ||
1657 1 : VSIFReadL(abyRecHeader, 1, L1B_NOAA15_HDR_REC_SIZE, fp) <
1658 : L1B_NOAA15_HDR_REC_SIZE)
1659 : {
1660 0 : CPLDebug("L1B", "Can't read NOAA-9/14 record header.");
1661 0 : return CE_Failure;
1662 : }
1663 :
1664 : // Fetch dataset name
1665 1 : memcpy(szDatasetName, abyRecHeader + L1B_NOAA15_HDR_REC_NAME_OFF,
1666 : L1B_DATASET_NAME_SIZE);
1667 1 : szDatasetName[L1B_DATASET_NAME_SIZE] = '\0';
1668 :
1669 : // Determine processing center where the dataset was created
1670 1 : if (STARTS_WITH_CI((const char *)abyRecHeader +
1671 : L1B_NOAA15_HDR_REC_SITE_OFF,
1672 : "CMS"))
1673 0 : eProcCenter = CMS;
1674 1 : else if (STARTS_WITH_CI((const char *)abyRecHeader +
1675 : L1B_NOAA15_HDR_REC_SITE_OFF,
1676 : "DSS"))
1677 0 : eProcCenter = DSS;
1678 1 : else if (STARTS_WITH_CI((const char *)abyRecHeader +
1679 : L1B_NOAA15_HDR_REC_SITE_OFF,
1680 : "NSS"))
1681 0 : eProcCenter = NSS;
1682 1 : else if (STARTS_WITH_CI((const char *)abyRecHeader +
1683 : L1B_NOAA15_HDR_REC_SITE_OFF,
1684 : "UKM"))
1685 0 : eProcCenter = UKM;
1686 : else
1687 1 : eProcCenter = UNKNOWN_CENTER;
1688 :
1689 : int nFormatVersionYear, nFormatVersionDayOfYear, nHeaderRecCount;
1690 :
1691 : /* Some products from NOAA-18 and NOAA-19 coming from 'ess' processing
1692 : * station */
1693 : /* have little-endian ordering. Try to detect it with some consistency
1694 : * checks */
1695 1 : int i = 0;
1696 1 : do
1697 : {
1698 2 : nFormatVersionYear = GetUInt16(
1699 : abyRecHeader + L1B_NOAA15_HDR_REC_FORMAT_VERSION_YEAR_OFF);
1700 2 : nFormatVersionDayOfYear = GetUInt16(
1701 : abyRecHeader + L1B_NOAA15_HDR_REC_FORMAT_VERSION_DAY_OFF);
1702 2 : nHeaderRecCount =
1703 2 : GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_HDR_REC_COUNT_OFF);
1704 2 : if (i == 2)
1705 0 : break;
1706 2 : if (!(nFormatVersionYear >= 1980 && nFormatVersionYear <= 2100) &&
1707 1 : !(nFormatVersionDayOfYear <= 366) && !(nHeaderRecCount == 1))
1708 : {
1709 1 : if (i == 0)
1710 1 : CPLDebug("L1B", "Trying little-endian ordering");
1711 : else
1712 0 : CPLDebug("L1B", "Not completely convincing... Returning to "
1713 : "big-endian order");
1714 1 : bByteSwap = !bByteSwap;
1715 : }
1716 : else
1717 : break;
1718 1 : i++;
1719 1 : } while (i <= 2);
1720 1 : nRecordSizeFromHeader =
1721 1 : GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_LOGICAL_REC_LENGTH_OFF);
1722 : int nFormatVersion =
1723 1 : GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_FORMAT_VERSION_OFF);
1724 1 : CPLDebug("L1B", "NOAA Level 1b Format Version Number = %d",
1725 : nFormatVersion);
1726 1 : CPLDebug("L1B", "Level 1b Format Version Year = %d",
1727 : nFormatVersionYear);
1728 1 : CPLDebug("L1B", "Level 1b Format Version Day of Year = %d",
1729 : nFormatVersionDayOfYear);
1730 1 : CPLDebug("L1B",
1731 : "Logical Record Length of source Level 1b data set prior to "
1732 : "processing = %d",
1733 : nRecordSizeFromHeader);
1734 : int nBlockSize =
1735 1 : GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_BLOCK_SIZE_OFF);
1736 1 : CPLDebug(
1737 : "L1B",
1738 : "Block Size of source Level 1b data set prior to processing = %d",
1739 : nBlockSize);
1740 1 : CPLDebug("L1B", "Count of Header Records in this Data Set = %d",
1741 : nHeaderRecCount);
1742 :
1743 : int nDataRecordCount =
1744 1 : GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_DATA_RECORD_COUNT_OFF);
1745 1 : CPLDebug("L1B", "Count of Data Records = %d", nDataRecordCount);
1746 :
1747 1 : int nCalibratedScanlineCount = GetUInt16(
1748 1 : abyRecHeader + L1B_NOAA15_HDR_REC_CALIBRATED_SCANLINE_COUNT_OFF);
1749 1 : CPLDebug("L1B", "Count of Calibrated, Earth Located Scan Lines = %d",
1750 : nCalibratedScanlineCount);
1751 :
1752 1 : int nMissingScanlineCount = GetUInt16(
1753 1 : abyRecHeader + L1B_NOAA15_HDR_REC_MISSING_SCANLINE_COUNT_OFF);
1754 1 : CPLDebug("L1B", "Count of Missing Scan Lines = %d",
1755 : nMissingScanlineCount);
1756 1 : if (nMissingScanlineCount != 0)
1757 1 : bExposeMaskBand = TRUE;
1758 :
1759 : char szEllipsoid[8 + 1];
1760 1 : memcpy(szEllipsoid, abyRecHeader + L1B_NOAA15_HDR_REC_ELLIPSOID_OFF, 8);
1761 1 : szEllipsoid[8] = '\0';
1762 1 : CPLDebug("L1B", "Reference Ellipsoid Model ID = '%s'", szEllipsoid);
1763 1 : if (EQUAL(szEllipsoid, "WGS-84 "))
1764 : {
1765 0 : m_oGCPSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
1766 : }
1767 1 : else if (EQUAL(szEllipsoid, " GRS 80"))
1768 : {
1769 1 : m_oGCPSRS.importFromWkt(
1770 : "GEOGCS[\"GRS 1980(IUGG, "
1771 : "1980)\",DATUM[\"unknown\",SPHEROID[\"GRS80\",6378137,298."
1772 : "257222101],TOWGS84[0,0,0,0,0,0,0]],PRIMEM[\"Greenwich\",0],"
1773 : "UNIT[\"degree\",0.0174532925199433]]");
1774 : }
1775 :
1776 : // Determine the spacecraft name
1777 : // See
1778 : // http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c8/sec83132-2.htm
1779 1 : int iWord = GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_ID_OFF);
1780 1 : switch (iWord)
1781 : {
1782 0 : case 2:
1783 0 : eSpacecraftID = NOAA16;
1784 0 : break;
1785 0 : case 4:
1786 0 : eSpacecraftID = NOAA15;
1787 0 : break;
1788 0 : case 6:
1789 0 : eSpacecraftID = NOAA17;
1790 0 : break;
1791 0 : case 7:
1792 0 : eSpacecraftID = NOAA18;
1793 0 : break;
1794 1 : case 8:
1795 1 : eSpacecraftID = NOAA19;
1796 1 : break;
1797 0 : case 11:
1798 0 : eSpacecraftID = METOP1;
1799 0 : break;
1800 0 : case 12:
1801 0 : eSpacecraftID = METOP2;
1802 0 : break;
1803 : // METOP3 is not documented yet
1804 0 : case 13:
1805 0 : eSpacecraftID = METOP3;
1806 0 : break;
1807 0 : case 14:
1808 0 : eSpacecraftID = METOP3;
1809 0 : break;
1810 0 : default:
1811 : #ifdef DEBUG
1812 0 : CPLDebug("L1B", "Unknown spacecraft ID \"%d\".", iWord);
1813 : #endif
1814 0 : return CE_Failure;
1815 : }
1816 :
1817 : // Determine the product data type
1818 1 : iWord = GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_PROD_OFF);
1819 1 : switch (iWord)
1820 : {
1821 0 : case 1:
1822 0 : eProductType = LAC;
1823 0 : break;
1824 0 : case 2:
1825 0 : eProductType = GAC;
1826 0 : break;
1827 1 : case 3:
1828 1 : eProductType = HRPT;
1829 1 : break;
1830 0 : case 4: // XXX: documentation specifies the code '4'
1831 : case 13: // for FRAC but real datasets contain '13 here.'
1832 0 : eProductType = FRAC;
1833 0 : break;
1834 0 : default:
1835 : #ifdef DEBUG
1836 0 : CPLDebug("L1B", "Unknown product type \"%d\".", iWord);
1837 : #endif
1838 0 : return CE_Failure;
1839 : }
1840 :
1841 : // Fetch instrument status. Helps to determine whether we have
1842 : // 3A or 3B channel in the dataset.
1843 1 : iInstrumentStatus =
1844 1 : GetUInt32(abyRecHeader + L1B_NOAA15_HDR_REC_STAT_OFF);
1845 :
1846 : // Determine receiving station name
1847 1 : iWord = GetUInt16(abyRecHeader + L1B_NOAA15_HDR_REC_SRC_OFF);
1848 1 : switch (iWord)
1849 : {
1850 0 : case 1:
1851 0 : eSource = GC;
1852 0 : break;
1853 0 : case 2:
1854 0 : eSource = WI;
1855 0 : break;
1856 0 : case 3:
1857 0 : eSource = SO;
1858 0 : break;
1859 0 : case 4:
1860 0 : eSource = SV;
1861 0 : break;
1862 0 : case 5:
1863 0 : eSource = MO;
1864 0 : break;
1865 1 : default:
1866 1 : eSource = UNKNOWN_STATION;
1867 1 : break;
1868 1 : }
1869 : }
1870 : else
1871 0 : return CE_Failure;
1872 :
1873 : /* -------------------------------------------------------------------- */
1874 : /* Set fetched information as metadata records */
1875 : /* -------------------------------------------------------------------- */
1876 :
1877 1 : SetMetadataItem("DATASET_NAME", szDatasetName);
1878 :
1879 1 : const char *pszText = nullptr;
1880 1 : switch (eSpacecraftID)
1881 : {
1882 0 : case TIROSN:
1883 0 : pszText = "TIROS-N";
1884 0 : break;
1885 0 : case NOAA6:
1886 0 : pszText = "NOAA-6(A)";
1887 0 : break;
1888 0 : case NOAAB:
1889 0 : pszText = "NOAA-B";
1890 0 : break;
1891 0 : case NOAA7:
1892 0 : pszText = "NOAA-7(C)";
1893 0 : break;
1894 0 : case NOAA8:
1895 0 : pszText = "NOAA-8(E)";
1896 0 : break;
1897 0 : case NOAA9_UNKNOWN:
1898 0 : pszText = "UNKNOWN";
1899 0 : break;
1900 0 : case NOAA9:
1901 0 : pszText = "NOAA-9(F)";
1902 0 : break;
1903 0 : case NOAA10:
1904 0 : pszText = "NOAA-10(G)";
1905 0 : break;
1906 0 : case NOAA11:
1907 0 : pszText = "NOAA-11(H)";
1908 0 : break;
1909 0 : case NOAA12:
1910 0 : pszText = "NOAA-12(D)";
1911 0 : break;
1912 0 : case NOAA13:
1913 0 : pszText = "NOAA-13(I)";
1914 0 : break;
1915 0 : case NOAA14:
1916 0 : pszText = "NOAA-14(J)";
1917 0 : break;
1918 0 : case NOAA15:
1919 0 : pszText = "NOAA-15(K)";
1920 0 : break;
1921 0 : case NOAA16:
1922 0 : pszText = "NOAA-16(L)";
1923 0 : break;
1924 0 : case NOAA17:
1925 0 : pszText = "NOAA-17(M)";
1926 0 : break;
1927 0 : case NOAA18:
1928 0 : pszText = "NOAA-18(N)";
1929 0 : break;
1930 1 : case NOAA19:
1931 1 : pszText = "NOAA-19(N')";
1932 1 : break;
1933 0 : case METOP2:
1934 0 : pszText = "METOP-A(2)";
1935 0 : break;
1936 0 : case METOP1:
1937 0 : pszText = "METOP-B(1)";
1938 0 : break;
1939 0 : case METOP3:
1940 0 : pszText = "METOP-C(3)";
1941 0 : break;
1942 0 : default:
1943 0 : pszText = "Unknown";
1944 0 : break;
1945 : }
1946 1 : SetMetadataItem("SATELLITE", pszText);
1947 :
1948 1 : switch (eProductType)
1949 : {
1950 0 : case LAC:
1951 0 : pszText = "AVHRR LAC";
1952 0 : break;
1953 1 : case HRPT:
1954 1 : pszText = "AVHRR HRPT";
1955 1 : break;
1956 0 : case GAC:
1957 0 : pszText = "AVHRR GAC";
1958 0 : break;
1959 0 : case FRAC:
1960 0 : pszText = "AVHRR FRAC";
1961 0 : break;
1962 0 : default:
1963 0 : pszText = "Unknown";
1964 0 : break;
1965 : }
1966 1 : SetMetadataItem("DATA_TYPE", pszText);
1967 :
1968 : // Get revolution number as string, we don't need this value for processing
1969 : char szRevolution[6];
1970 1 : memcpy(szRevolution, szDatasetName + 32, 5);
1971 1 : szRevolution[5] = '\0';
1972 1 : SetMetadataItem("REVOLUTION", szRevolution);
1973 :
1974 1 : switch (eSource)
1975 : {
1976 0 : case DU:
1977 0 : pszText = "Dundee, Scotland, UK";
1978 0 : break;
1979 0 : case GC:
1980 0 : pszText = "Fairbanks, Alaska, USA (formerly Gilmore Creek)";
1981 0 : break;
1982 0 : case HO:
1983 0 : pszText = "Honolulu, Hawaii, USA";
1984 0 : break;
1985 0 : case MO:
1986 0 : pszText = "Monterey, California, USA";
1987 0 : break;
1988 0 : case WE:
1989 0 : pszText = "Western Europe CDA, Lannion, France";
1990 0 : break;
1991 0 : case SO:
1992 0 : pszText = "SOCC (Satellite Operations Control Center), Suitland, "
1993 : "Maryland, USA";
1994 0 : break;
1995 0 : case WI:
1996 0 : pszText = "Wallops Island, Virginia, USA";
1997 0 : break;
1998 1 : default:
1999 1 : pszText = "Unknown receiving station";
2000 1 : break;
2001 : }
2002 1 : SetMetadataItem("SOURCE", pszText);
2003 :
2004 1 : switch (eProcCenter)
2005 : {
2006 0 : case CMS:
2007 0 : pszText = "Centre de Meteorologie Spatiale - Lannion, France";
2008 0 : break;
2009 0 : case DSS:
2010 0 : pszText =
2011 : "Dundee Satellite Receiving Station - Dundee, Scotland, UK";
2012 0 : break;
2013 0 : case NSS:
2014 0 : pszText = "NOAA/NESDIS - Suitland, Maryland, USA";
2015 0 : break;
2016 0 : case UKM:
2017 0 : pszText =
2018 : "United Kingdom Meteorological Office - Bracknell, England, UK";
2019 0 : break;
2020 1 : default:
2021 1 : pszText = "Unknown processing center";
2022 1 : break;
2023 : }
2024 1 : SetMetadataItem("PROCESSING_CENTER", pszText);
2025 :
2026 1 : return CE_None;
2027 : }
2028 :
2029 : /************************************************************************/
2030 : /* ComputeFileOffsets() */
2031 : /************************************************************************/
2032 :
2033 1 : int L1BDataset::ComputeFileOffsets()
2034 : {
2035 1 : CPLDebug("L1B", "Data format = %s",
2036 1 : (iDataFormat == PACKED10BIT) ? "Packed 10 bit"
2037 1 : : (iDataFormat == UNPACKED16BIT) ? "Unpacked 16 bit"
2038 : : "Unpacked 8 bit");
2039 :
2040 1 : switch (eProductType)
2041 : {
2042 1 : case HRPT:
2043 : case LAC:
2044 : case FRAC:
2045 1 : nRasterXSize = 2048;
2046 1 : nBufferSize = 20484;
2047 1 : iGCPStart =
2048 : 25 -
2049 : 1; /* http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c2/sec2-4.htm
2050 : */
2051 1 : iGCPStep = 40;
2052 1 : nGCPsPerLine = 51;
2053 1 : if (eL1BFormat == L1B_NOAA9)
2054 : {
2055 0 : if (iDataFormat == PACKED10BIT)
2056 : {
2057 0 : nRecordSize = 14800;
2058 0 : nRecordDataEnd = 14104;
2059 : }
2060 0 : else if (iDataFormat == UNPACKED16BIT)
2061 : {
2062 0 : switch (nBands)
2063 : {
2064 0 : case 1:
2065 0 : nRecordSize = 4544;
2066 0 : nRecordDataEnd = 4544;
2067 0 : break;
2068 0 : case 2:
2069 0 : nRecordSize = 8640;
2070 0 : nRecordDataEnd = 8640;
2071 0 : break;
2072 0 : case 3:
2073 0 : nRecordSize = 12736;
2074 0 : nRecordDataEnd = 12736;
2075 0 : break;
2076 0 : case 4:
2077 0 : nRecordSize = 16832;
2078 0 : nRecordDataEnd = 16832;
2079 0 : break;
2080 0 : case 5:
2081 0 : nRecordSize = 20928;
2082 0 : nRecordDataEnd = 20928;
2083 0 : break;
2084 : }
2085 : }
2086 : else // UNPACKED8BIT
2087 : {
2088 0 : switch (nBands)
2089 : {
2090 0 : case 1:
2091 0 : nRecordSize = 2496;
2092 0 : nRecordDataEnd = 2496;
2093 0 : break;
2094 0 : case 2:
2095 0 : nRecordSize = 4544;
2096 0 : nRecordDataEnd = 4544;
2097 0 : break;
2098 0 : case 3:
2099 0 : nRecordSize = 6592;
2100 0 : nRecordDataEnd = 6592;
2101 0 : break;
2102 0 : case 4:
2103 0 : nRecordSize = 8640;
2104 0 : nRecordDataEnd = 8640;
2105 0 : break;
2106 0 : case 5:
2107 0 : nRecordSize = 10688;
2108 0 : nRecordDataEnd = 10688;
2109 0 : break;
2110 : }
2111 : }
2112 0 : nDataStartOffset = nRecordSize + L1B_NOAA9_HEADER_SIZE;
2113 0 : nRecordDataStart = 448;
2114 0 : iGCPCodeOffset = 52;
2115 0 : iGCPOffset = 104;
2116 : }
2117 :
2118 1 : else if (eL1BFormat == L1B_NOAA15 || eL1BFormat == L1B_NOAA15_NOHDR)
2119 : {
2120 1 : if (iDataFormat == PACKED10BIT)
2121 : {
2122 0 : nRecordSize = 15872;
2123 0 : nRecordDataEnd = 14920;
2124 0 : iCLAVRStart = 14984;
2125 : }
2126 1 : else if (iDataFormat == UNPACKED16BIT)
2127 : { /* Table 8.3.1.3.3.1-3 */
2128 1 : switch (nBands)
2129 : {
2130 0 : case 1:
2131 0 : nRecordSize = 6144;
2132 0 : nRecordDataEnd = 5360;
2133 0 : iCLAVRStart =
2134 : 5368 + 56; /* guessed but not verified */
2135 0 : break;
2136 0 : case 2:
2137 0 : nRecordSize = 10240;
2138 0 : nRecordDataEnd = 9456;
2139 0 : iCLAVRStart =
2140 : 9464 + 56; /* guessed but not verified */
2141 0 : break;
2142 0 : case 3:
2143 0 : nRecordSize = 14336;
2144 0 : nRecordDataEnd = 13552;
2145 0 : iCLAVRStart =
2146 : 13560 + 56; /* guessed but not verified */
2147 0 : break;
2148 0 : case 4:
2149 0 : nRecordSize = 18432;
2150 0 : nRecordDataEnd = 17648;
2151 0 : iCLAVRStart =
2152 : 17656 + 56; /* guessed but not verified */
2153 0 : break;
2154 1 : case 5:
2155 1 : nRecordSize = 22528;
2156 1 : nRecordDataEnd = 21744;
2157 1 : iCLAVRStart = 21752 + 56;
2158 1 : break;
2159 : }
2160 : }
2161 : else // UNPACKED8BIT
2162 : { /* Table 8.3.1.3.3.1-2 */
2163 0 : switch (nBands)
2164 : {
2165 0 : case 1:
2166 0 : nRecordSize = 4096;
2167 0 : nRecordDataEnd = 3312;
2168 0 : iCLAVRStart =
2169 : 3320 + 56; /* guessed but not verified */
2170 0 : break;
2171 0 : case 2:
2172 0 : nRecordSize = 6144;
2173 0 : nRecordDataEnd = 5360;
2174 0 : iCLAVRStart =
2175 : 5368 + 56; /* guessed but not verified */
2176 0 : break;
2177 0 : case 3:
2178 0 : nRecordSize = 8192;
2179 0 : nRecordDataEnd = 7408;
2180 0 : iCLAVRStart =
2181 : 7416 + 56; /* guessed but not verified */
2182 0 : break;
2183 0 : case 4:
2184 0 : nRecordSize = 10240;
2185 0 : nRecordDataEnd = 9456;
2186 0 : iCLAVRStart =
2187 : 9464 + 56; /* guessed but not verified */
2188 0 : break;
2189 0 : case 5:
2190 0 : nRecordSize = 12288;
2191 0 : nRecordDataEnd = 11504;
2192 0 : iCLAVRStart =
2193 : 11512 + 56; /* guessed but not verified */
2194 0 : break;
2195 : }
2196 : }
2197 2 : nDataStartOffset = (eL1BFormat == L1B_NOAA15_NOHDR)
2198 1 : ? nRecordDataEnd
2199 0 : : nRecordSize + L1B_NOAA15_HEADER_SIZE;
2200 1 : nRecordDataStart = 1264;
2201 1 : iGCPCodeOffset = 0; // XXX: not exist for NOAA15?
2202 1 : iGCPOffset = 640;
2203 : }
2204 : else
2205 0 : return 0;
2206 1 : break;
2207 :
2208 0 : case GAC:
2209 0 : nRasterXSize = 409;
2210 0 : nBufferSize = 4092;
2211 0 : iGCPStart = 5 - 1; // FIXME: depends of scan direction
2212 0 : iGCPStep = 8;
2213 0 : nGCPsPerLine = 51;
2214 0 : if (eL1BFormat == L1B_NOAA9)
2215 : {
2216 0 : if (iDataFormat == PACKED10BIT)
2217 : {
2218 0 : nRecordSize = 3220;
2219 0 : nRecordDataEnd = 3176;
2220 : }
2221 0 : else if (iDataFormat == UNPACKED16BIT)
2222 0 : switch (nBands)
2223 : {
2224 0 : case 1:
2225 0 : nRecordSize = 1268;
2226 0 : nRecordDataEnd = 1266;
2227 0 : break;
2228 0 : case 2:
2229 0 : nRecordSize = 2084;
2230 0 : nRecordDataEnd = 2084;
2231 0 : break;
2232 0 : case 3:
2233 0 : nRecordSize = 2904;
2234 0 : nRecordDataEnd = 2902;
2235 0 : break;
2236 0 : case 4:
2237 0 : nRecordSize = 3720;
2238 0 : nRecordDataEnd = 3720;
2239 0 : break;
2240 0 : case 5:
2241 0 : nRecordSize = 4540;
2242 0 : nRecordDataEnd = 4538;
2243 0 : break;
2244 : }
2245 : else // UNPACKED8BIT
2246 : {
2247 0 : switch (nBands)
2248 : {
2249 0 : case 1:
2250 0 : nRecordSize = 860;
2251 0 : nRecordDataEnd = 858;
2252 0 : break;
2253 0 : case 2:
2254 0 : nRecordSize = 1268;
2255 0 : nRecordDataEnd = 1266;
2256 0 : break;
2257 0 : case 3:
2258 0 : nRecordSize = 1676;
2259 0 : nRecordDataEnd = 1676;
2260 0 : break;
2261 0 : case 4:
2262 0 : nRecordSize = 2084;
2263 0 : nRecordDataEnd = 2084;
2264 0 : break;
2265 0 : case 5:
2266 0 : nRecordSize = 2496;
2267 0 : nRecordDataEnd = 2494;
2268 0 : break;
2269 : }
2270 : }
2271 0 : nDataStartOffset = nRecordSize * 2 + L1B_NOAA9_HEADER_SIZE;
2272 0 : nRecordDataStart = 448;
2273 0 : iGCPCodeOffset = 52;
2274 0 : iGCPOffset = 104;
2275 : }
2276 :
2277 0 : else if (eL1BFormat == L1B_NOAA15 || eL1BFormat == L1B_NOAA15_NOHDR)
2278 : {
2279 0 : if (iDataFormat == PACKED10BIT)
2280 : {
2281 0 : nRecordSize = 4608;
2282 0 : nRecordDataEnd = 3992;
2283 0 : iCLAVRStart = 4056;
2284 : }
2285 0 : else if (iDataFormat == UNPACKED16BIT)
2286 : { /* Table 8.3.1.4.3.1-3 */
2287 0 : switch (nBands)
2288 : {
2289 0 : case 1:
2290 0 : nRecordSize = 2360;
2291 0 : nRecordDataEnd = 2082;
2292 0 : iCLAVRStart =
2293 : 2088 + 56; /* guessed but not verified */
2294 0 : break;
2295 0 : case 2:
2296 0 : nRecordSize = 3176;
2297 0 : nRecordDataEnd = 2900;
2298 0 : iCLAVRStart =
2299 : 2904 + 56; /* guessed but not verified */
2300 0 : break;
2301 0 : case 3:
2302 0 : nRecordSize = 3992;
2303 0 : nRecordDataEnd = 3718;
2304 0 : iCLAVRStart =
2305 : 3720 + 56; /* guessed but not verified */
2306 0 : break;
2307 0 : case 4:
2308 0 : nRecordSize = 4816;
2309 0 : nRecordDataEnd = 4536;
2310 0 : iCLAVRStart =
2311 : 4544 + 56; /* guessed but not verified */
2312 0 : break;
2313 0 : case 5:
2314 0 : nRecordSize = 5632;
2315 0 : nRecordDataEnd = 5354;
2316 0 : iCLAVRStart = 5360 + 56;
2317 0 : break;
2318 : }
2319 : }
2320 : else // UNPACKED8BIT
2321 : { /* Table 8.3.1.4.3.1-2 but record length is wrong in the table
2322 : ! */
2323 0 : switch (nBands)
2324 : {
2325 0 : case 1:
2326 0 : nRecordSize = 1952;
2327 0 : nRecordDataEnd = 1673;
2328 0 : iCLAVRStart =
2329 : 1680 + 56; /* guessed but not verified */
2330 0 : break;
2331 0 : case 2:
2332 0 : nRecordSize = 2360;
2333 0 : nRecordDataEnd = 2082;
2334 0 : iCLAVRStart =
2335 : 2088 + 56; /* guessed but not verified */
2336 0 : break;
2337 0 : case 3:
2338 0 : nRecordSize = 2768;
2339 0 : nRecordDataEnd = 2491;
2340 0 : iCLAVRStart =
2341 : 2496 + 56; /* guessed but not verified */
2342 0 : break;
2343 0 : case 4:
2344 0 : nRecordSize = 3176;
2345 0 : nRecordDataEnd = 2900;
2346 0 : iCLAVRStart =
2347 : 2904 + 56; /* guessed but not verified */
2348 0 : break;
2349 0 : case 5:
2350 0 : nRecordSize = 3584;
2351 0 : nRecordDataEnd = 3309;
2352 0 : iCLAVRStart =
2353 : 3312 + 56; /* guessed but not verified */
2354 0 : break;
2355 : }
2356 : }
2357 0 : nDataStartOffset = (eL1BFormat == L1B_NOAA15_NOHDR)
2358 0 : ? nRecordDataEnd
2359 0 : : nRecordSize + L1B_NOAA15_HEADER_SIZE;
2360 0 : nRecordDataStart = 1264;
2361 0 : iGCPCodeOffset = 0; // XXX: not exist for NOAA15?
2362 0 : iGCPOffset = 640;
2363 : }
2364 : else
2365 0 : return 0;
2366 0 : break;
2367 0 : default:
2368 0 : return 0;
2369 : }
2370 :
2371 1 : return 1;
2372 : }
2373 :
2374 : /************************************************************************/
2375 : /* L1BGeolocDataset */
2376 : /************************************************************************/
2377 :
2378 : class L1BGeolocDataset final : public GDALDataset
2379 : {
2380 : friend class L1BGeolocRasterBand;
2381 :
2382 : L1BDataset *poL1BDS;
2383 : int bInterpolGeolocationDS;
2384 :
2385 : public:
2386 : L1BGeolocDataset(L1BDataset *poMainDS, int bInterpolGeolocationDS);
2387 : virtual ~L1BGeolocDataset();
2388 :
2389 : static GDALDataset *CreateGeolocationDS(L1BDataset *poL1BDS,
2390 : int bInterpolGeolocationDS);
2391 : };
2392 :
2393 : /************************************************************************/
2394 : /* L1BGeolocRasterBand */
2395 : /************************************************************************/
2396 :
2397 : class L1BGeolocRasterBand final : public GDALRasterBand
2398 : {
2399 : public:
2400 : L1BGeolocRasterBand(L1BGeolocDataset *poDS, int nBand);
2401 :
2402 : virtual CPLErr IReadBlock(int, int, void *) override;
2403 : virtual double GetNoDataValue(int *pbSuccess = nullptr) override;
2404 : };
2405 :
2406 : /************************************************************************/
2407 : /* L1BGeolocDataset() */
2408 : /************************************************************************/
2409 :
2410 0 : L1BGeolocDataset::L1BGeolocDataset(L1BDataset *poL1BDSIn,
2411 0 : int bInterpolGeolocationDSIn)
2412 0 : : poL1BDS(poL1BDSIn), bInterpolGeolocationDS(bInterpolGeolocationDSIn)
2413 : {
2414 0 : if (bInterpolGeolocationDS)
2415 0 : nRasterXSize = poL1BDS->nRasterXSize;
2416 : else
2417 0 : nRasterXSize = poL1BDS->nGCPsPerLine;
2418 0 : nRasterYSize = poL1BDS->nRasterYSize;
2419 0 : }
2420 :
2421 : /************************************************************************/
2422 : /* ~L1BGeolocDataset() */
2423 : /************************************************************************/
2424 :
2425 0 : L1BGeolocDataset::~L1BGeolocDataset()
2426 : {
2427 0 : delete poL1BDS;
2428 0 : }
2429 :
2430 : /************************************************************************/
2431 : /* L1BGeolocRasterBand() */
2432 : /************************************************************************/
2433 :
2434 0 : L1BGeolocRasterBand::L1BGeolocRasterBand(L1BGeolocDataset *poDSIn, int nBandIn)
2435 : {
2436 0 : poDS = poDSIn;
2437 0 : nBand = nBandIn;
2438 0 : nRasterXSize = poDSIn->nRasterXSize;
2439 0 : nRasterYSize = poDSIn->nRasterYSize;
2440 0 : eDataType = GDT_Float64;
2441 0 : nBlockXSize = nRasterXSize;
2442 0 : nBlockYSize = 1;
2443 0 : if (nBand == 1)
2444 0 : SetDescription("GEOLOC X");
2445 : else
2446 0 : SetDescription("GEOLOC Y");
2447 0 : }
2448 :
2449 : /************************************************************************/
2450 : /* LagrangeInterpol() */
2451 : /************************************************************************/
2452 :
2453 : /* ----------------------------------------------------------------------------
2454 : * Perform a Lagrangian interpolation through the given x,y coordinates
2455 : * and return the interpolated y value for the given x value.
2456 : * The array size and thus the polynomial order is defined by numpt.
2457 : * Input: x[] and y[] are of size numpt,
2458 : * x0 is the x value for which we calculate the corresponding y
2459 : * Returns: y value calculated for given x0.
2460 : */
2461 0 : static double LagrangeInterpol(const double x[], const double y[], double x0,
2462 : int numpt)
2463 : {
2464 : int i, j;
2465 : double L;
2466 0 : double y0 = 0;
2467 :
2468 0 : for (i = 0; i < numpt; i++)
2469 : {
2470 0 : L = 1.0;
2471 0 : for (j = 0; j < numpt; j++)
2472 : {
2473 0 : if (i == j)
2474 0 : continue;
2475 0 : L = L * (x0 - x[j]) / (x[i] - x[j]);
2476 : }
2477 0 : y0 = y0 + L * y[i];
2478 : }
2479 0 : return y0;
2480 : }
2481 :
2482 : /************************************************************************/
2483 : /* L1BInterpol() */
2484 : /************************************************************************/
2485 :
2486 : /* ----------------------------------------------------------------------------
2487 : * Interpolate an array of size numPoints where the only values set on input are
2488 : * at knownFirst, and intervals of knownStep thereafter.
2489 : * On return all the rest from 0..numPoints-1 will be filled in.
2490 : * Uses the LagrangeInterpol() function to do the interpolation; 5-point for the
2491 : * beginning and end of the array and 4-point for the rest.
2492 : * To use this function for NOAA level 1B data extract the 51 latitude values
2493 : * into their appropriate places in the vals array then call L1BInterpol to
2494 : * calculate the rest of the values. Do similarly for longitudes, solar zenith
2495 : * angles, and any others which are present in the file.
2496 : * Reference:
2497 : * http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/html/c2/sec2-4.htm
2498 : */
2499 :
2500 : #define MIDDLE_INTERP_ORDER 4
2501 : #define END_INTERP_ORDER 5 /* Ensure this is an odd number, 5 is suitable.*/
2502 :
2503 : /* Convert number of known point to its index in full array */
2504 : #define IDX(N) ((N)*knownStep + knownFirst)
2505 :
2506 0 : static void L1BInterpol(
2507 : double vals[], int numKnown, /* Number of known points (typically 51) */
2508 : int knownFirst, /* Index in full array of first known point (24) */
2509 : int knownStep, /* Interval to next and subsequent known points (40) */
2510 : int numPoints /* Number of points in whole array (2048) */)
2511 : {
2512 : int i, j;
2513 : double x[END_INTERP_ORDER];
2514 : double y[END_INTERP_ORDER];
2515 :
2516 : /* First extrapolate first 24 points */
2517 0 : for (i = 0; i < END_INTERP_ORDER; i++)
2518 : {
2519 0 : x[i] = IDX(i);
2520 0 : y[i] = vals[IDX(i)];
2521 : }
2522 0 : for (i = 0; i < knownFirst; i++)
2523 : {
2524 0 : vals[i] = LagrangeInterpol(x, y, i, END_INTERP_ORDER);
2525 : }
2526 :
2527 : /* Next extrapolate last 23 points */
2528 0 : for (i = 0; i < END_INTERP_ORDER; i++)
2529 : {
2530 0 : x[i] = IDX(numKnown - END_INTERP_ORDER + i);
2531 0 : y[i] = vals[IDX(numKnown - END_INTERP_ORDER + i)];
2532 : }
2533 0 : for (i = IDX(numKnown - 1); i < numPoints; i++)
2534 : {
2535 0 : vals[i] = LagrangeInterpol(x, y, i, END_INTERP_ORDER);
2536 : }
2537 :
2538 : /* Interpolate all intermediate points using two before and two after */
2539 0 : for (i = knownFirst; i < IDX(numKnown - 1); i++)
2540 : {
2541 : double x2[MIDDLE_INTERP_ORDER];
2542 : double y2[MIDDLE_INTERP_ORDER];
2543 : int startpt;
2544 :
2545 : /* Find a suitable set of two known points before and two after */
2546 0 : startpt = (i / knownStep) - MIDDLE_INTERP_ORDER / 2;
2547 0 : if (startpt < 0)
2548 0 : startpt = 0;
2549 0 : if (startpt + MIDDLE_INTERP_ORDER - 1 >= numKnown)
2550 0 : startpt = numKnown - MIDDLE_INTERP_ORDER;
2551 0 : for (j = 0; j < MIDDLE_INTERP_ORDER; j++)
2552 : {
2553 0 : x2[j] = IDX(startpt + j);
2554 0 : y2[j] = vals[IDX(startpt + j)];
2555 : }
2556 0 : vals[i] = LagrangeInterpol(x2, y2, i, MIDDLE_INTERP_ORDER);
2557 : }
2558 0 : }
2559 :
2560 : /************************************************************************/
2561 : /* IReadBlock() */
2562 : /************************************************************************/
2563 :
2564 0 : CPLErr L1BGeolocRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
2565 : int nBlockYOff, void *pData)
2566 : {
2567 0 : L1BGeolocDataset *poGDS = (L1BGeolocDataset *)poDS;
2568 0 : L1BDataset *poL1BDS = poGDS->poL1BDS;
2569 : GDAL_GCP *pasGCPList =
2570 0 : (GDAL_GCP *)CPLCalloc(poL1BDS->nGCPsPerLine, sizeof(GDAL_GCP));
2571 0 : GDALInitGCPs(poL1BDS->nGCPsPerLine, pasGCPList);
2572 :
2573 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);
2574 :
2575 : /* -------------------------------------------------------------------- */
2576 : /* Seek to data. */
2577 : /* -------------------------------------------------------------------- */
2578 0 : CPL_IGNORE_RET_VAL(
2579 0 : VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));
2580 :
2581 0 : CPL_IGNORE_RET_VAL(
2582 0 : VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordDataStart, poL1BDS->fp));
2583 :
2584 : /* Fetch the GCPs for the row */
2585 0 : int nGotGCPs = poL1BDS->FetchGCPs(pasGCPList, pabyRecordHeader, nBlockYOff);
2586 0 : double *padfData = (double *)pData;
2587 : int i;
2588 0 : if (poGDS->bInterpolGeolocationDS)
2589 : {
2590 : /* Fill the known position */
2591 0 : for (i = 0; i < nGotGCPs; i++)
2592 : {
2593 0 : double dfVal =
2594 0 : (nBand == 1) ? pasGCPList[i].dfGCPX : pasGCPList[i].dfGCPY;
2595 0 : padfData[poL1BDS->iGCPStart + i * poL1BDS->iGCPStep] = dfVal;
2596 : }
2597 :
2598 0 : if (nGotGCPs == poL1BDS->nGCPsPerLine)
2599 : {
2600 : /* And do Lagangian interpolation to fill the holes */
2601 0 : L1BInterpol(padfData, poL1BDS->nGCPsPerLine, poL1BDS->iGCPStart,
2602 : poL1BDS->iGCPStep, nRasterXSize);
2603 : }
2604 : else
2605 : {
2606 0 : int iFirstNonValid = 0;
2607 0 : if (nGotGCPs > 5)
2608 0 : iFirstNonValid = poL1BDS->iGCPStart +
2609 0 : nGotGCPs * poL1BDS->iGCPStep +
2610 0 : poL1BDS->iGCPStep / 2;
2611 0 : for (i = iFirstNonValid; i < nRasterXSize; i++)
2612 : {
2613 0 : padfData[i] = GetNoDataValue(nullptr);
2614 : }
2615 0 : if (iFirstNonValid > 0)
2616 : {
2617 0 : L1BInterpol(padfData, poL1BDS->nGCPsPerLine, poL1BDS->iGCPStart,
2618 : poL1BDS->iGCPStep, iFirstNonValid);
2619 : }
2620 : }
2621 : }
2622 : else
2623 : {
2624 0 : for (i = 0; i < nGotGCPs; i++)
2625 : {
2626 0 : padfData[i] =
2627 0 : (nBand == 1) ? pasGCPList[i].dfGCPX : pasGCPList[i].dfGCPY;
2628 : }
2629 0 : for (i = nGotGCPs; i < nRasterXSize; i++)
2630 0 : padfData[i] = GetNoDataValue(nullptr);
2631 : }
2632 :
2633 0 : if (poL1BDS->eLocationIndicator == ASCEND)
2634 : {
2635 0 : for (i = 0; i < nRasterXSize / 2; i++)
2636 : {
2637 0 : double dfTmp = padfData[i];
2638 0 : padfData[i] = padfData[nRasterXSize - 1 - i];
2639 0 : padfData[nRasterXSize - 1 - i] = dfTmp;
2640 : }
2641 : }
2642 :
2643 0 : CPLFree(pabyRecordHeader);
2644 0 : GDALDeinitGCPs(poL1BDS->nGCPsPerLine, pasGCPList);
2645 0 : CPLFree(pasGCPList);
2646 :
2647 0 : return CE_None;
2648 : }
2649 :
2650 : /************************************************************************/
2651 : /* GetNoDataValue() */
2652 : /************************************************************************/
2653 :
2654 0 : double L1BGeolocRasterBand::GetNoDataValue(int *pbSuccess)
2655 : {
2656 0 : if (pbSuccess)
2657 0 : *pbSuccess = TRUE;
2658 0 : return -200.0;
2659 : }
2660 :
2661 : /************************************************************************/
2662 : /* CreateGeolocationDS() */
2663 : /************************************************************************/
2664 :
2665 0 : GDALDataset *L1BGeolocDataset::CreateGeolocationDS(L1BDataset *poL1BDS,
2666 : int bInterpolGeolocationDS)
2667 : {
2668 : L1BGeolocDataset *poGeolocDS =
2669 0 : new L1BGeolocDataset(poL1BDS, bInterpolGeolocationDS);
2670 0 : for (int i = 1; i <= 2; i++)
2671 : {
2672 0 : poGeolocDS->SetBand(i, new L1BGeolocRasterBand(poGeolocDS, i));
2673 : }
2674 0 : return poGeolocDS;
2675 : }
2676 :
2677 : /************************************************************************/
2678 : /* L1BSolarZenithAnglesDataset */
2679 : /************************************************************************/
2680 :
2681 : class L1BSolarZenithAnglesDataset final : public GDALDataset
2682 : {
2683 : friend class L1BSolarZenithAnglesRasterBand;
2684 :
2685 : L1BDataset *poL1BDS;
2686 :
2687 : public:
2688 : explicit L1BSolarZenithAnglesDataset(L1BDataset *poMainDS);
2689 : virtual ~L1BSolarZenithAnglesDataset();
2690 :
2691 : static GDALDataset *CreateSolarZenithAnglesDS(L1BDataset *poL1BDS);
2692 : };
2693 :
2694 : /************************************************************************/
2695 : /* L1BSolarZenithAnglesRasterBand */
2696 : /************************************************************************/
2697 :
2698 : class L1BSolarZenithAnglesRasterBand final : public GDALRasterBand
2699 : {
2700 : public:
2701 : L1BSolarZenithAnglesRasterBand(L1BSolarZenithAnglesDataset *poDS,
2702 : int nBand);
2703 :
2704 : virtual CPLErr IReadBlock(int, int, void *) override;
2705 : virtual double GetNoDataValue(int *pbSuccess = nullptr) override;
2706 : };
2707 :
2708 : /************************************************************************/
2709 : /* L1BSolarZenithAnglesDataset() */
2710 : /************************************************************************/
2711 :
2712 0 : L1BSolarZenithAnglesDataset::L1BSolarZenithAnglesDataset(L1BDataset *poL1BDSIn)
2713 : {
2714 0 : poL1BDS = poL1BDSIn;
2715 0 : nRasterXSize = 51;
2716 0 : nRasterYSize = poL1BDSIn->nRasterYSize;
2717 0 : }
2718 :
2719 : /************************************************************************/
2720 : /* ~L1BSolarZenithAnglesDataset() */
2721 : /************************************************************************/
2722 :
2723 0 : L1BSolarZenithAnglesDataset::~L1BSolarZenithAnglesDataset()
2724 : {
2725 0 : delete poL1BDS;
2726 0 : }
2727 :
2728 : /************************************************************************/
2729 : /* L1BSolarZenithAnglesRasterBand() */
2730 : /************************************************************************/
2731 :
2732 0 : L1BSolarZenithAnglesRasterBand::L1BSolarZenithAnglesRasterBand(
2733 0 : L1BSolarZenithAnglesDataset *poDSIn, int nBandIn)
2734 : {
2735 0 : poDS = poDSIn;
2736 0 : nBand = nBandIn;
2737 0 : nRasterXSize = poDSIn->nRasterXSize;
2738 0 : nRasterYSize = poDSIn->nRasterYSize;
2739 0 : eDataType = GDT_Float32;
2740 0 : nBlockXSize = nRasterXSize;
2741 0 : nBlockYSize = 1;
2742 0 : }
2743 :
2744 : /************************************************************************/
2745 : /* IReadBlock() */
2746 : /************************************************************************/
2747 :
2748 0 : CPLErr L1BSolarZenithAnglesRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
2749 : int nBlockYOff, void *pData)
2750 : {
2751 0 : L1BSolarZenithAnglesDataset *poGDS = (L1BSolarZenithAnglesDataset *)poDS;
2752 0 : L1BDataset *poL1BDS = poGDS->poL1BDS;
2753 : int i;
2754 :
2755 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);
2756 :
2757 : /* -------------------------------------------------------------------- */
2758 : /* Seek to data. */
2759 : /* -------------------------------------------------------------------- */
2760 0 : CPL_IGNORE_RET_VAL(
2761 0 : VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));
2762 :
2763 0 : CPL_IGNORE_RET_VAL(
2764 0 : VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordSize, poL1BDS->fp));
2765 :
2766 : const int nValidValues =
2767 0 : std::min(nRasterXSize,
2768 0 : static_cast<int>(pabyRecordHeader[poL1BDS->iGCPCodeOffset]));
2769 0 : float *pafData = (float *)pData;
2770 :
2771 0 : int bHasFractional = (poL1BDS->nRecordDataEnd + 20 <= poL1BDS->nRecordSize);
2772 :
2773 : #ifdef notdef
2774 : if (bHasFractional)
2775 : {
2776 : for (i = 0; i < 20; i++)
2777 : {
2778 : GByte val = pabyRecordHeader[poL1BDS->nRecordDataEnd + i];
2779 : for (int j = 0; j < 8; j++)
2780 : fprintf(stderr, "%c", /*ok*/
2781 : ((val >> (7 - j)) & 1) ? '1' : '0');
2782 : fprintf(stderr, " "); /*ok*/
2783 : }
2784 : fprintf(stderr, "\n"); /*ok*/
2785 : }
2786 : #endif
2787 :
2788 0 : for (i = 0; i < nValidValues; i++)
2789 : {
2790 0 : pafData[i] = pabyRecordHeader[poL1BDS->iGCPCodeOffset + 1 + i] / 2.0f;
2791 :
2792 0 : if (bHasFractional)
2793 : {
2794 : /* Cf
2795 : * http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/l/app-l.htm#notl-2
2796 : */
2797 : /* This is not very clear on if bits must be counted from MSB or LSB
2798 : */
2799 : /* but when testing on n12gac10bit.l1b, it appears that the 3 bits
2800 : * for i=0 are the 3 MSB bits */
2801 0 : int nAddBitStart = i * 3;
2802 : int nFractional;
2803 :
2804 : #if 1
2805 0 : if ((nAddBitStart % 8) + 3 <= 8)
2806 : {
2807 0 : nFractional = (pabyRecordHeader[poL1BDS->nRecordDataEnd +
2808 0 : (nAddBitStart / 8)] >>
2809 0 : (8 - ((nAddBitStart % 8) + 3))) &
2810 : 0x7;
2811 : }
2812 : else
2813 : {
2814 0 : nFractional = (((pabyRecordHeader[poL1BDS->nRecordDataEnd +
2815 0 : (nAddBitStart / 8)]
2816 0 : << 8) |
2817 0 : pabyRecordHeader[poL1BDS->nRecordDataEnd +
2818 0 : (nAddBitStart / 8) + 1]) >>
2819 0 : (16 - ((nAddBitStart % 8) + 3))) &
2820 : 0x7;
2821 : }
2822 : #else
2823 : nFractional = (pabyRecordHeader[poL1BDS->nRecordDataEnd +
2824 : (nAddBitStart / 8)] >>
2825 : (nAddBitStart % 8)) &
2826 : 0x7;
2827 : if ((nAddBitStart % 8) + 3 > 8)
2828 : nFractional |= (pabyRecordHeader[poL1BDS->nRecordDataEnd +
2829 : (nAddBitStart / 8) + 1] &
2830 : ((1 << (((nAddBitStart % 8) + 3 - 8))) - 1))
2831 : << (3 - ((((nAddBitStart % 8) + 3 - 8))));
2832 : * /
2833 : #endif
2834 0 : if (nFractional > 4)
2835 : {
2836 0 : CPLDebug("L1B",
2837 : "For nBlockYOff=%d, i=%d, wrong fractional value : %d",
2838 : nBlockYOff, i, nFractional);
2839 : }
2840 :
2841 0 : pafData[i] += nFractional / 10.0f;
2842 : }
2843 : }
2844 :
2845 0 : for (; i < nRasterXSize; i++)
2846 0 : pafData[i] = static_cast<float>(GetNoDataValue(nullptr));
2847 :
2848 0 : if (poL1BDS->eLocationIndicator == ASCEND)
2849 : {
2850 0 : for (i = 0; i < nRasterXSize / 2; i++)
2851 : {
2852 0 : float fTmp = pafData[i];
2853 0 : pafData[i] = pafData[nRasterXSize - 1 - i];
2854 0 : pafData[nRasterXSize - 1 - i] = fTmp;
2855 : }
2856 : }
2857 :
2858 0 : CPLFree(pabyRecordHeader);
2859 :
2860 0 : return CE_None;
2861 : }
2862 :
2863 : /************************************************************************/
2864 : /* GetNoDataValue() */
2865 : /************************************************************************/
2866 :
2867 0 : double L1BSolarZenithAnglesRasterBand::GetNoDataValue(int *pbSuccess)
2868 : {
2869 0 : if (pbSuccess)
2870 0 : *pbSuccess = TRUE;
2871 0 : return -200.0;
2872 : }
2873 :
2874 : /************************************************************************/
2875 : /* CreateSolarZenithAnglesDS() */
2876 : /************************************************************************/
2877 :
2878 : GDALDataset *
2879 0 : L1BSolarZenithAnglesDataset::CreateSolarZenithAnglesDS(L1BDataset *poL1BDS)
2880 : {
2881 : L1BSolarZenithAnglesDataset *poGeolocDS =
2882 0 : new L1BSolarZenithAnglesDataset(poL1BDS);
2883 0 : for (int i = 1; i <= 1; i++)
2884 : {
2885 0 : poGeolocDS->SetBand(i,
2886 0 : new L1BSolarZenithAnglesRasterBand(poGeolocDS, i));
2887 : }
2888 0 : return poGeolocDS;
2889 : }
2890 :
2891 : /************************************************************************/
2892 : /* L1BNOAA15AnglesDataset */
2893 : /************************************************************************/
2894 :
2895 : class L1BNOAA15AnglesDataset final : public GDALDataset
2896 : {
2897 : friend class L1BNOAA15AnglesRasterBand;
2898 :
2899 : L1BDataset *poL1BDS;
2900 :
2901 : public:
2902 : explicit L1BNOAA15AnglesDataset(L1BDataset *poMainDS);
2903 : virtual ~L1BNOAA15AnglesDataset();
2904 :
2905 : static GDALDataset *CreateAnglesDS(L1BDataset *poL1BDS);
2906 : };
2907 :
2908 : /************************************************************************/
2909 : /* L1BNOAA15AnglesRasterBand */
2910 : /************************************************************************/
2911 :
2912 : class L1BNOAA15AnglesRasterBand final : public GDALRasterBand
2913 : {
2914 : public:
2915 : L1BNOAA15AnglesRasterBand(L1BNOAA15AnglesDataset *poDS, int nBand);
2916 :
2917 : virtual CPLErr IReadBlock(int, int, void *) override;
2918 : };
2919 :
2920 : /************************************************************************/
2921 : /* L1BNOAA15AnglesDataset() */
2922 : /************************************************************************/
2923 :
2924 0 : L1BNOAA15AnglesDataset::L1BNOAA15AnglesDataset(L1BDataset *poL1BDSIn)
2925 0 : : poL1BDS(poL1BDSIn)
2926 : {
2927 0 : nRasterXSize = 51;
2928 0 : nRasterYSize = poL1BDS->nRasterYSize;
2929 0 : }
2930 :
2931 : /************************************************************************/
2932 : /* ~L1BNOAA15AnglesDataset() */
2933 : /************************************************************************/
2934 :
2935 0 : L1BNOAA15AnglesDataset::~L1BNOAA15AnglesDataset()
2936 : {
2937 0 : delete poL1BDS;
2938 0 : }
2939 :
2940 : /************************************************************************/
2941 : /* L1BNOAA15AnglesRasterBand() */
2942 : /************************************************************************/
2943 :
2944 0 : L1BNOAA15AnglesRasterBand::L1BNOAA15AnglesRasterBand(
2945 0 : L1BNOAA15AnglesDataset *poDSIn, int nBandIn)
2946 : {
2947 0 : poDS = poDSIn;
2948 0 : nBand = nBandIn;
2949 0 : nRasterXSize = poDSIn->nRasterXSize;
2950 0 : nRasterYSize = poDSIn->nRasterYSize;
2951 0 : eDataType = GDT_Float32;
2952 0 : nBlockXSize = nRasterXSize;
2953 0 : nBlockYSize = 1;
2954 0 : if (nBand == 1)
2955 0 : SetDescription("Solar zenith angles");
2956 0 : else if (nBand == 2)
2957 0 : SetDescription("Satellite zenith angles");
2958 : else
2959 0 : SetDescription("Relative azimuth angles");
2960 0 : }
2961 :
2962 : /************************************************************************/
2963 : /* IReadBlock() */
2964 : /************************************************************************/
2965 :
2966 0 : CPLErr L1BNOAA15AnglesRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
2967 : int nBlockYOff, void *pData)
2968 : {
2969 0 : L1BNOAA15AnglesDataset *poGDS = (L1BNOAA15AnglesDataset *)poDS;
2970 0 : L1BDataset *poL1BDS = poGDS->poL1BDS;
2971 : int i;
2972 :
2973 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);
2974 :
2975 : /* -------------------------------------------------------------------- */
2976 : /* Seek to data. */
2977 : /* -------------------------------------------------------------------- */
2978 0 : CPL_IGNORE_RET_VAL(
2979 0 : VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));
2980 :
2981 0 : CPL_IGNORE_RET_VAL(
2982 0 : VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordSize, poL1BDS->fp));
2983 :
2984 0 : float *pafData = (float *)pData;
2985 :
2986 0 : for (i = 0; i < nRasterXSize; i++)
2987 : {
2988 : GInt16 i16 =
2989 0 : poL1BDS->GetInt16(pabyRecordHeader + 328 + 6 * i + 2 * (nBand - 1));
2990 0 : pafData[i] = i16 / 100.0f;
2991 : }
2992 :
2993 0 : if (poL1BDS->eLocationIndicator == ASCEND)
2994 : {
2995 0 : for (i = 0; i < nRasterXSize / 2; i++)
2996 : {
2997 0 : float fTmp = pafData[i];
2998 0 : pafData[i] = pafData[nRasterXSize - 1 - i];
2999 0 : pafData[nRasterXSize - 1 - i] = fTmp;
3000 : }
3001 : }
3002 :
3003 0 : CPLFree(pabyRecordHeader);
3004 :
3005 0 : return CE_None;
3006 : }
3007 :
3008 : /************************************************************************/
3009 : /* CreateAnglesDS() */
3010 : /************************************************************************/
3011 :
3012 0 : GDALDataset *L1BNOAA15AnglesDataset::CreateAnglesDS(L1BDataset *poL1BDS)
3013 : {
3014 0 : L1BNOAA15AnglesDataset *poGeolocDS = new L1BNOAA15AnglesDataset(poL1BDS);
3015 0 : for (int i = 1; i <= 3; i++)
3016 : {
3017 0 : poGeolocDS->SetBand(i, new L1BNOAA15AnglesRasterBand(poGeolocDS, i));
3018 : }
3019 0 : return poGeolocDS;
3020 : }
3021 :
3022 : /************************************************************************/
3023 : /* L1BCloudsDataset */
3024 : /************************************************************************/
3025 :
3026 : class L1BCloudsDataset final : public GDALDataset
3027 : {
3028 : friend class L1BCloudsRasterBand;
3029 :
3030 : L1BDataset *poL1BDS;
3031 :
3032 : public:
3033 : explicit L1BCloudsDataset(L1BDataset *poMainDS);
3034 : virtual ~L1BCloudsDataset();
3035 :
3036 : static GDALDataset *CreateCloudsDS(L1BDataset *poL1BDS);
3037 : };
3038 :
3039 : /************************************************************************/
3040 : /* L1BCloudsRasterBand */
3041 : /************************************************************************/
3042 :
3043 : class L1BCloudsRasterBand final : public GDALRasterBand
3044 : {
3045 : public:
3046 : L1BCloudsRasterBand(L1BCloudsDataset *poDS, int nBand);
3047 :
3048 : virtual CPLErr IReadBlock(int, int, void *) override;
3049 : };
3050 :
3051 : /************************************************************************/
3052 : /* L1BCloudsDataset() */
3053 : /************************************************************************/
3054 :
3055 0 : L1BCloudsDataset::L1BCloudsDataset(L1BDataset *poL1BDSIn) : poL1BDS(poL1BDSIn)
3056 : {
3057 0 : nRasterXSize = poL1BDSIn->nRasterXSize;
3058 0 : nRasterYSize = poL1BDSIn->nRasterYSize;
3059 0 : }
3060 :
3061 : /************************************************************************/
3062 : /* ~L1BCloudsDataset() */
3063 : /************************************************************************/
3064 :
3065 0 : L1BCloudsDataset::~L1BCloudsDataset()
3066 : {
3067 0 : delete poL1BDS;
3068 0 : }
3069 :
3070 : /************************************************************************/
3071 : /* L1BCloudsRasterBand() */
3072 : /************************************************************************/
3073 :
3074 0 : L1BCloudsRasterBand::L1BCloudsRasterBand(L1BCloudsDataset *poDSIn, int nBandIn)
3075 : {
3076 0 : poDS = poDSIn;
3077 0 : nBand = nBandIn;
3078 0 : nRasterXSize = poDSIn->nRasterXSize;
3079 0 : nRasterYSize = poDSIn->nRasterYSize;
3080 0 : eDataType = GDT_Byte;
3081 0 : nBlockXSize = nRasterXSize;
3082 0 : nBlockYSize = 1;
3083 0 : }
3084 :
3085 : /************************************************************************/
3086 : /* IReadBlock() */
3087 : /************************************************************************/
3088 :
3089 0 : CPLErr L1BCloudsRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff,
3090 : int nBlockYOff, void *pData)
3091 : {
3092 0 : L1BCloudsDataset *poGDS = (L1BCloudsDataset *)poDS;
3093 0 : L1BDataset *poL1BDS = poGDS->poL1BDS;
3094 : int i;
3095 :
3096 0 : GByte *pabyRecordHeader = (GByte *)CPLMalloc(poL1BDS->nRecordSize);
3097 :
3098 : /* -------------------------------------------------------------------- */
3099 : /* Seek to data. */
3100 : /* -------------------------------------------------------------------- */
3101 0 : CPL_IGNORE_RET_VAL(
3102 0 : VSIFSeekL(poL1BDS->fp, poL1BDS->GetLineOffset(nBlockYOff), SEEK_SET));
3103 :
3104 0 : CPL_IGNORE_RET_VAL(
3105 0 : VSIFReadL(pabyRecordHeader, 1, poL1BDS->nRecordSize, poL1BDS->fp));
3106 :
3107 0 : GByte *pabyData = (GByte *)pData;
3108 :
3109 0 : for (i = 0; i < nRasterXSize; i++)
3110 : {
3111 0 : pabyData[i] = ((pabyRecordHeader[poL1BDS->iCLAVRStart + (i / 4)] >>
3112 0 : (8 - ((i % 4) * 2 + 2))) &
3113 : 0x3);
3114 : }
3115 :
3116 0 : if (poL1BDS->eLocationIndicator == ASCEND)
3117 : {
3118 0 : for (i = 0; i < nRasterXSize / 2; i++)
3119 : {
3120 0 : GByte byTmp = pabyData[i];
3121 0 : pabyData[i] = pabyData[nRasterXSize - 1 - i];
3122 0 : pabyData[nRasterXSize - 1 - i] = byTmp;
3123 : }
3124 : }
3125 :
3126 0 : CPLFree(pabyRecordHeader);
3127 :
3128 0 : return CE_None;
3129 : }
3130 :
3131 : /************************************************************************/
3132 : /* CreateCloudsDS() */
3133 : /************************************************************************/
3134 :
3135 0 : GDALDataset *L1BCloudsDataset::CreateCloudsDS(L1BDataset *poL1BDS)
3136 : {
3137 0 : L1BCloudsDataset *poGeolocDS = new L1BCloudsDataset(poL1BDS);
3138 0 : for (int i = 1; i <= 1; i++)
3139 : {
3140 0 : poGeolocDS->SetBand(i, new L1BCloudsRasterBand(poGeolocDS, i));
3141 : }
3142 0 : return poGeolocDS;
3143 : }
3144 :
3145 : /************************************************************************/
3146 : /* DetectFormat() */
3147 : /************************************************************************/
3148 :
3149 53733 : L1BFileFormat L1BDataset::DetectFormat(const char *pszFilename,
3150 : const GByte *pabyHeader,
3151 : int nHeaderBytes)
3152 :
3153 : {
3154 53733 : if (pabyHeader == nullptr || nHeaderBytes < L1B_NOAA9_HEADER_SIZE)
3155 49185 : return L1B_NONE;
3156 :
3157 : // try NOAA-18 formats
3158 4548 : if (*(pabyHeader + 0) == '\0' && *(pabyHeader + 1) == '\0' &&
3159 759 : *(pabyHeader + 2) == '\0' && *(pabyHeader + 3) == '\0' &&
3160 140 : *(pabyHeader + 4) == '\0' && *(pabyHeader + 5) == '\0' &&
3161 130 : EQUALN((const char *)(pabyHeader + 22), "/N1BD/N18/", 10))
3162 0 : return L1B_NOAA15_NOHDR;
3163 :
3164 : // We will try the NOAA-15 and later formats first
3165 4548 : if (nHeaderBytes > L1B_NOAA15_HEADER_SIZE + 61 &&
3166 3503 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 25) == '.' &&
3167 20 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 30) == '.' &&
3168 0 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 33) == '.' &&
3169 0 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 40) == '.' &&
3170 0 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 46) == '.' &&
3171 0 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 52) == '.' &&
3172 0 : *(pabyHeader + L1B_NOAA15_HEADER_SIZE + 61) == '.')
3173 0 : return L1B_NOAA15;
3174 :
3175 : // Next try the NOAA-9/14 formats
3176 4548 : if (*(pabyHeader + 8 + 25) == '.' && *(pabyHeader + 8 + 30) == '.' &&
3177 0 : *(pabyHeader + 8 + 33) == '.' && *(pabyHeader + 8 + 40) == '.' &&
3178 0 : *(pabyHeader + 8 + 46) == '.' && *(pabyHeader + 8 + 52) == '.' &&
3179 0 : *(pabyHeader + 8 + 61) == '.')
3180 0 : return L1B_NOAA9;
3181 :
3182 : // Next try the NOAA-9/14 formats with dataset name in EBCDIC
3183 4548 : if (*(pabyHeader + 8 + 25) == 'K' && *(pabyHeader + 8 + 30) == 'K' &&
3184 0 : *(pabyHeader + 8 + 33) == 'K' && *(pabyHeader + 8 + 40) == 'K' &&
3185 0 : *(pabyHeader + 8 + 46) == 'K' && *(pabyHeader + 8 + 52) == 'K' &&
3186 0 : *(pabyHeader + 8 + 61) == 'K')
3187 0 : return L1B_NOAA9;
3188 :
3189 : // Finally try the AAPP formats
3190 4548 : if (*(pabyHeader + 25) == '.' && *(pabyHeader + 30) == '.' &&
3191 2 : *(pabyHeader + 33) == '.' && *(pabyHeader + 40) == '.' &&
3192 2 : *(pabyHeader + 46) == '.' && *(pabyHeader + 52) == '.' &&
3193 2 : *(pabyHeader + 61) == '.')
3194 2 : return L1B_NOAA15_NOHDR;
3195 :
3196 : // A few NOAA <= 9 datasets with no dataset name in TBM header
3197 4546 : if (strlen(pszFilename) == L1B_DATASET_NAME_SIZE && pszFilename[3] == '.' &&
3198 1 : pszFilename[8] == '.' && pszFilename[11] == '.' &&
3199 0 : pszFilename[18] == '.' && pszFilename[24] == '.' &&
3200 0 : pszFilename[30] == '.' && pszFilename[39] == '.' &&
3201 0 : memcmp(pabyHeader + 30,
3202 : "\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"
3203 : "\0\0\0\0\0\0\0\0\0\0\0",
3204 0 : L1B_DATASET_NAME_SIZE) == 0 &&
3205 0 : (pabyHeader[75] == '+' || pabyHeader[75] == '-') &&
3206 0 : (pabyHeader[78] == '+' || pabyHeader[78] == '-') &&
3207 0 : (pabyHeader[81] == '+' || pabyHeader[81] == '-') &&
3208 0 : (pabyHeader[85] == '+' || pabyHeader[85] == '-'))
3209 0 : return L1B_NOAA9;
3210 :
3211 4546 : return L1B_NONE;
3212 : }
3213 :
3214 : /************************************************************************/
3215 : /* Identify() */
3216 : /************************************************************************/
3217 :
3218 53732 : int L1BDataset::Identify(GDALOpenInfo *poOpenInfo)
3219 :
3220 : {
3221 53732 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS:"))
3222 0 : return TRUE;
3223 53732 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS_INTERPOL:"))
3224 0 : return TRUE;
3225 53732 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_SOLAR_ZENITH_ANGLES:"))
3226 0 : return TRUE;
3227 53732 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_ANGLES:"))
3228 0 : return TRUE;
3229 53732 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_CLOUDS:"))
3230 0 : return TRUE;
3231 :
3232 53732 : if (DetectFormat(CPLGetFilename(poOpenInfo->pszFilename),
3233 53732 : poOpenInfo->pabyHeader,
3234 53732 : poOpenInfo->nHeaderBytes) == L1B_NONE)
3235 53731 : return FALSE;
3236 :
3237 1 : return TRUE;
3238 : }
3239 :
3240 : /************************************************************************/
3241 : /* Open() */
3242 : /************************************************************************/
3243 :
3244 1 : GDALDataset *L1BDataset::Open(GDALOpenInfo *poOpenInfo)
3245 :
3246 : {
3247 1 : GDALDataset *poOutDS = nullptr;
3248 1 : VSILFILE *fp = nullptr;
3249 2 : CPLString osFilename = poOpenInfo->pszFilename;
3250 1 : int bAskGeolocationDS = FALSE;
3251 1 : int bInterpolGeolocationDS = FALSE;
3252 1 : int bAskSolarZenithAnglesDS = FALSE;
3253 1 : int bAskAnglesDS = FALSE;
3254 1 : int bAskCloudsDS = FALSE;
3255 : L1BFileFormat eL1BFormat;
3256 :
3257 1 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS:") ||
3258 1 : STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS_INTERPOL:") ||
3259 1 : STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_SOLAR_ZENITH_ANGLES:") ||
3260 1 : STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_ANGLES:") ||
3261 1 : STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_CLOUDS:"))
3262 : {
3263 : GByte abyHeader[1024];
3264 0 : const char *pszFilename = nullptr;
3265 0 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS_INTERPOL:"))
3266 : {
3267 0 : bAskGeolocationDS = TRUE;
3268 0 : bInterpolGeolocationDS = TRUE;
3269 0 : pszFilename = poOpenInfo->pszFilename + strlen("L1BGCPS_INTERPOL:");
3270 : }
3271 0 : else if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1BGCPS:"))
3272 : {
3273 0 : bAskGeolocationDS = TRUE;
3274 0 : pszFilename = poOpenInfo->pszFilename + strlen("L1BGCPS:");
3275 : }
3276 0 : else if (STARTS_WITH_CI(poOpenInfo->pszFilename,
3277 : "L1B_SOLAR_ZENITH_ANGLES:"))
3278 : {
3279 0 : bAskSolarZenithAnglesDS = TRUE;
3280 0 : pszFilename =
3281 0 : poOpenInfo->pszFilename + strlen("L1B_SOLAR_ZENITH_ANGLES:");
3282 : }
3283 0 : else if (STARTS_WITH_CI(poOpenInfo->pszFilename, "L1B_ANGLES:"))
3284 : {
3285 0 : bAskAnglesDS = TRUE;
3286 0 : pszFilename = poOpenInfo->pszFilename + strlen("L1B_ANGLES:");
3287 : }
3288 : else
3289 : {
3290 0 : bAskCloudsDS = TRUE;
3291 0 : pszFilename = poOpenInfo->pszFilename + strlen("L1B_CLOUDS:");
3292 : }
3293 0 : if (pszFilename[0] == '"')
3294 0 : pszFilename++;
3295 0 : osFilename = pszFilename;
3296 0 : if (!osFilename.empty() && osFilename.back() == '"')
3297 0 : osFilename.pop_back();
3298 0 : fp = VSIFOpenL(osFilename, "rb");
3299 0 : if (!fp)
3300 : {
3301 0 : CPLError(CE_Failure, CPLE_AppDefined, "Can't open file \"%s\".",
3302 : osFilename.c_str());
3303 0 : return nullptr;
3304 : }
3305 0 : CPL_IGNORE_RET_VAL(VSIFReadL(abyHeader, 1, sizeof(abyHeader) - 1, fp));
3306 0 : abyHeader[sizeof(abyHeader) - 1] = '\0';
3307 0 : eL1BFormat = DetectFormat(CPLGetFilename(osFilename), abyHeader,
3308 0 : sizeof(abyHeader));
3309 : }
3310 : else
3311 : eL1BFormat =
3312 1 : DetectFormat(CPLGetFilename(osFilename), poOpenInfo->pabyHeader,
3313 : poOpenInfo->nHeaderBytes);
3314 :
3315 1 : if (eL1BFormat == L1B_NONE)
3316 : {
3317 0 : if (fp != nullptr)
3318 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
3319 0 : return nullptr;
3320 : }
3321 :
3322 : /* -------------------------------------------------------------------- */
3323 : /* Confirm the requested access is supported. */
3324 : /* -------------------------------------------------------------------- */
3325 1 : if (poOpenInfo->eAccess == GA_Update)
3326 : {
3327 0 : CPLError(CE_Failure, CPLE_NotSupported,
3328 : "The L1B driver does not support update access to existing"
3329 : " datasets.\n");
3330 0 : if (fp != nullptr)
3331 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
3332 0 : return nullptr;
3333 : }
3334 :
3335 : /* -------------------------------------------------------------------- */
3336 : /* Create a corresponding GDALDataset. */
3337 : /* -------------------------------------------------------------------- */
3338 : L1BDataset *poDS;
3339 : VSIStatBufL sStat;
3340 :
3341 1 : poDS = new L1BDataset(eL1BFormat);
3342 :
3343 1 : if (fp == nullptr)
3344 1 : fp = VSIFOpenL(osFilename, "rb");
3345 1 : poDS->fp = fp;
3346 1 : if (!poDS->fp || VSIStatL(osFilename, &sStat) != 0)
3347 : {
3348 0 : CPLDebug("L1B", "Can't open file \"%s\".", osFilename.c_str());
3349 0 : goto bad;
3350 : }
3351 :
3352 : /* -------------------------------------------------------------------- */
3353 : /* Read the header. */
3354 : /* -------------------------------------------------------------------- */
3355 1 : if (poDS->ProcessDatasetHeader(CPLGetFilename(osFilename)) != CE_None)
3356 : {
3357 0 : CPLDebug("L1B", "Error reading L1B record header.");
3358 0 : goto bad;
3359 : }
3360 :
3361 1 : if (poDS->eL1BFormat == L1B_NOAA15_NOHDR &&
3362 1 : poDS->nRecordSizeFromHeader == 22016 &&
3363 1 : (sStat.st_size % 22016 /* poDS->nRecordSizeFromHeader*/) == 0)
3364 : {
3365 1 : poDS->iDataFormat = UNPACKED16BIT;
3366 1 : poDS->ComputeFileOffsets();
3367 1 : poDS->nDataStartOffset = poDS->nRecordSizeFromHeader;
3368 1 : poDS->nRecordSize = poDS->nRecordSizeFromHeader;
3369 1 : poDS->iCLAVRStart = 0;
3370 : }
3371 0 : else if (poDS->bGuessDataFormat)
3372 : {
3373 : int nTempYSize;
3374 : GUInt16 nScanlineNumber;
3375 : int j;
3376 :
3377 : /* If the data format is unspecified, try each one of the 3 known data
3378 : * formats */
3379 : /* It is considered valid when the spacing between the first 5 scanline
3380 : * numbers */
3381 : /* is a constant */
3382 :
3383 0 : for (j = 0; j < 3; j++)
3384 : {
3385 0 : poDS->iDataFormat = (L1BDataFormat)(PACKED10BIT + j);
3386 0 : if (!poDS->ComputeFileOffsets())
3387 0 : goto bad;
3388 :
3389 0 : nTempYSize = static_cast<int>(
3390 0 : (sStat.st_size - poDS->nDataStartOffset) / poDS->nRecordSize);
3391 0 : if (nTempYSize < 5)
3392 0 : continue;
3393 :
3394 0 : int nLastScanlineNumber = 0;
3395 0 : int nDiffLine = 0;
3396 : int i;
3397 0 : for (i = 0; i < 5; i++)
3398 : {
3399 0 : nScanlineNumber = 0;
3400 :
3401 0 : CPL_IGNORE_RET_VAL(VSIFSeekL(
3402 0 : poDS->fp, poDS->nDataStartOffset + i * poDS->nRecordSize,
3403 : SEEK_SET));
3404 0 : CPL_IGNORE_RET_VAL(VSIFReadL(&nScanlineNumber, 1, 2, poDS->fp));
3405 0 : nScanlineNumber = poDS->GetUInt16(&nScanlineNumber);
3406 :
3407 0 : if (i == 1)
3408 : {
3409 0 : nDiffLine = nScanlineNumber - nLastScanlineNumber;
3410 0 : if (nDiffLine == 0)
3411 0 : break;
3412 : }
3413 0 : else if (i > 1)
3414 : {
3415 0 : if (nDiffLine != nScanlineNumber - nLastScanlineNumber)
3416 0 : break;
3417 : }
3418 :
3419 0 : nLastScanlineNumber = nScanlineNumber;
3420 : }
3421 :
3422 0 : if (i == 5)
3423 : {
3424 0 : CPLDebug("L1B", "Guessed data format : %s",
3425 0 : (poDS->iDataFormat == PACKED10BIT) ? "10"
3426 0 : : (poDS->iDataFormat == UNPACKED8BIT) ? "08"
3427 : : "16");
3428 0 : break;
3429 : }
3430 : }
3431 :
3432 0 : if (j == 3)
3433 : {
3434 0 : CPLError(CE_Failure, CPLE_AppDefined,
3435 : "Could not guess data format of L1B product");
3436 0 : goto bad;
3437 : }
3438 : }
3439 : else
3440 : {
3441 0 : if (!poDS->ComputeFileOffsets())
3442 0 : goto bad;
3443 : }
3444 :
3445 1 : CPLDebug("L1B", "nRecordDataStart = %d", poDS->nRecordDataStart);
3446 1 : CPLDebug("L1B", "nRecordDataEnd = %d", poDS->nRecordDataEnd);
3447 1 : CPLDebug("L1B", "nDataStartOffset = %d", poDS->nDataStartOffset);
3448 1 : CPLDebug("L1B", "iCLAVRStart = %d", poDS->iCLAVRStart);
3449 1 : CPLDebug("L1B", "nRecordSize = %d", poDS->nRecordSize);
3450 :
3451 : // Compute number of lines dynamically, so we can read partially
3452 : // downloaded files.
3453 1 : if (poDS->nDataStartOffset > sStat.st_size)
3454 0 : goto bad;
3455 1 : poDS->nRasterYSize = static_cast<int>(
3456 1 : (sStat.st_size - poDS->nDataStartOffset) / poDS->nRecordSize);
3457 :
3458 : /* -------------------------------------------------------------------- */
3459 : /* Deal with GCPs */
3460 : /* -------------------------------------------------------------------- */
3461 1 : poDS->ProcessRecordHeaders();
3462 :
3463 1 : if (bAskGeolocationDS)
3464 : {
3465 0 : return L1BGeolocDataset::CreateGeolocationDS(poDS,
3466 0 : bInterpolGeolocationDS);
3467 : }
3468 1 : else if (bAskSolarZenithAnglesDS)
3469 : {
3470 0 : if (eL1BFormat == L1B_NOAA9)
3471 0 : return L1BSolarZenithAnglesDataset::CreateSolarZenithAnglesDS(poDS);
3472 : else
3473 : {
3474 0 : delete poDS;
3475 0 : return nullptr;
3476 : }
3477 : }
3478 1 : else if (bAskAnglesDS)
3479 : {
3480 0 : if (eL1BFormat != L1B_NOAA9)
3481 0 : return L1BNOAA15AnglesDataset::CreateAnglesDS(poDS);
3482 : else
3483 : {
3484 0 : delete poDS;
3485 0 : return nullptr;
3486 : }
3487 : }
3488 1 : else if (bAskCloudsDS)
3489 : {
3490 0 : if (poDS->iCLAVRStart > 0)
3491 0 : poOutDS = L1BCloudsDataset::CreateCloudsDS(poDS);
3492 : else
3493 : {
3494 0 : delete poDS;
3495 0 : return nullptr;
3496 : }
3497 : }
3498 : else
3499 : {
3500 1 : poOutDS = poDS;
3501 : }
3502 :
3503 : {
3504 2 : CPLString osTMP;
3505 : int bInterpol =
3506 1 : CPLTestBool(CPLGetConfigOption("L1B_INTERPOL_GCPS", "TRUE"));
3507 :
3508 1 : char *pszWKT = nullptr;
3509 1 : poDS->m_oGCPSRS.exportToWkt(&pszWKT);
3510 1 : poOutDS->SetMetadataItem("SRS", pszWKT,
3511 1 : "GEOLOCATION"); /* unused by gdalgeoloc.cpp */
3512 1 : CPLFree(pszWKT);
3513 :
3514 1 : if (bInterpol)
3515 1 : osTMP.Printf("L1BGCPS_INTERPOL:\"%s\"", osFilename.c_str());
3516 : else
3517 0 : osTMP.Printf("L1BGCPS:\"%s\"", osFilename.c_str());
3518 1 : poOutDS->SetMetadataItem("X_DATASET", osTMP, "GEOLOCATION");
3519 1 : poOutDS->SetMetadataItem("X_BAND", "1", "GEOLOCATION");
3520 1 : poOutDS->SetMetadataItem("Y_DATASET", osTMP, "GEOLOCATION");
3521 1 : poOutDS->SetMetadataItem("Y_BAND", "2", "GEOLOCATION");
3522 :
3523 1 : if (bInterpol)
3524 : {
3525 1 : poOutDS->SetMetadataItem("PIXEL_OFFSET", "0", "GEOLOCATION");
3526 1 : poOutDS->SetMetadataItem("PIXEL_STEP", "1", "GEOLOCATION");
3527 : }
3528 : else
3529 : {
3530 0 : osTMP.Printf("%d", poDS->iGCPStart);
3531 0 : poOutDS->SetMetadataItem("PIXEL_OFFSET", osTMP, "GEOLOCATION");
3532 0 : osTMP.Printf("%d", poDS->iGCPStep);
3533 0 : poOutDS->SetMetadataItem("PIXEL_STEP", osTMP, "GEOLOCATION");
3534 : }
3535 :
3536 1 : poOutDS->SetMetadataItem("LINE_OFFSET", "0", "GEOLOCATION");
3537 1 : poOutDS->SetMetadataItem("LINE_STEP", "1", "GEOLOCATION");
3538 : }
3539 :
3540 1 : if (poOutDS != poDS)
3541 0 : return poOutDS;
3542 :
3543 1 : if (eL1BFormat == L1B_NOAA9)
3544 : {
3545 0 : char **papszSubdatasets = nullptr;
3546 0 : papszSubdatasets = CSLSetNameValue(
3547 : papszSubdatasets, "SUBDATASET_1_NAME",
3548 : CPLSPrintf("L1B_SOLAR_ZENITH_ANGLES:\"%s\"", osFilename.c_str()));
3549 0 : papszSubdatasets = CSLSetNameValue(
3550 : papszSubdatasets, "SUBDATASET_1_DESC", "Solar zenith angles");
3551 0 : poDS->SetMetadata(papszSubdatasets, "SUBDATASETS");
3552 0 : CSLDestroy(papszSubdatasets);
3553 : }
3554 : else
3555 : {
3556 1 : char **papszSubdatasets = nullptr;
3557 1 : papszSubdatasets = CSLSetNameValue(
3558 : papszSubdatasets, "SUBDATASET_1_NAME",
3559 : CPLSPrintf("L1B_ANGLES:\"%s\"", osFilename.c_str()));
3560 : papszSubdatasets =
3561 1 : CSLSetNameValue(papszSubdatasets, "SUBDATASET_1_DESC",
3562 : "Solar zenith angles, satellite zenith angles and "
3563 : "relative azimuth angles");
3564 :
3565 1 : if (poDS->iCLAVRStart > 0)
3566 : {
3567 0 : papszSubdatasets = CSLSetNameValue(
3568 : papszSubdatasets, "SUBDATASET_2_NAME",
3569 : CPLSPrintf("L1B_CLOUDS:\"%s\"", osFilename.c_str()));
3570 : papszSubdatasets =
3571 0 : CSLSetNameValue(papszSubdatasets, "SUBDATASET_2_DESC",
3572 : "Clouds from AVHRR (CLAVR)");
3573 : }
3574 :
3575 1 : poDS->SetMetadata(papszSubdatasets, "SUBDATASETS");
3576 1 : CSLDestroy(papszSubdatasets);
3577 : }
3578 :
3579 : /* -------------------------------------------------------------------- */
3580 : /* Create band information objects. */
3581 : /* -------------------------------------------------------------------- */
3582 : int iBand, i;
3583 :
3584 6 : for (iBand = 1, i = 0; iBand <= poDS->nBands; iBand++)
3585 : {
3586 5 : poDS->SetBand(iBand, new L1BRasterBand(poDS, iBand));
3587 :
3588 : // Channels descriptions
3589 5 : if (poDS->eSpacecraftID >= NOAA6 && poDS->eSpacecraftID <= METOP3)
3590 : {
3591 5 : if (!(i & 0x01) && poDS->iChannelsMask & 0x01)
3592 : {
3593 1 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[0]);
3594 1 : i |= 0x01;
3595 1 : continue;
3596 : }
3597 4 : if (!(i & 0x02) && poDS->iChannelsMask & 0x02)
3598 : {
3599 1 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[1]);
3600 1 : i |= 0x02;
3601 1 : continue;
3602 : }
3603 3 : if (!(i & 0x04) && poDS->iChannelsMask & 0x04)
3604 : {
3605 1 : if (poDS->eSpacecraftID >= NOAA15 &&
3606 1 : poDS->eSpacecraftID <= METOP3)
3607 1 : if (poDS->iInstrumentStatus & 0x0400)
3608 1 : poDS->GetRasterBand(iBand)->SetDescription(
3609 1 : apszBandDesc[7]);
3610 : else
3611 0 : poDS->GetRasterBand(iBand)->SetDescription(
3612 0 : apszBandDesc[6]);
3613 : else
3614 0 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[2]);
3615 1 : i |= 0x04;
3616 1 : continue;
3617 : }
3618 2 : if (!(i & 0x08) && poDS->iChannelsMask & 0x08)
3619 : {
3620 1 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[3]);
3621 1 : i |= 0x08;
3622 1 : continue;
3623 : }
3624 1 : if (!(i & 0x10) && poDS->iChannelsMask & 0x10)
3625 : {
3626 1 : if (poDS->eSpacecraftID == NOAA13) // 5 NOAA-13
3627 0 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[5]);
3628 1 : else if (poDS->eSpacecraftID == NOAA6 ||
3629 1 : poDS->eSpacecraftID == NOAA8 ||
3630 1 : poDS->eSpacecraftID == NOAA10) // 4 NOAA-6,-8,-10
3631 0 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[3]);
3632 : else
3633 1 : poDS->GetRasterBand(iBand)->SetDescription(apszBandDesc[4]);
3634 1 : i |= 0x10;
3635 1 : continue;
3636 : }
3637 : }
3638 : }
3639 :
3640 1 : if (poDS->bExposeMaskBand)
3641 1 : poDS->poMaskBand = new L1BMaskBand(poDS);
3642 :
3643 : /* -------------------------------------------------------------------- */
3644 : /* Initialize any PAM information. */
3645 : /* -------------------------------------------------------------------- */
3646 1 : poDS->SetDescription(poOpenInfo->pszFilename);
3647 1 : poDS->TryLoadXML();
3648 :
3649 : /* -------------------------------------------------------------------- */
3650 : /* Check for external overviews. */
3651 : /* -------------------------------------------------------------------- */
3652 2 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
3653 1 : poOpenInfo->GetSiblingFiles());
3654 :
3655 : /* -------------------------------------------------------------------- */
3656 : /* Fetch metadata in CSV file */
3657 : /* -------------------------------------------------------------------- */
3658 1 : if (CPLTestBool(CPLGetConfigOption("L1B_FETCH_METADATA", "NO")))
3659 : {
3660 0 : poDS->FetchMetadata();
3661 : }
3662 :
3663 1 : return poDS;
3664 :
3665 0 : bad:
3666 0 : delete poDS;
3667 0 : return nullptr;
3668 : }
3669 :
3670 : /************************************************************************/
3671 : /* GDALRegister_L1B() */
3672 : /************************************************************************/
3673 :
3674 1595 : void GDALRegister_L1B()
3675 :
3676 : {
3677 1595 : if (GDALGetDriverByName("L1B") != nullptr)
3678 302 : return;
3679 :
3680 1293 : GDALDriver *poDriver = new GDALDriver();
3681 :
3682 1293 : poDriver->SetDescription("L1B");
3683 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
3684 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
3685 1293 : "NOAA Polar Orbiter Level 1b Data Set");
3686 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/l1b.html");
3687 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
3688 1293 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
3689 :
3690 1293 : poDriver->pfnOpen = L1BDataset::Open;
3691 1293 : poDriver->pfnIdentify = L1BDataset::Identify;
3692 :
3693 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
3694 : }
|
