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