1 /******************************************************************************
2 *
3 * Project: APP ENVISAT Support
4 * Purpose: Reader for ENVISAT format image data.
5 * Author: Frank Warmerdam, warmerdam@pobox.com
6 *
7 ******************************************************************************
8 * Copyright (c) 2001, Atlantis Scientific, Inc.
9 * Copyright (c) 2009-2013, Even Rouault <even dot rouault at spatialys.com>
10 *
11 * Permission is hereby granted, free of charge, to any person obtaining a
12 * copy of this software and associated documentation files (the "Software"),
13 * to deal in the Software without restriction, including without limitation
14 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 * and/or sell copies of the Software, and to permit persons to whom the
16 * Software is furnished to do so, subject to the following conditions:
17 *
18 * The above copyright notice and this permission notice shall be included
19 * in all copies or substantial portions of the Software.
20 *
21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 * DEALINGS IN THE SOFTWARE.
28 ****************************************************************************/
29
30 #include "adsrange.hpp"
31 #include "rawdataset.h"
32 #include "cpl_string.h"
33 #include "gdal_frmts.h"
34 #include "ogr_srs_api.h"
35 #include "timedelta.hpp"
36
37 CPL_CVSID("$Id: envisatdataset.cpp fa752ad6eabafaf630a704e1892a9d837d683cb3 2021-03-06 17:04:38 +0100 Even Rouault $")
38
39 CPL_C_START
40 #include "EnvisatFile.h"
41 #include "records.h"
42 CPL_C_END
43
44 /************************************************************************/
45 /* ==================================================================== */
46 /* MerisL2FlagBand */
47 /* ==================================================================== */
48 /************************************************************************/
49 class MerisL2FlagBand final: public GDALPamRasterBand
50 {
51 public:
52 MerisL2FlagBand( GDALDataset *, int, VSILFILE*, vsi_l_offset, int );
53 ~MerisL2FlagBand() override;
54 virtual CPLErr IReadBlock( int, int, void * ) override;
55
56 private:
57 vsi_l_offset nImgOffset;
58 int nPrefixBytes;
59 size_t nBytePerPixel;
60 size_t nRecordSize;
61 size_t nDataSize;
62 GByte *pReadBuf;
63 VSILFILE *fpImage;
64 };
65
66 /************************************************************************/
67 /* MerisL2FlagBand() */
68 /************************************************************************/
MerisL2FlagBand(GDALDataset * poDSIn,int nBandIn,VSILFILE * fpImageIn,vsi_l_offset nImgOffsetIn,int nPrefixBytesIn)69 MerisL2FlagBand::MerisL2FlagBand( GDALDataset *poDSIn, int nBandIn,
70 VSILFILE* fpImageIn,
71 vsi_l_offset nImgOffsetIn,
72 int nPrefixBytesIn ) :
73 nImgOffset(nImgOffsetIn),
74 nPrefixBytes(nPrefixBytesIn),
75 nBytePerPixel(3),
76 nRecordSize(0),
77 nDataSize(0),
78 pReadBuf(nullptr)
79 {
80 poDS = poDSIn;
81 nBand = nBandIn;
82
83 fpImage = fpImageIn;
84
85 eDataType = GDT_UInt32;
86
87 nBlockXSize = poDS->GetRasterXSize();
88 nBlockYSize = 1;
89 nRecordSize = nPrefixBytesIn + nBlockXSize * nBytePerPixel;
90 nDataSize = nBlockXSize * nBytePerPixel;
91 pReadBuf = static_cast<GByte *>(CPLMalloc(nRecordSize));
92 }
93
94 /************************************************************************/
95 /* ~MerisL2FlagBand() */
96 /************************************************************************/
~MerisL2FlagBand()97 MerisL2FlagBand::~MerisL2FlagBand()
98 {
99 CPLFree( pReadBuf );
100 }
101
102 /************************************************************************/
103 /* IReadBlock() */
104 /************************************************************************/
IReadBlock(CPL_UNUSED int nBlockXOff,int nBlockYOff,void * pImage)105 CPLErr MerisL2FlagBand::IReadBlock( CPL_UNUSED int nBlockXOff,
106 int nBlockYOff,
107 void * pImage )
108 {
109 CPLAssert( nBlockXOff == 0 );
110 CPLAssert( pReadBuf != nullptr );
111
112 vsi_l_offset nOffset = nImgOffset + nPrefixBytes +
113 nBlockYOff * nBlockYSize * nRecordSize;
114
115 if ( VSIFSeekL( fpImage, nOffset, SEEK_SET ) != 0 )
116 {
117 CPLError( CE_Failure, CPLE_FileIO,
118 "Seek to %d for scanline %d failed.\n",
119 (int)nOffset, nBlockYOff );
120 return CE_Failure;
121 }
122
123 if ( VSIFReadL( pReadBuf, 1, nDataSize, fpImage ) != nDataSize )
124 {
125 CPLError( CE_Failure, CPLE_FileIO,
126 "Read of %d bytes for scanline %d failed.\n",
127 (int)nDataSize, nBlockYOff );
128 return CE_Failure;
129 }
130
131 const unsigned int nUInt32Size = 4;
132 for( unsigned iImg = 0, iBuf = 0;
133 iImg < nBlockXSize * nUInt32Size;
134 iImg += nUInt32Size, iBuf += (unsigned)nBytePerPixel )
135 {
136 #ifdef CPL_LSB
137 ((GByte*) pImage)[iImg] = pReadBuf[iBuf + 2];
138 ((GByte*) pImage)[iImg + 1] = pReadBuf[iBuf + 1];
139 ((GByte*) pImage)[iImg + 2] = pReadBuf[iBuf];
140 ((GByte*) pImage)[iImg + 3] = 0;
141 #else
142 ((GByte*) pImage)[iImg] = 0;
143 ((GByte*) pImage)[iImg + 1] = pReadBuf[iBuf];
144 ((GByte*) pImage)[iImg + 2] = pReadBuf[iBuf + 1];
145 ((GByte*) pImage)[iImg + 3] = pReadBuf[iBuf + 2];
146 #endif
147 }
148
149 return CE_None;
150 }
151
152 /************************************************************************/
153 /* ==================================================================== */
154 /* EnvisatDataset */
155 /* ==================================================================== */
156 /************************************************************************/
157
158 class EnvisatDataset final: public RawDataset
159 {
160 EnvisatFile *hEnvisatFile;
161 VSILFILE *fpImage;
162
163 int nGCPCount;
164 GDAL_GCP *pasGCPList;
165
166 char **papszTempMD;
167
168 void ScanForGCPs_ASAR();
169 void ScanForGCPs_MERIS();
170
171 void UnwrapGCPs();
172
173 void CollectMetadata( EnvisatFile_HeaderFlag );
174 void CollectDSDMetadata();
175 void CollectADSMetadata();
176
177 public:
178 EnvisatDataset();
179 virtual ~EnvisatDataset();
180
181 virtual int GetGCPCount() override;
182 virtual const char *_GetGCPProjection() override;
GetGCPSpatialRef() const183 const OGRSpatialReference* GetGCPSpatialRef() const override {
184 return GetGCPSpatialRefFromOldGetGCPProjection();
185 }
186 virtual const GDAL_GCP *GetGCPs() override;
187 virtual char **GetMetadataDomainList() override;
188 virtual char **GetMetadata( const char * pszDomain ) override;
189
190 static GDALDataset *Open( GDALOpenInfo * );
191 };
192
193 /************************************************************************/
194 /* ==================================================================== */
195 /* EnvisatDataset */
196 /* ==================================================================== */
197 /************************************************************************/
198
199 /************************************************************************/
200 /* EnvisatDataset() */
201 /************************************************************************/
202
EnvisatDataset()203 EnvisatDataset::EnvisatDataset() :
204 hEnvisatFile(nullptr),
205 fpImage(nullptr),
206 nGCPCount(0),
207 pasGCPList(nullptr),
208 papszTempMD(nullptr)
209 {}
210
211 /************************************************************************/
212 /* ~EnvisatDataset() */
213 /************************************************************************/
214
~EnvisatDataset()215 EnvisatDataset::~EnvisatDataset()
216
217 {
218 FlushCache();
219
220 if( hEnvisatFile != nullptr )
221 EnvisatFile_Close( hEnvisatFile );
222
223 if( fpImage != nullptr )
224 CPL_IGNORE_RET_VAL(VSIFCloseL( fpImage ));
225
226 if( nGCPCount > 0 )
227 {
228 GDALDeinitGCPs( nGCPCount, pasGCPList );
229 CPLFree( pasGCPList );
230 }
231
232 CSLDestroy( papszTempMD );
233 }
234
235 /************************************************************************/
236 /* GetGCPCount() */
237 /************************************************************************/
238
GetGCPCount()239 int EnvisatDataset::GetGCPCount()
240
241 {
242 return nGCPCount;
243 }
244
245 /************************************************************************/
246 /* GetGCPProjection() */
247 /************************************************************************/
248
_GetGCPProjection()249 const char *EnvisatDataset::_GetGCPProjection()
250
251 {
252 if( nGCPCount > 0 )
253 return SRS_WKT_WGS84_LAT_LONG;
254
255 return "";
256 }
257
258 /************************************************************************/
259 /* GetGCP() */
260 /************************************************************************/
261
GetGCPs()262 const GDAL_GCP *EnvisatDataset::GetGCPs()
263
264 {
265 return pasGCPList;
266 }
267
268 /************************************************************************/
269 /* UnwrapGCPs() */
270 /************************************************************************/
271
272 /* external C++ implementation of the in-place unwrapper */
273 void EnvisatUnwrapGCPs( int nGCPCount, GDAL_GCP *pasGCPList ) ;
274
UnwrapGCPs()275 void EnvisatDataset::UnwrapGCPs()
276 {
277 EnvisatUnwrapGCPs( nGCPCount, pasGCPList ) ;
278 }
279
280 /************************************************************************/
281 /* ScanForGCPs_ASAR() */
282 /************************************************************************/
283
ScanForGCPs_ASAR()284 void EnvisatDataset::ScanForGCPs_ASAR()
285
286 {
287 /* -------------------------------------------------------------------- */
288 /* Do we have a meaningful geolocation grid? */
289 /* -------------------------------------------------------------------- */
290 int nDatasetIndex = EnvisatFile_GetDatasetIndex( hEnvisatFile,
291 "GEOLOCATION GRID ADS" );
292 if( nDatasetIndex == -1 )
293 return;
294
295 int nNumDSR, nDSRSize;
296 if( EnvisatFile_GetDatasetInfo( hEnvisatFile, nDatasetIndex,
297 nullptr, nullptr, nullptr, nullptr, nullptr,
298 &nNumDSR, &nDSRSize ) != SUCCESS )
299 return;
300
301 if( nNumDSR == 0 || nDSRSize != 521 )
302 return;
303
304 /* -------------------------------------------------------------------- */
305 /* Collect the first GCP set from each record. */
306 /* -------------------------------------------------------------------- */
307 GByte abyRecord[521];
308 int nRange = 0;
309 int nRangeOffset = 0;
310 GUInt32 unValue;
311
312 nGCPCount = 0;
313 pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),(nNumDSR+1) * 11);
314
315 for( int iRecord = 0; iRecord < nNumDSR; iRecord++ )
316 {
317 if( EnvisatFile_ReadDatasetRecord( hEnvisatFile, nDatasetIndex,
318 iRecord, abyRecord ) != SUCCESS )
319 continue;
320
321 memcpy( &unValue, abyRecord + 13, 4 );
322 nRange = CPL_MSBWORD32( unValue ) + nRangeOffset;
323
324 if((iRecord>1) && (int(pasGCPList[nGCPCount-1].dfGCPLine+0.5) > nRange))
325 {
326 int delta = (int) (pasGCPList[nGCPCount-1].dfGCPLine -
327 pasGCPList[nGCPCount-12].dfGCPLine);
328 nRange = int(pasGCPList[nGCPCount-1].dfGCPLine+0.5) + delta;
329 nRangeOffset = nRange-1;
330 }
331
332 for( int iGCP = 0; iGCP < 11; iGCP++ )
333 {
334 GDALInitGCPs( 1, pasGCPList + nGCPCount );
335
336 CPLFree( pasGCPList[nGCPCount].pszId );
337
338 char szId[128];
339 snprintf( szId, sizeof(szId), "%d", nGCPCount+1 );
340 pasGCPList[nGCPCount].pszId = CPLStrdup( szId );
341
342 memcpy( &unValue, abyRecord + 25 + iGCP*4, 4 );
343 int nSample = CPL_MSBWORD32(unValue);
344
345 memcpy( &unValue, abyRecord + 25 + 176 + iGCP*4, 4 );
346 pasGCPList[nGCPCount].dfGCPX = ((int)CPL_MSBWORD32(unValue))*0.000001;
347
348 memcpy( &unValue, abyRecord + 25 + 132 + iGCP*4, 4 );
349 pasGCPList[nGCPCount].dfGCPY = ((int)CPL_MSBWORD32(unValue))*0.000001;
350
351 pasGCPList[nGCPCount].dfGCPZ = 0.0;
352
353 pasGCPList[nGCPCount].dfGCPLine = nRange - 0.5;
354 pasGCPList[nGCPCount].dfGCPPixel = nSample - 0.5;
355
356 nGCPCount++;
357 }
358 }
359
360 /* -------------------------------------------------------------------- */
361 /* We also collect the bottom GCPs from the last granule. */
362 /* -------------------------------------------------------------------- */
363 memcpy( &unValue, abyRecord + 17, 4 );
364 nRange = nRange + CPL_MSBWORD32( unValue ) - 1;
365
366 for( int iGCP = 0; iGCP < 11; iGCP++ )
367 {
368 GDALInitGCPs( 1, pasGCPList + nGCPCount );
369
370 CPLFree( pasGCPList[nGCPCount].pszId );
371
372 char szId[128];
373 snprintf( szId, sizeof(szId), "%d", nGCPCount+1 );
374 pasGCPList[nGCPCount].pszId = CPLStrdup( szId );
375
376 memcpy( &unValue, abyRecord + 279 + iGCP*4, 4 );
377 int nSample = CPL_MSBWORD32(unValue);
378
379 memcpy( &unValue, abyRecord + 279 + 176 + iGCP*4, 4 );
380 pasGCPList[nGCPCount].dfGCPX = ((int)CPL_MSBWORD32(unValue))*0.000001;
381
382 memcpy( &unValue, abyRecord + 279 + 132 + iGCP*4, 4 );
383 pasGCPList[nGCPCount].dfGCPY = ((int)CPL_MSBWORD32(unValue))*0.000001;
384
385 pasGCPList[nGCPCount].dfGCPZ = 0.0;
386
387 pasGCPList[nGCPCount].dfGCPLine = nRange - 0.5;
388 pasGCPList[nGCPCount].dfGCPPixel = nSample - 0.5;
389
390 nGCPCount++;
391 }
392 }
393
394 /************************************************************************/
395 /* ScanForGCPs_MERIS() */
396 /************************************************************************/
397
ScanForGCPs_MERIS()398 void EnvisatDataset::ScanForGCPs_MERIS()
399
400 {
401 /* -------------------------------------------------------------------- */
402 /* Do we have a meaningful geolocation grid? Search for a */
403 /* DS_TYPE=A and a name containing "geolocation" or "tie */
404 /* points". */
405 /* -------------------------------------------------------------------- */
406 int nDatasetIndex = EnvisatFile_GetDatasetIndex( hEnvisatFile,
407 "Tie points ADS" );
408 if( nDatasetIndex == -1 )
409 return;
410
411 int nNumDSR, nDSRSize;
412 if( EnvisatFile_GetDatasetInfo( hEnvisatFile, nDatasetIndex,
413 nullptr, nullptr, nullptr, nullptr, nullptr,
414 &nNumDSR, &nDSRSize ) != SUCCESS )
415 return;
416
417 if( nNumDSR == 0 )
418 return;
419
420 /* -------------------------------------------------------------------- */
421 /* Figure out the tiepoint space, and how many we have. */
422 /* -------------------------------------------------------------------- */
423 int nLinesPerTiePoint =
424 EnvisatFile_GetKeyValueAsInt( hEnvisatFile, SPH,
425 "LINES_PER_TIE_PT", 0 );
426 int nSamplesPerTiePoint =
427 EnvisatFile_GetKeyValueAsInt( hEnvisatFile, SPH,
428 "SAMPLES_PER_TIE_PT", 0 );
429
430 if( nLinesPerTiePoint == 0 || nSamplesPerTiePoint == 0 )
431 return;
432
433 int nTPPerColumn = nNumDSR;
434 int nTPPerLine = (GetRasterXSize() + nSamplesPerTiePoint - 1)
435 / nSamplesPerTiePoint;
436
437 /* -------------------------------------------------------------------- */
438 /* Find a measurement type dataset to use as a reference raster */
439 /* band. */
440 /* -------------------------------------------------------------------- */
441
442 int nMDSIndex = 0;
443
444 for( ; true; nMDSIndex++ )
445 {
446 char *pszDSType = nullptr;
447 if( EnvisatFile_GetDatasetInfo( hEnvisatFile, nMDSIndex,
448 nullptr, &pszDSType, nullptr, nullptr, nullptr, nullptr, nullptr ) == FAILURE )
449 {
450 CPLDebug("EnvisatDataset",
451 "Unable to find MDS in Envisat file.") ;
452 return ;
453 }
454 if( EQUAL(pszDSType,"M") ) break;
455 }
456
457 /* -------------------------------------------------------------------- */
458 /* Get subset of TP ADS records matching the MDS records */
459 /* -------------------------------------------------------------------- */
460
461 /* get the MDS line sampling time interval */
462 TimeDelta tdMDSSamplingInterval( 0 , 0 ,
463 EnvisatFile_GetKeyValueAsInt( hEnvisatFile, SPH,
464 "LINE_TIME_INTERVAL", 0 ) );
465
466 /* get range of TiePoint ADS records matching the measurements */
467 ADSRangeLastAfter arTP( *hEnvisatFile , nDatasetIndex,
468 nMDSIndex , tdMDSSamplingInterval ) ;
469
470 /* check if there are any TPs to be used */
471 if ( arTP.getDSRCount() <= 0 )
472 {
473 CPLDebug( "EnvisatDataset" , "No tiepoint covering "
474 "the measurement records." ) ;
475 return; /* No TPs - no extraction. */
476 }
477
478 /* check if TPs cover the whole range of MDSRs */
479 if(( arTP.getFirstOffset() < 0 )||( arTP.getLastOffset() < 0 ))
480 {
481 CPLDebug( "EnvisatDataset" , "The tiepoints do not cover "
482 "whole range of measurement records." ) ;
483 /* Not good but we can still extract some of the TPS, can we? */
484 }
485
486 /* Check TP record spacing */
487 if ((1+(arTP.getFirstOffset()+arTP.getLastOffset()+GetRasterYSize()-1)
488 / nLinesPerTiePoint ) != arTP.getDSRCount() )
489 {
490 CPLDebug( "EnvisatDataset", "Not enough tiepoints per column! "
491 "received=%d expected=%d", nTPPerColumn ,
492 1 + (arTP.getFirstOffset()+arTP.getLastOffset()+
493 GetRasterYSize()-1) / nLinesPerTiePoint ) ;
494 return; // That is far more serious - we risk misplacing TPs.
495 }
496
497 bool isBrowseProduct;
498 if ( 50*nTPPerLine + 13 == nDSRSize ) /* regular product */
499 {
500 isBrowseProduct = false ;
501 }
502 else if ( 8*nTPPerLine + 13 == nDSRSize ) /* browse product */
503 {
504 /* although BPs are rare there is no reason not to support them */
505 isBrowseProduct = true ;
506 }
507 else
508 {
509 CPLDebug( "EnvisatDataset", "Unexpected size of 'Tie points ADS' !"
510 " received=%d expected=%d or %d" , nDSRSize ,
511 50*nTPPerLine+13, 8*nTPPerLine+13 ) ;
512 return;
513 }
514
515 /* -------------------------------------------------------------------- */
516 /* Collect the first GCP set from each record. */
517 /* -------------------------------------------------------------------- */
518
519 GByte *pabyRecord = (GByte *) CPLMalloc(nDSRSize-13);
520
521 GUInt32 *tpLat = ((GUInt32*)pabyRecord) + nTPPerLine*0 ; /* latitude */
522 GUInt32 *tpLon = ((GUInt32*)pabyRecord) + nTPPerLine*1 ; /* longitude */
523 GUInt32 *tpLtc = ((GUInt32*)pabyRecord) + nTPPerLine*4 ; /* lat. DEM correction */
524 GUInt32 *tpLnc = ((GUInt32*)pabyRecord) + nTPPerLine*5 ; /* lon. DEM correction */
525
526 nGCPCount = 0;
527 pasGCPList = (GDAL_GCP *) CPLCalloc( sizeof(GDAL_GCP),
528 arTP.getDSRCount() * nTPPerLine );
529
530 for( int ir = 0 ; ir < arTP.getDSRCount() ; ir++ )
531 {
532 int iRecord = ir + arTP.getFirstIndex() ;
533
534 double dfGCPLine = 0.5 +
535 ( iRecord*nLinesPerTiePoint - arTP.getFirstOffset() ) ;
536
537 if( EnvisatFile_ReadDatasetRecordChunk( hEnvisatFile, nDatasetIndex,
538 iRecord , pabyRecord, 13 , -1 ) != SUCCESS )
539 continue;
540
541 for( int iGCP = 0; iGCP < nTPPerLine; iGCP++ )
542 {
543 GDALInitGCPs( 1, pasGCPList + nGCPCount );
544
545 CPLFree( pasGCPList[nGCPCount].pszId );
546
547 char szId[128];
548 snprintf( szId, sizeof(szId), "%d", nGCPCount+1 );
549 pasGCPList[nGCPCount].pszId = CPLStrdup( szId );
550
551 #define INT32(x) ((GInt32)CPL_MSBWORD32(x))
552
553 pasGCPList[nGCPCount].dfGCPX = 1e-6*INT32(tpLon[iGCP]) ;
554 pasGCPList[nGCPCount].dfGCPY = 1e-6*INT32(tpLat[iGCP]) ;
555 pasGCPList[nGCPCount].dfGCPZ = 0.0;
556
557 if( !isBrowseProduct ) /* add DEM corrections */
558 {
559 pasGCPList[nGCPCount].dfGCPX += 1e-6*INT32(tpLnc[iGCP]) ;
560 pasGCPList[nGCPCount].dfGCPY += 1e-6*INT32(tpLtc[iGCP]) ;
561 }
562
563 #undef INT32
564
565 pasGCPList[nGCPCount].dfGCPLine = dfGCPLine ;
566 pasGCPList[nGCPCount].dfGCPPixel = iGCP*nSamplesPerTiePoint + 0.5;
567
568 nGCPCount++;
569 }
570 }
571 CPLFree( pabyRecord );
572 }
573
574 /************************************************************************/
575 /* GetMetadataDomainList() */
576 /************************************************************************/
577
GetMetadataDomainList()578 char **EnvisatDataset::GetMetadataDomainList()
579 {
580 return CSLAddString(GDALDataset::GetMetadataDomainList(), "envisat-ds-*-*");
581 }
582
583 /************************************************************************/
584 /* GetMetadata() */
585 /************************************************************************/
586
GetMetadata(const char * pszDomain)587 char **EnvisatDataset::GetMetadata( const char * pszDomain )
588
589 {
590 if( pszDomain == nullptr || !STARTS_WITH_CI(pszDomain, "envisat-ds-") )
591 return GDALDataset::GetMetadata( pszDomain );
592
593 /* -------------------------------------------------------------------- */
594 /* Get the dataset name and record number. */
595 /* -------------------------------------------------------------------- */
596 char szDSName[128];
597 strncpy( szDSName, pszDomain+11, sizeof(szDSName) );
598 szDSName[sizeof(szDSName)-1] = 0;
599
600 int nRecord = -1;
601 for( int i = 0; i < (int) sizeof(szDSName)-1; i++ )
602 {
603 if( szDSName[i] == '-' )
604 {
605 szDSName[i] = '\0';
606 nRecord = atoi(szDSName+1);
607 break;
608 }
609 }
610
611 if( nRecord == -1 )
612 return nullptr;
613
614 /* -------------------------------------------------------------------- */
615 /* Get the dataset index and info. */
616 /* -------------------------------------------------------------------- */
617 int nDSIndex = EnvisatFile_GetDatasetIndex( hEnvisatFile, szDSName );
618 if( nDSIndex == -1 )
619 return nullptr;
620
621 int nDSRSize, nNumDSR;
622 EnvisatFile_GetDatasetInfo( hEnvisatFile, nDSIndex, nullptr, nullptr, nullptr,
623 nullptr, nullptr, &nNumDSR, &nDSRSize );
624
625 if( nDSRSize == -1 || nRecord < 0 || nRecord >= nNumDSR )
626 return nullptr;
627
628 /* -------------------------------------------------------------------- */
629 /* Read the requested record. */
630 /* -------------------------------------------------------------------- */
631 char *pszRecord = (char *) CPLMalloc(nDSRSize+1);
632
633 if( EnvisatFile_ReadDatasetRecord( hEnvisatFile, nDSIndex, nRecord,
634 pszRecord ) == FAILURE )
635 {
636 CPLFree( pszRecord );
637 return nullptr;
638 }
639
640 /* -------------------------------------------------------------------- */
641 /* Massage the data into a safe textual format. For now we */
642 /* just turn zero bytes into spaces. */
643 /* -------------------------------------------------------------------- */
644 CSLDestroy( papszTempMD );
645
646 char *pszEscapedRecord
647 = CPLEscapeString( pszRecord, nDSRSize, CPLES_BackslashQuotable );
648 papszTempMD = CSLSetNameValue( nullptr, "EscapedRecord", pszEscapedRecord );
649 CPLFree( pszEscapedRecord );
650
651 for( int i = 0; i < nDSRSize; i++ )
652 if( pszRecord[i] == '\0' )
653 pszRecord[i] = ' ';
654
655 papszTempMD = CSLSetNameValue( papszTempMD, "RawRecord", pszRecord );
656
657 CPLFree( pszRecord );
658
659 return papszTempMD;
660 }
661
662 /************************************************************************/
663 /* CollectDSDMetadata() */
664 /* */
665 /* Collect metadata based on any DSD entries with filenames */
666 /* associated. */
667 /************************************************************************/
668
CollectDSDMetadata()669 void EnvisatDataset::CollectDSDMetadata()
670
671 {
672 char *pszDSName, *pszFilename;
673
674 for( int iDSD = 0;
675 EnvisatFile_GetDatasetInfo( hEnvisatFile, iDSD, &pszDSName, nullptr,
676 &pszFilename, nullptr, nullptr, nullptr, nullptr ) == SUCCESS;
677 iDSD++ )
678 {
679 if( pszFilename == nullptr
680 || strlen(pszFilename) == 0
681 || STARTS_WITH_CI(pszFilename, "NOT USED")
682 || STARTS_WITH_CI(pszFilename, " "))
683 continue;
684
685 const int max_len = 128;
686 char szKey[max_len];
687
688 strcpy( szKey, "DS_");
689 strncat( szKey, pszDSName, max_len - strlen(szKey) - 1 );
690
691 // strip trailing spaces.
692 for( int i = static_cast<int>(strlen(szKey))-1; i && szKey[i] == ' '; i-- )
693 szKey[i] = '\0';
694
695 // convert spaces into underscores.
696 for( int i = 0; szKey[i] != '\0'; i++ )
697 {
698 if( szKey[i] == ' ' )
699 szKey[i] = '_';
700 }
701
702 strcat( szKey, "_NAME" );
703
704 char szTrimmedName[max_len];
705 strcpy( szTrimmedName, pszFilename );
706 for( int i = static_cast<int>(strlen(szTrimmedName))-1; i && szTrimmedName[i] == ' '; i--)
707 szTrimmedName[i] = '\0';
708
709 SetMetadataItem( szKey, szTrimmedName );
710 }
711 }
712
713 /************************************************************************/
714 /* CollectADSMetadata() */
715 /* */
716 /* Collect metadata from envisat ADS and GADS. */
717 /************************************************************************/
718
CollectADSMetadata()719 void EnvisatDataset::CollectADSMetadata()
720 {
721 int nNumDsr, nDSRSize;
722 const char *pszDSName, *pszDSType, *pszDSFilename;
723
724 const char *pszProduct
725 = EnvisatFile_GetKeyValueAsString( hEnvisatFile, MPH,
726 "PRODUCT", "" );
727
728 for( int nDSIndex = 0;
729 EnvisatFile_GetDatasetInfo( hEnvisatFile, nDSIndex,
730 (char **) &pszDSName,
731 (char **) &pszDSType,
732 (char **) &pszDSFilename,
733 nullptr, nullptr,
734 &nNumDsr, &nDSRSize ) == SUCCESS;
735 ++nDSIndex )
736 {
737 if( STARTS_WITH_CI(pszDSFilename, "NOT USED") || (nNumDsr <= 0) )
738 continue;
739 if( !EQUAL(pszDSType,"A") && !EQUAL(pszDSType,"G") )
740 continue;
741
742 for ( int nRecord = 0; nRecord < nNumDsr; ++nRecord )
743 {
744 char szPrefix[128];
745 strncpy( szPrefix, pszDSName, sizeof(szPrefix) - 1);
746 szPrefix[sizeof(szPrefix) - 1] = '\0';
747
748 // strip trailing spaces
749 for( int i = static_cast<int>(strlen(szPrefix))-1; i && szPrefix[i] == ' '; --i )
750 szPrefix[i] = '\0';
751
752 // convert spaces into underscores
753 for( int i = 0; szPrefix[i] != '\0'; i++ )
754 {
755 if( szPrefix[i] == ' ' )
756 szPrefix[i] = '_';
757 }
758
759 char *pszRecord = (char *) CPLMalloc(nDSRSize+1);
760
761 if( EnvisatFile_ReadDatasetRecord( hEnvisatFile, nDSIndex, nRecord,
762 pszRecord ) == FAILURE )
763 {
764 CPLFree( pszRecord );
765 return;
766 }
767
768 const EnvisatRecordDescr *pRecordDescr
769 = EnvisatFile_GetRecordDescriptor(pszProduct, pszDSName);
770 if (pRecordDescr)
771 {
772 const EnvisatFieldDescr *pField = pRecordDescr->pFields;
773 while ( pField && pField->szName )
774 {
775 char szValue[1024];
776 if ( CE_None == EnvisatFile_GetFieldAsString(pszRecord, nDSRSize,
777 pField, szValue, sizeof(szValue)) )
778 {
779 char szKey[256];
780 if (nNumDsr == 1)
781 snprintf( szKey, sizeof(szKey), "%s_%s", szPrefix, pField->szName);
782 else
783 // sprintf(szKey, "%s_%02d_%s", szPrefix, nRecord,
784 snprintf( szKey, sizeof(szKey), "%s_%d_%s", szPrefix, nRecord,
785 pField->szName);
786 SetMetadataItem(szKey, szValue, "RECORDS");
787 }
788 // silently ignore conversion errors
789
790 ++pField;
791 }
792 }
793 CPLFree( pszRecord );
794 }
795 }
796 }
797
798 /************************************************************************/
799 /* CollectMetadata() */
800 /* */
801 /* Collect metadata from the SPH or MPH header fields. */
802 /************************************************************************/
803
CollectMetadata(EnvisatFile_HeaderFlag eMPHOrSPH)804 void EnvisatDataset::CollectMetadata( EnvisatFile_HeaderFlag eMPHOrSPH )
805
806 {
807 for( int iKey = 0; true; iKey++ )
808 {
809 const char *pszKey
810 = EnvisatFile_GetKeyByIndex(hEnvisatFile, eMPHOrSPH, iKey);
811 if( pszKey == nullptr )
812 break;
813
814 const char *pszValue
815 = EnvisatFile_GetKeyValueAsString( hEnvisatFile, eMPHOrSPH,
816 pszKey, nullptr );
817
818 if( pszValue == nullptr )
819 continue;
820
821 // skip some uninteresting structural information.
822 if( EQUAL(pszKey,"TOT_SIZE")
823 || EQUAL(pszKey,"SPH_SIZE")
824 || EQUAL(pszKey,"NUM_DSD")
825 || EQUAL(pszKey,"DSD_SIZE")
826 || EQUAL(pszKey,"NUM_DATA_SETS") )
827 continue;
828
829 char szHeaderKey[128];
830 if( eMPHOrSPH == MPH )
831 snprintf( szHeaderKey, sizeof(szHeaderKey), "MPH_%s", pszKey );
832 else
833 snprintf( szHeaderKey, sizeof(szHeaderKey), "SPH_%s", pszKey );
834
835 SetMetadataItem( szHeaderKey, pszValue );
836 }
837 }
838
839 /************************************************************************/
840 /* Open() */
841 /************************************************************************/
842
Open(GDALOpenInfo * poOpenInfo)843 GDALDataset *EnvisatDataset::Open( GDALOpenInfo * poOpenInfo )
844
845 {
846 /* -------------------------------------------------------------------- */
847 /* Check the header. */
848 /* -------------------------------------------------------------------- */
849 if( poOpenInfo->nHeaderBytes < 8 || poOpenInfo->fpL == nullptr )
850 return nullptr;
851
852 if( !STARTS_WITH_CI((const char *) poOpenInfo->pabyHeader, "PRODUCT=") )
853 return nullptr;
854
855 /* -------------------------------------------------------------------- */
856 /* Try opening the dataset. */
857 /* -------------------------------------------------------------------- */
858 EnvisatFile *hEnvisatFile = nullptr;
859 if( EnvisatFile_Open( &hEnvisatFile, poOpenInfo->pszFilename, "r" )
860 == FAILURE )
861 return nullptr;
862
863 /* -------------------------------------------------------------------- */
864 /* Find a measurement type dataset to use as our reference */
865 /* raster band. */
866 /* -------------------------------------------------------------------- */
867 int dsr_size, num_dsr, ds_offset;
868 char *pszDSType = nullptr;
869
870 int ds_index = 0;
871 for( ; true; ds_index++ )
872 {
873 if( EnvisatFile_GetDatasetInfo( hEnvisatFile, ds_index,
874 nullptr, &pszDSType, nullptr,
875 &ds_offset, nullptr,
876 &num_dsr, &dsr_size ) == FAILURE )
877 {
878 CPLError( CE_Failure, CPLE_AppDefined,
879 "Unable to find \"MDS1\" measurement dataset in "
880 "Envisat file." );
881 EnvisatFile_Close( hEnvisatFile );
882 return nullptr;
883 }
884
885 /* Have we found what we are looking for? A Measurement ds. */
886 if( EQUAL(pszDSType,"M") )
887 break;
888 }
889
890 /* -------------------------------------------------------------------- */
891 /* Confirm the requested access is supported. */
892 /* -------------------------------------------------------------------- */
893 if( poOpenInfo->eAccess == GA_Update )
894 {
895 EnvisatFile_Close( hEnvisatFile );
896 CPLError( CE_Failure, CPLE_NotSupported,
897 "The ENVISAT driver does not support update access to existing"
898 " datasets.\n" );
899 return nullptr;
900 }
901 /* -------------------------------------------------------------------- */
902 /* Create a corresponding GDALDataset. */
903 /* -------------------------------------------------------------------- */
904 EnvisatDataset *poDS = new EnvisatDataset();
905
906 poDS->hEnvisatFile = hEnvisatFile;
907
908 /* -------------------------------------------------------------------- */
909 /* Setup image definition. */
910 /* -------------------------------------------------------------------- */
911 EnvisatFile_GetDatasetInfo( hEnvisatFile, ds_index,
912 nullptr, nullptr, nullptr, &ds_offset, nullptr,
913 &num_dsr, &dsr_size );
914
915 poDS->nRasterXSize = EnvisatFile_GetKeyValueAsInt( hEnvisatFile, SPH,
916 "LINE_LENGTH", 0 );
917 poDS->nRasterYSize = num_dsr;
918 poDS->eAccess = GA_ReadOnly;
919
920 const char *pszProduct
921 = EnvisatFile_GetKeyValueAsString( hEnvisatFile, MPH, "PRODUCT", "" );
922 const char *pszDataType
923 = EnvisatFile_GetKeyValueAsString( hEnvisatFile, SPH, "DATA_TYPE",
924 "" );
925 const char *pszSampleType
926 = EnvisatFile_GetKeyValueAsString( hEnvisatFile, SPH, "SAMPLE_TYPE",
927 "" );
928
929 GDALDataType eDataType;
930 if( EQUAL(pszDataType,"FLT32") && STARTS_WITH_CI(pszSampleType, "COMPLEX"))
931 eDataType = GDT_CFloat32;
932 else if( EQUAL(pszDataType,"FLT32") )
933 eDataType = GDT_Float32;
934 else if( EQUAL(pszDataType,"UWORD") )
935 eDataType = GDT_UInt16;
936 else if( EQUAL(pszDataType,"SWORD") && STARTS_WITH_CI(pszSampleType, "COMPLEX") )
937 eDataType = GDT_CInt16;
938 else if( EQUAL(pszDataType,"SWORD") )
939 eDataType = GDT_Int16;
940 else if( STARTS_WITH_CI(pszProduct,"ATS_TOA_1") )
941 {
942 /* all 16bit data, no line length provided */
943 eDataType = GDT_Int16;
944 poDS->nRasterXSize = (dsr_size - 20) / 2;
945 }
946 else if( poDS->nRasterXSize == 0 )
947 {
948 CPLError( CE_Warning, CPLE_AppDefined,
949 "Envisat product format not recognised. Assuming 8bit\n"
950 "with no per-record prefix data. Results may be useless!" );
951 eDataType = GDT_Byte;
952 poDS->nRasterXSize = dsr_size;
953 }
954 else
955 {
956 if( dsr_size >= 2 * poDS->nRasterXSize )
957 eDataType = GDT_UInt16;
958 else
959 eDataType = GDT_Byte;
960 }
961
962 const int bNative =
963 #ifdef CPL_LSB
964 FALSE
965 #else
966 TRUE
967 #endif
968 ;
969
970 int nPrefixBytes = dsr_size -
971 ((GDALGetDataTypeSize(eDataType) / 8) * poDS->nRasterXSize);
972
973 /* -------------------------------------------------------------------- */
974 /* Fail out if we didn't get non-zero sizes. */
975 /* -------------------------------------------------------------------- */
976 if( poDS->nRasterXSize < 1 || poDS->nRasterYSize < 1 )
977 {
978 CPLError( CE_Failure, CPLE_AppDefined,
979 "Unable to determine organization of dataset. It would\n"
980 "appear this is an Envisat dataset, but an unsupported\n"
981 "data product. Unable to utilize." );
982 delete poDS;
983 return nullptr;
984 }
985
986 poDS->fpImage = poOpenInfo->fpL;
987 poOpenInfo->fpL = nullptr;
988
989 /* -------------------------------------------------------------------- */
990 /* Try to collect GCPs. */
991 /* -------------------------------------------------------------------- */
992
993 /* -------------------------------------------------------------------- */
994 /* Scan for all datasets matching the reference dataset. */
995 /* -------------------------------------------------------------------- */
996 int num_dsr2, dsr_size2, iBand = 0;
997 const char *pszDSName = nullptr;
998 char szBandName[128];
999 bool bMiltiChannel;
1000
1001 for( ds_index = 0;
1002 EnvisatFile_GetDatasetInfo( hEnvisatFile, ds_index,
1003 (char **) &pszDSName, nullptr, nullptr,
1004 &ds_offset, nullptr,
1005 &num_dsr2, &dsr_size2 ) == SUCCESS;
1006 ds_index++ )
1007 {
1008 if( !EQUAL(pszDSType,"M") || num_dsr2 != num_dsr )
1009 continue;
1010
1011 if( STARTS_WITH_CI(pszProduct, "MER") && (pszProduct[8] == '2') &&
1012 ( (strstr(pszDSName, "MDS(16)") != nullptr) ||
1013 (strstr(pszDSName, "MDS(19)") != nullptr)) )
1014 bMiltiChannel = true;
1015 else
1016 bMiltiChannel = false;
1017
1018 if( (dsr_size2 == dsr_size) && !bMiltiChannel )
1019 {
1020 poDS->SetBand( iBand+1,
1021 new RawRasterBand( poDS, iBand+1, poDS->fpImage,
1022 ds_offset + nPrefixBytes,
1023 GDALGetDataTypeSize(eDataType) / 8,
1024 dsr_size,
1025 eDataType, bNative,
1026 RawRasterBand::OwnFP::NO ) );
1027 iBand++;
1028
1029 poDS->GetRasterBand(iBand)->SetDescription( pszDSName );
1030 }
1031 /* -------------------------------------------------------------------- */
1032 /* Handle MERIS Level 2 datasets with data type different from */
1033 /* the one declared in the SPH */
1034 /* -------------------------------------------------------------------- */
1035 else if( STARTS_WITH_CI(pszProduct, "MER") &&
1036 (strstr(pszDSName, "Flags") != nullptr) )
1037 {
1038 if (pszProduct[8] == '1')
1039 {
1040 // Flags
1041 poDS->SetBand( iBand+1,
1042 new RawRasterBand( poDS, iBand+1, poDS->fpImage,
1043 ds_offset + nPrefixBytes, 3,
1044 dsr_size, GDT_Byte, bNative,
1045 RawRasterBand::OwnFP::NO ) );
1046 iBand++;
1047
1048 poDS->GetRasterBand(iBand)->SetDescription( pszDSName );
1049
1050 // Detector indices
1051 poDS->SetBand( iBand+1,
1052 new RawRasterBand( poDS, iBand+1, poDS->fpImage,
1053 ds_offset + nPrefixBytes + 1,
1054 3, dsr_size, GDT_Int16,
1055 bNative,
1056 RawRasterBand::OwnFP::NO ) );
1057 iBand++;
1058
1059 const char *pszSuffix = strstr( pszDSName, "MDS" );
1060 if ( pszSuffix != nullptr)
1061 snprintf( szBandName, sizeof(szBandName), "Detector index %s", pszSuffix );
1062 else
1063 snprintf( szBandName, sizeof(szBandName), "%s", "Detector index" );
1064 poDS->GetRasterBand(iBand)->SetDescription( szBandName );
1065 }
1066 else if ( (pszProduct[8] == '2') &&
1067 (dsr_size2 >= 3 * poDS->nRasterXSize ) )
1068 {
1069 int nFlagPrefixBytes = dsr_size2 - 3 * poDS->nRasterXSize;
1070
1071 poDS->SetBand( iBand+1,
1072 new MerisL2FlagBand( poDS, iBand+1, poDS->fpImage,
1073 ds_offset, nFlagPrefixBytes ) );
1074 iBand++;
1075
1076 poDS->GetRasterBand(iBand)->SetDescription( pszDSName );
1077 }
1078 }
1079 else if( STARTS_WITH_CI(pszProduct, "MER") && (pszProduct[8] == '2') )
1080 {
1081 int nPrefixBytes2, nSubBands, nSubBandIdx, nSubBandOffset;
1082
1083 int nPixelSize = 1;
1084 GDALDataType eDataType2 = GDT_Byte;
1085
1086 nSubBands = dsr_size2 / poDS->nRasterXSize;
1087 if( (nSubBands < 1) || (nSubBands > 3) )
1088 nSubBands = 0;
1089
1090 nPrefixBytes2 = dsr_size2 -
1091 (nSubBands * nPixelSize * poDS->nRasterXSize);
1092
1093 for (nSubBandIdx = 0; nSubBandIdx < nSubBands; ++nSubBandIdx)
1094 {
1095 nSubBandOffset =
1096 ds_offset + nPrefixBytes2 + nSubBandIdx * nPixelSize;
1097 poDS->SetBand( iBand+1,
1098 new RawRasterBand( poDS, iBand+1, poDS->fpImage,
1099 nSubBandOffset,
1100 nPixelSize * nSubBands,
1101 dsr_size2, eDataType2, bNative,
1102 RawRasterBand::OwnFP::NO ) );
1103 iBand++;
1104
1105 if (nSubBands > 1)
1106 {
1107 snprintf( szBandName, sizeof(szBandName), "%s (%d)", pszDSName, nSubBandIdx );
1108 poDS->GetRasterBand(iBand)->SetDescription( szBandName );
1109 }
1110 else
1111 poDS->GetRasterBand(iBand)->SetDescription( pszDSName );
1112 }
1113 }
1114 }
1115
1116 /* -------------------------------------------------------------------- */
1117 /* Collect metadata. */
1118 /* -------------------------------------------------------------------- */
1119 poDS->CollectMetadata( MPH );
1120 poDS->CollectMetadata( SPH );
1121 poDS->CollectDSDMetadata();
1122 poDS->CollectADSMetadata();
1123
1124 if( STARTS_WITH_CI(pszProduct, "MER") )
1125 poDS->ScanForGCPs_MERIS();
1126 else
1127 poDS->ScanForGCPs_ASAR();
1128
1129 /* unwrap GCPs for products crossing date border */
1130 poDS->UnwrapGCPs();
1131
1132 /* -------------------------------------------------------------------- */
1133 /* Initialize any PAM information. */
1134 /* -------------------------------------------------------------------- */
1135 poDS->SetDescription( poOpenInfo->pszFilename );
1136 poDS->TryLoadXML();
1137
1138 /* -------------------------------------------------------------------- */
1139 /* Check for overviews. */
1140 /* -------------------------------------------------------------------- */
1141 poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
1142
1143 return poDS;
1144 }
1145
1146 /************************************************************************/
1147 /* GDALRegister_Envisat() */
1148 /************************************************************************/
1149
GDALRegister_Envisat()1150 void GDALRegister_Envisat()
1151
1152 {
1153 if( GDALGetDriverByName( "ESAT" ) != nullptr )
1154 return;
1155
1156 GDALDriver *poDriver = new GDALDriver();
1157
1158 poDriver->SetDescription( "ESAT" );
1159 poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
1160 poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "Envisat Image Format" );
1161 poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1162 "drivers/raster/esat.html" );
1163 poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "n1" );
1164 poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
1165
1166 poDriver->pfnOpen = EnvisatDataset::Open;
1167
1168 GetGDALDriverManager()->RegisterDriver( poDriver );
1169 }
1170