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