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