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