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