1 /******************************************************************************
2 *
3 * Project: ISIS Version 3 Driver
4 * Purpose: Implementation of ISIS3Dataset
5 * Author: Trent Hare (thare@usgs.gov)
6 * Frank Warmerdam (warmerdam@pobox.com)
7 * Even Rouault (even.rouault at spatialys.com)
8 *
9 * NOTE: Original code authored by Trent and placed in the public domain as
10 * per US government policy. I have (within my rights) appropriated it and
11 * placed it under the following license. This is not intended to diminish
12 * Trents contribution.
13 ******************************************************************************
14 * Copyright (c) 2007, Frank Warmerdam <warmerdam@pobox.com>
15 * Copyright (c) 2009-2010, Even Rouault <even.rouault at spatialys.com>
16 * Copyright (c) 2017 Hobu Inc
17 * Copyright (c) 2017, Dmitry Baryshnikov <polimax@mail.ru>
18 * Copyright (c) 2017, NextGIS <info@nextgis.com>
19 *
20 * Permission is hereby granted, free of charge, to any person obtaining a
21 * copy of this software and associated documentation files (the "Software"),
22 * to deal in the Software without restriction, including without limitation
23 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
24 * and/or sell copies of the Software, and to permit persons to whom the
25 * Software is furnished to do so, subject to the following conditions:
26 *
27 * The above copyright notice and this permission notice shall be included
28 * in all copies or substantial portions of the Software.
29 *
30 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
31 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
32 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
33 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
34 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
35 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
36 * DEALINGS IN THE SOFTWARE.
37 ****************************************************************************/
38
39 #include "cpl_json.h"
40 #include "cpl_string.h"
41 #include "cpl_time.h"
42 #include "cpl_vsi_error.h"
43 #include "gdal_frmts.h"
44 #include "gdal_proxy.h"
45 #include "nasakeywordhandler.h"
46 #include "ogrgeojsonreader.h"
47 #include "ogr_spatialref.h"
48 #include "rawdataset.h"
49 #include "vrtdataset.h"
50 #include "cpl_safemaths.hpp"
51
52 // For gethostname()
53 #ifdef _WIN32
54 #include <winsock2.h>
55 #else
56 #include <unistd.h>
57 #endif
58
59 #include <algorithm>
60 #include <map>
61 #include <utility> // pair
62 #include <vector>
63
64 // Constants coming from ISIS3 source code
65 // in isis/src/base/objs/SpecialPixel/SpecialPixel.h
66
67 //There are several types of special pixels
68 // * Isis::Null Pixel has no data available
69 // * Isis::Lis Pixel was saturated on the instrument
70 // * Isis::His Pixel was saturated on the instrument
71 // * Isis::Lrs Pixel was saturated during a computation
72 // * Isis::Hrs Pixel was saturated during a computation
73
74 // 1-byte special pixel values
75 const unsigned char NULL1 = 0;
76 const unsigned char LOW_REPR_SAT1 = 0;
77 const unsigned char LOW_INSTR_SAT1 = 0;
78 const unsigned char HIGH_INSTR_SAT1 = 255;
79 const unsigned char HIGH_REPR_SAT1 = 255;
80
81 // 2-byte unsigned special pixel values
82 const unsigned short NULLU2 = 0;
83 const unsigned short LOW_REPR_SATU2 = 1;
84 const unsigned short LOW_INSTR_SATU2 = 2;
85 const unsigned short HIGH_INSTR_SATU2 = 65534;
86 const unsigned short HIGH_REPR_SATU2 = 65535;
87
88 // 2-byte signed special pixel values
89 const short NULL2 = -32768;
90 const short LOW_REPR_SAT2 = -32767;
91 const short LOW_INSTR_SAT2 = -32766;
92 const short HIGH_INSTR_SAT2 = -32765;
93 const short HIGH_REPR_SAT2 = -32764;
94
95 // Define 4-byte special pixel values for IEEE floating point
96 const float NULL4 = -3.4028226550889045e+38f; // 0xFF7FFFFB;
97 const float LOW_REPR_SAT4 = -3.4028228579130005e+38f; // 0xFF7FFFFC;
98 const float LOW_INSTR_SAT4 = -3.4028230607370965e+38f; // 0xFF7FFFFD;
99 const float HIGH_INSTR_SAT4 = -3.4028232635611926e+38f; // 0xFF7FFFFE;
100 const float HIGH_REPR_SAT4 = -3.4028234663852886e+38f; // 0xFF7FFFFF;
101
102 // Must be large enough to hold an integer
103 static const char* const pszSTARTBYTE_PLACEHOLDER = "!*^STARTBYTE^*!";
104 // Must be large enough to hold an integer
105 static const char* const pszLABEL_BYTES_PLACEHOLDER = "!*^LABEL_BYTES^*!";
106 // Must be large enough to hold an integer
107 static const char* const pszHISTORY_STARTBYTE_PLACEHOLDER =
108 "!*^HISTORY_STARTBYTE^*!";
109
110 CPL_CVSID("$Id: isis3dataset.cpp c04e0ea92cfb521061e84b9c3ba75b0e30345ffd 2020-07-02 22:27:16 +0200 Even Rouault $")
111
112 /************************************************************************/
113 /* ==================================================================== */
114 /* ISISDataset */
115 /* ==================================================================== */
116 /************************************************************************/
117
118 class ISIS3Dataset final: public RawDataset
119 {
120 friend class ISIS3RawRasterBand;
121 friend class ISISTiledBand;
122 friend class ISIS3WrapperRasterBand;
123
124 class NonPixelSection
125 {
126 public:
127 CPLString osSrcFilename;
128 CPLString osDstFilename; // empty for same file
129 vsi_l_offset nSrcOffset;
130 vsi_l_offset nSize;
131 CPLString osPlaceHolder; // empty if not same file
132 };
133
134 VSILFILE *m_fpLabel; // label file (only used for writing)
135 VSILFILE *m_fpImage; // image data file. May be == fpLabel
136 GDALDataset *m_poExternalDS; // external dataset (GeoTIFF)
137 bool m_bGeoTIFFAsRegularExternal; // creation only
138 bool m_bGeoTIFFInitDone; // creation only
139
140 CPLString m_osExternalFilename;
141 bool m_bIsLabelWritten; // creation only
142
143 bool m_bIsTiled;
144 bool m_bInitToNodata; // creation only
145
146 NASAKeywordHandler m_oKeywords;
147
148 bool m_bGotTransform;
149 double m_adfGeoTransform[6];
150
151 bool m_bHasSrcNoData; // creation only
152 double m_dfSrcNoData; // creation only
153
154 OGRSpatialReference m_oSRS;
155
156 // creation only variables
157 CPLString m_osComment;
158 CPLString m_osLatitudeType;
159 CPLString m_osLongitudeDirection;
160 CPLString m_osTargetName;
161 bool m_bForce360;
162 bool m_bWriteBoundingDegrees;
163 CPLString m_osBoundingDegrees;
164
165 CPLJSONObject m_oJSonLabel;
166 CPLString m_osHistory; // creation only
167 bool m_bUseSrcLabel; // creation only
168 bool m_bUseSrcMapping; // creation only
169 bool m_bUseSrcHistory; // creation only
170 bool m_bAddGDALHistory; // creation only
171 CPLString m_osGDALHistory; // creation only
172 std::vector<NonPixelSection> m_aoNonPixelSections; // creation only
173 CPLJSONObject m_oSrcJSonLabel; // creation only
174 CPLStringList m_aosISIS3MD;
175 CPLStringList m_aosAdditionalFiles;
176 CPLString m_osFromFilename; // creation only
177
178 RawBinaryLayout m_sLayout{};
179
180 const char *GetKeyword( const char *pszPath,
181 const char *pszDefault = "");
182
183 double FixLong( double dfLong );
184 void BuildLabel();
185 void BuildHistory();
186 void WriteLabel();
187 void InvalidateLabel();
188
189 static CPLString SerializeAsPDL( const CPLJSONObject& oObj );
190 static void SerializeAsPDL( VSILFILE* fp, const CPLJSONObject& oObj,
191 int nDepth = 0 );
192
193 public:
194 ISIS3Dataset();
195 virtual ~ISIS3Dataset();
196
197 virtual int CloseDependentDatasets() override;
198
199 virtual CPLErr GetGeoTransform( double * padfTransform ) override;
200 virtual CPLErr SetGeoTransform( double * padfTransform ) override;
201
202 const OGRSpatialReference* GetSpatialRef() const override;
203 CPLErr SetSpatialRef(const OGRSpatialReference* poSRS) override;
204
205 virtual char **GetFileList() override;
206
207 virtual char **GetMetadataDomainList() override;
208 virtual char **GetMetadata( const char* pszDomain = "" ) override;
209 virtual CPLErr SetMetadata( char** papszMD, const char* pszDomain = "" )
210 override;
211
212 bool GetRawBinaryLayout(GDALDataset::RawBinaryLayout&) override;
213
214 static int Identify( GDALOpenInfo * );
215 static GDALDataset *Open( GDALOpenInfo * );
216 static GDALDataset *Create( const char * pszFilename,
217 int nXSize, int nYSize, int nBands,
218 GDALDataType eType, char ** papszOptions );
219 static GDALDataset* CreateCopy( const char *pszFilename,
220 GDALDataset *poSrcDS,
221 int bStrict,
222 char ** papszOptions,
223 GDALProgressFunc pfnProgress,
224 void * pProgressData );
225 };
226
227 /************************************************************************/
228 /* ==================================================================== */
229 /* ISISTiledBand */
230 /* ==================================================================== */
231 /************************************************************************/
232
233 class ISISTiledBand final: public GDALPamRasterBand
234 {
235 friend class ISIS3Dataset;
236
237 VSILFILE *m_fpVSIL;
238 GIntBig m_nFirstTileOffset;
239 GIntBig m_nXTileOffset;
240 GIntBig m_nYTileOffset;
241 int m_bNativeOrder;
242 bool m_bHasOffset;
243 bool m_bHasScale;
244 double m_dfOffset;
245 double m_dfScale;
246 double m_dfNoData;
247
248 public:
249
250 ISISTiledBand( GDALDataset *poDS, VSILFILE *fpVSIL,
251 int nBand, GDALDataType eDT,
252 int nTileXSize, int nTileYSize,
253 GIntBig nFirstTileOffset,
254 GIntBig nXTileOffset,
255 GIntBig nYTileOffset,
256 int bNativeOrder );
~ISISTiledBand()257 virtual ~ISISTiledBand() {}
258
259 virtual CPLErr IReadBlock( int, int, void * ) override;
260 virtual CPLErr IWriteBlock( int, int, void * ) override;
261
262 virtual double GetOffset( int *pbSuccess = nullptr ) override;
263 virtual double GetScale( int *pbSuccess = nullptr ) override;
264 virtual CPLErr SetOffset( double dfNewOffset ) override;
265 virtual CPLErr SetScale( double dfNewScale ) override;
266 virtual double GetNoDataValue( int *pbSuccess = nullptr ) override;
267 virtual CPLErr SetNoDataValue( double dfNewNoData ) override;
268
269 void SetMaskBand(GDALRasterBand* poMaskBand);
270 };
271
272 /************************************************************************/
273 /* ==================================================================== */
274 /* ISIS3RawRasterBand */
275 /* ==================================================================== */
276 /************************************************************************/
277
278 class ISIS3RawRasterBand final: public RawRasterBand
279 {
280 friend class ISIS3Dataset;
281
282 bool m_bHasOffset;
283 bool m_bHasScale;
284 double m_dfOffset;
285 double m_dfScale;
286 double m_dfNoData;
287
288 public:
289 ISIS3RawRasterBand( GDALDataset *l_poDS, int l_nBand,
290 VSILFILE * l_fpRaw,
291 vsi_l_offset l_nImgOffset,
292 int l_nPixelOffset,
293 int l_nLineOffset,
294 GDALDataType l_eDataType,
295 int l_bNativeOrder );
~ISIS3RawRasterBand()296 virtual ~ISIS3RawRasterBand() {}
297
298 virtual CPLErr IReadBlock( int, int, void * ) override;
299 virtual CPLErr IWriteBlock( int, int, void * ) override;
300
301 virtual CPLErr IRasterIO( GDALRWFlag, int, int, int, int,
302 void *, int, int, GDALDataType,
303 GSpacing nPixelSpace, GSpacing nLineSpace,
304 GDALRasterIOExtraArg* psExtraArg ) override;
305
306 virtual double GetOffset( int *pbSuccess = nullptr ) override;
307 virtual double GetScale( int *pbSuccess = nullptr ) override;
308 virtual CPLErr SetOffset( double dfNewOffset ) override;
309 virtual CPLErr SetScale( double dfNewScale ) override;
310 virtual double GetNoDataValue( int *pbSuccess = nullptr ) override;
311 virtual CPLErr SetNoDataValue( double dfNewNoData ) override;
312
313 void SetMaskBand(GDALRasterBand* poMaskBand);
314 };
315
316 /************************************************************************/
317 /* ==================================================================== */
318 /* ISIS3WrapperRasterBand */
319 /* */
320 /* proxy for bands stored in other formats. */
321 /* ==================================================================== */
322 /************************************************************************/
323 class ISIS3WrapperRasterBand final: public GDALProxyRasterBand
324 {
325 friend class ISIS3Dataset;
326
327 GDALRasterBand* m_poBaseBand;
328 bool m_bHasOffset;
329 bool m_bHasScale;
330 double m_dfOffset;
331 double m_dfScale;
332 double m_dfNoData;
333
334 protected:
RefUnderlyingRasterBand()335 virtual GDALRasterBand* RefUnderlyingRasterBand() override
336 { return m_poBaseBand; }
337
338 public:
339 explicit ISIS3WrapperRasterBand( GDALRasterBand* poBaseBandIn );
~ISIS3WrapperRasterBand()340 ~ISIS3WrapperRasterBand() {}
341
342 void InitFile();
343
344 virtual CPLErr Fill(double dfRealValue, double dfImaginaryValue = 0) override;
345 virtual CPLErr IWriteBlock( int, int, void * ) override;
346
347 virtual CPLErr IRasterIO( GDALRWFlag, int, int, int, int,
348 void *, int, int, GDALDataType,
349 GSpacing nPixelSpace, GSpacing nLineSpace,
350 GDALRasterIOExtraArg* psExtraArg ) override;
351
352 virtual double GetOffset( int *pbSuccess = nullptr ) override;
353 virtual double GetScale( int *pbSuccess = nullptr ) override;
354 virtual CPLErr SetOffset( double dfNewOffset ) override;
355 virtual CPLErr SetScale( double dfNewScale ) override;
356 virtual double GetNoDataValue( int *pbSuccess = nullptr ) override;
357 virtual CPLErr SetNoDataValue( double dfNewNoData ) override;
358
GetMaskFlags()359 int GetMaskFlags() override { return nMaskFlags; }
GetMaskBand()360 GDALRasterBand* GetMaskBand() override { return poMask; }
361 void SetMaskBand(GDALRasterBand* poMaskBand);
362 };
363
364 /************************************************************************/
365 /* ==================================================================== */
366 /* ISISMaskBand */
367 /* ==================================================================== */
368
369 class ISISMaskBand final: public GDALRasterBand
370 {
371 GDALRasterBand *m_poBaseBand;
372 void *m_pBuffer;
373
374 public:
375
376 explicit ISISMaskBand( GDALRasterBand* poBaseBand );
377 ~ISISMaskBand();
378
379 virtual CPLErr IReadBlock( int, int, void * ) override;
380
381 };
382
383 /************************************************************************/
384 /* ISISTiledBand() */
385 /************************************************************************/
386
ISISTiledBand(GDALDataset * poDSIn,VSILFILE * fpVSILIn,int nBandIn,GDALDataType eDT,int nTileXSize,int nTileYSize,GIntBig nFirstTileOffsetIn,GIntBig nXTileOffsetIn,GIntBig nYTileOffsetIn,int bNativeOrderIn)387 ISISTiledBand::ISISTiledBand( GDALDataset *poDSIn, VSILFILE *fpVSILIn,
388 int nBandIn, GDALDataType eDT,
389 int nTileXSize, int nTileYSize,
390 GIntBig nFirstTileOffsetIn,
391 GIntBig nXTileOffsetIn,
392 GIntBig nYTileOffsetIn,
393 int bNativeOrderIn ) :
394 m_fpVSIL(fpVSILIn),
395 m_nFirstTileOffset(0),
396 m_nXTileOffset(nXTileOffsetIn),
397 m_nYTileOffset(nYTileOffsetIn),
398 m_bNativeOrder(bNativeOrderIn),
399 m_bHasOffset(false),
400 m_bHasScale(false),
401 m_dfOffset(0.0),
402 m_dfScale(1.0),
403 m_dfNoData(0.0)
404 {
405 poDS = poDSIn;
406 nBand = nBandIn;
407 eDataType = eDT;
408 nBlockXSize = nTileXSize;
409 nBlockYSize = nTileYSize;
410 nRasterXSize = poDSIn->GetRasterXSize();
411 nRasterYSize = poDSIn->GetRasterYSize();
412
413 const int l_nBlocksPerRow = DIV_ROUND_UP(nRasterXSize, nBlockXSize);
414 const int l_nBlocksPerColumn = DIV_ROUND_UP(nRasterYSize, nBlockYSize);
415
416 if( m_nXTileOffset == 0 && m_nYTileOffset == 0 )
417 {
418 m_nXTileOffset =
419 static_cast<GIntBig>(GDALGetDataTypeSizeBytes(eDT)) *
420 nTileXSize;
421 if( m_nXTileOffset > GINTBIG_MAX / nTileYSize )
422 {
423 CPLError(CE_Failure, CPLE_AppDefined, "Integer overflow");
424 return;
425 }
426 m_nXTileOffset *= nTileYSize;
427
428 if( m_nXTileOffset > GINTBIG_MAX / l_nBlocksPerRow )
429 {
430 CPLError(CE_Failure, CPLE_AppDefined, "Integer overflow");
431 return;
432 }
433 m_nYTileOffset = m_nXTileOffset * l_nBlocksPerRow;
434 }
435
436 m_nFirstTileOffset = nFirstTileOffsetIn;
437 if( nBand > 1 )
438 {
439 if( m_nYTileOffset > GINTBIG_MAX / (nBand - 1) ||
440 (nBand-1) * m_nYTileOffset > GINTBIG_MAX / l_nBlocksPerColumn ||
441 m_nFirstTileOffset > GINTBIG_MAX - (nBand-1) * m_nYTileOffset * l_nBlocksPerColumn )
442 {
443 CPLError(CE_Failure, CPLE_AppDefined, "Integer overflow");
444 return;
445 }
446 m_nFirstTileOffset += (nBand-1) * m_nYTileOffset * l_nBlocksPerColumn;
447 }
448 }
449
450 /************************************************************************/
451 /* IReadBlock() */
452 /************************************************************************/
453
IReadBlock(int nXBlock,int nYBlock,void * pImage)454 CPLErr ISISTiledBand::IReadBlock( int nXBlock, int nYBlock, void *pImage )
455
456 {
457 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
458 if( poGDS->m_osExternalFilename.empty() )
459 {
460 if( !poGDS->m_bIsLabelWritten )
461 poGDS->WriteLabel();
462 }
463
464 const GIntBig nOffset = m_nFirstTileOffset +
465 nXBlock * m_nXTileOffset + nYBlock * m_nYTileOffset;
466 const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
467 const size_t nBlockSize = static_cast<size_t>(nDTSize)
468 * nBlockXSize * nBlockYSize;
469
470 if( VSIFSeekL( m_fpVSIL, nOffset, SEEK_SET ) != 0 )
471 {
472 CPLError( CE_Failure, CPLE_FileIO,
473 "Failed to seek to offset %d to read tile %d,%d.",
474 static_cast<int>( nOffset ), nXBlock, nYBlock );
475 return CE_Failure;
476 }
477
478 if( VSIFReadL( pImage, 1, nBlockSize, m_fpVSIL ) != nBlockSize )
479 {
480 CPLError( CE_Failure, CPLE_FileIO,
481 "Failed to read %d bytes for tile %d,%d.",
482 static_cast<int>( nBlockSize ), nXBlock, nYBlock );
483 return CE_Failure;
484 }
485
486 if( !m_bNativeOrder && eDataType != GDT_Byte )
487 GDALSwapWords( pImage, nDTSize,
488 nBlockXSize*nBlockYSize,
489 nDTSize );
490
491 return CE_None;
492 }
493
494 /************************************************************************/
495 /* RemapNoDataT() */
496 /************************************************************************/
497
RemapNoDataT(T * pBuffer,int nItems,T srcNoData,T dstNoData)498 template<class T> static void RemapNoDataT( T* pBuffer, int nItems,
499 T srcNoData, T dstNoData )
500 {
501 for( int i = 0; i < nItems; i++ )
502 {
503 if( pBuffer[i] == srcNoData )
504 pBuffer[i] = dstNoData;
505 }
506 }
507
508 /************************************************************************/
509 /* RemapNoData() */
510 /************************************************************************/
511
RemapNoData(GDALDataType eDataType,void * pBuffer,int nItems,double dfSrcNoData,double dfDstNoData)512 static void RemapNoData( GDALDataType eDataType,
513 void* pBuffer, int nItems, double dfSrcNoData,
514 double dfDstNoData )
515 {
516 if( eDataType == GDT_Byte )
517 {
518 RemapNoDataT( reinterpret_cast<GByte*>(pBuffer),
519 nItems,
520 static_cast<GByte>(dfSrcNoData),
521 static_cast<GByte>(dfDstNoData) );
522 }
523 else if( eDataType == GDT_UInt16 )
524 {
525 RemapNoDataT( reinterpret_cast<GUInt16*>(pBuffer),
526 nItems,
527 static_cast<GUInt16>(dfSrcNoData),
528 static_cast<GUInt16>(dfDstNoData) );
529 }
530 else if( eDataType == GDT_Int16)
531 {
532 RemapNoDataT( reinterpret_cast<GInt16*>(pBuffer),
533 nItems,
534 static_cast<GInt16>(dfSrcNoData),
535 static_cast<GInt16>(dfDstNoData) );
536 }
537 else
538 {
539 CPLAssert( eDataType == GDT_Float32 );
540 RemapNoDataT( reinterpret_cast<float*>(pBuffer),
541 nItems,
542 static_cast<float>(dfSrcNoData),
543 static_cast<float>(dfDstNoData) );
544 }
545 }
546
547 /**
548 * Get or create CPLJSONObject.
549 * @param oParent Parent CPLJSONObject.
550 * @param osKey Key name.
551 * @return CPLJSONObject class instance.
552 */
GetOrCreateJSONObject(CPLJSONObject & oParent,const std::string & osKey)553 static CPLJSONObject GetOrCreateJSONObject(CPLJSONObject &oParent,
554 const std::string &osKey)
555 {
556 CPLJSONObject oChild = oParent[osKey];
557 if( oChild.IsValid() && oChild.GetType() != CPLJSONObject::Type::Object )
558 {
559 oParent.Delete( osKey );
560 oChild.Deinit();
561 }
562
563 if( !oChild.IsValid() )
564 {
565 oChild = CPLJSONObject();
566 oParent.Add( osKey, oChild );
567 }
568 return oChild;
569 }
570
571 /************************************************************************/
572 /* IReadBlock() */
573 /************************************************************************/
574
IWriteBlock(int nXBlock,int nYBlock,void * pImage)575 CPLErr ISISTiledBand::IWriteBlock( int nXBlock, int nYBlock, void *pImage )
576
577 {
578 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
579 if( poGDS->m_osExternalFilename.empty() )
580 {
581 if( !poGDS->m_bIsLabelWritten )
582 poGDS->WriteLabel();
583 }
584
585 if( poGDS->m_bHasSrcNoData && poGDS->m_dfSrcNoData != m_dfNoData )
586 {
587 RemapNoData( eDataType, pImage, nBlockXSize * nBlockYSize,
588 poGDS->m_dfSrcNoData, m_dfNoData );
589 }
590
591 const GIntBig nOffset = m_nFirstTileOffset +
592 nXBlock * m_nXTileOffset + nYBlock * m_nYTileOffset;
593 const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
594 const size_t nBlockSize = static_cast<size_t>(nDTSize)
595 * nBlockXSize * nBlockYSize;
596
597 const int l_nBlocksPerRow = DIV_ROUND_UP(nRasterXSize, nBlockXSize);
598 const int l_nBlocksPerColumn = DIV_ROUND_UP(nRasterYSize, nBlockYSize);
599
600 // Pad partial blocks to nodata value
601 if( nXBlock == l_nBlocksPerRow - 1 && (nRasterXSize % nBlockXSize) != 0 )
602 {
603 GByte* pabyImage = static_cast<GByte*>(pImage);
604 int nXStart = nRasterXSize % nBlockXSize;
605 for( int iY = 0; iY < nBlockYSize; iY++ )
606 {
607 GDALCopyWords( &m_dfNoData, GDT_Float64, 0,
608 pabyImage + (iY * nBlockXSize + nXStart) * nDTSize,
609 eDataType, nDTSize,
610 nBlockXSize - nXStart );
611 }
612 }
613 if( nYBlock == l_nBlocksPerColumn - 1 &&
614 (nRasterYSize % nBlockYSize) != 0 )
615 {
616 GByte* pabyImage = static_cast<GByte*>(pImage);
617 for( int iY = nRasterYSize % nBlockYSize; iY < nBlockYSize; iY++ )
618 {
619 GDALCopyWords( &m_dfNoData, GDT_Float64, 0,
620 pabyImage + iY * nBlockXSize * nDTSize,
621 eDataType, nDTSize,
622 nBlockXSize );
623 }
624 }
625
626 if( VSIFSeekL( m_fpVSIL, nOffset, SEEK_SET ) != 0 )
627 {
628 CPLError( CE_Failure, CPLE_FileIO,
629 "Failed to seek to offset %d to read tile %d,%d.",
630 static_cast<int>( nOffset ), nXBlock, nYBlock );
631 return CE_Failure;
632 }
633
634 if( !m_bNativeOrder && eDataType != GDT_Byte )
635 GDALSwapWords( pImage, nDTSize,
636 nBlockXSize*nBlockYSize,
637 nDTSize );
638
639 if( VSIFWriteL( pImage, 1, nBlockSize, m_fpVSIL ) != nBlockSize )
640 {
641 CPLError( CE_Failure, CPLE_FileIO,
642 "Failed to write %d bytes for tile %d,%d.",
643 static_cast<int>( nBlockSize ), nXBlock, nYBlock );
644 return CE_Failure;
645 }
646
647 if( !m_bNativeOrder && eDataType != GDT_Byte )
648 GDALSwapWords( pImage, nDTSize,
649 nBlockXSize*nBlockYSize,
650 nDTSize );
651
652 return CE_None;
653 }
654
655 /************************************************************************/
656 /* SetMaskBand() */
657 /************************************************************************/
658
SetMaskBand(GDALRasterBand * poMaskBand)659 void ISISTiledBand::SetMaskBand(GDALRasterBand* poMaskBand)
660 {
661 bOwnMask = true;
662 poMask = poMaskBand;
663 nMaskFlags = 0;
664 }
665
666 /************************************************************************/
667 /* GetOffset() */
668 /************************************************************************/
669
GetOffset(int * pbSuccess)670 double ISISTiledBand::GetOffset( int *pbSuccess )
671 {
672 if( pbSuccess )
673 *pbSuccess = m_bHasOffset;
674 return m_dfOffset;
675 }
676
677 /************************************************************************/
678 /* GetScale() */
679 /************************************************************************/
680
GetScale(int * pbSuccess)681 double ISISTiledBand::GetScale( int *pbSuccess )
682 {
683 if( pbSuccess )
684 *pbSuccess = m_bHasScale;
685 return m_dfScale;
686 }
687
688 /************************************************************************/
689 /* SetOffset() */
690 /************************************************************************/
691
SetOffset(double dfNewOffset)692 CPLErr ISISTiledBand::SetOffset( double dfNewOffset )
693 {
694 m_dfOffset = dfNewOffset;
695 m_bHasOffset = true;
696 return CE_None;
697 }
698
699 /************************************************************************/
700 /* SetScale() */
701 /************************************************************************/
702
SetScale(double dfNewScale)703 CPLErr ISISTiledBand::SetScale( double dfNewScale )
704 {
705 m_dfScale = dfNewScale;
706 m_bHasScale = true;
707 return CE_None;
708 }
709
710 /************************************************************************/
711 /* GetNoDataValue() */
712 /************************************************************************/
713
GetNoDataValue(int * pbSuccess)714 double ISISTiledBand::GetNoDataValue( int *pbSuccess )
715 {
716 if( pbSuccess )
717 *pbSuccess = true;
718 return m_dfNoData;
719 }
720
721 /************************************************************************/
722 /* SetNoDataValue() */
723 /************************************************************************/
724
SetNoDataValue(double dfNewNoData)725 CPLErr ISISTiledBand::SetNoDataValue( double dfNewNoData )
726 {
727 m_dfNoData = dfNewNoData;
728 return CE_None;
729 }
730
731 /************************************************************************/
732 /* ISIS3RawRasterBand() */
733 /************************************************************************/
734
ISIS3RawRasterBand(GDALDataset * l_poDS,int l_nBand,VSILFILE * l_fpRaw,vsi_l_offset l_nImgOffset,int l_nPixelOffset,int l_nLineOffset,GDALDataType l_eDataType,int l_bNativeOrder)735 ISIS3RawRasterBand::ISIS3RawRasterBand( GDALDataset *l_poDS, int l_nBand,
736 VSILFILE * l_fpRaw,
737 vsi_l_offset l_nImgOffset,
738 int l_nPixelOffset,
739 int l_nLineOffset,
740 GDALDataType l_eDataType,
741 int l_bNativeOrder )
742 : RawRasterBand(l_poDS, l_nBand, l_fpRaw, l_nImgOffset, l_nPixelOffset,
743 l_nLineOffset,
744 l_eDataType, l_bNativeOrder, RawRasterBand::OwnFP::NO),
745 m_bHasOffset(false),
746 m_bHasScale(false),
747 m_dfOffset(0.0),
748 m_dfScale(1.0),
749 m_dfNoData(0.0)
750 {
751 }
752
753 /************************************************************************/
754 /* IReadBlock() */
755 /************************************************************************/
756
IReadBlock(int nXBlock,int nYBlock,void * pImage)757 CPLErr ISIS3RawRasterBand::IReadBlock( int nXBlock, int nYBlock, void *pImage )
758
759 {
760 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
761 if( poGDS->m_osExternalFilename.empty() )
762 {
763 if( !poGDS->m_bIsLabelWritten )
764 poGDS->WriteLabel();
765 }
766 return RawRasterBand::IReadBlock( nXBlock, nYBlock, pImage );
767 }
768
769 /************************************************************************/
770 /* IWriteBlock() */
771 /************************************************************************/
772
IWriteBlock(int nXBlock,int nYBlock,void * pImage)773 CPLErr ISIS3RawRasterBand::IWriteBlock( int nXBlock, int nYBlock,
774 void *pImage )
775
776 {
777 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
778 if( poGDS->m_osExternalFilename.empty() )
779 {
780 if( !poGDS->m_bIsLabelWritten )
781 poGDS->WriteLabel();
782 }
783
784 if( poGDS->m_bHasSrcNoData && poGDS->m_dfSrcNoData != m_dfNoData )
785 {
786 RemapNoData( eDataType, pImage, nBlockXSize * nBlockYSize,
787 poGDS->m_dfSrcNoData, m_dfNoData );
788 }
789
790 return RawRasterBand::IWriteBlock( nXBlock, nYBlock, pImage );
791 }
792
793 /************************************************************************/
794 /* IRasterIO() */
795 /************************************************************************/
796
IRasterIO(GDALRWFlag eRWFlag,int nXOff,int nYOff,int nXSize,int nYSize,void * pData,int nBufXSize,int nBufYSize,GDALDataType eBufType,GSpacing nPixelSpace,GSpacing nLineSpace,GDALRasterIOExtraArg * psExtraArg)797 CPLErr ISIS3RawRasterBand::IRasterIO( GDALRWFlag eRWFlag,
798 int nXOff, int nYOff, int nXSize, int nYSize,
799 void * pData, int nBufXSize, int nBufYSize,
800 GDALDataType eBufType,
801 GSpacing nPixelSpace, GSpacing nLineSpace,
802 GDALRasterIOExtraArg* psExtraArg )
803
804 {
805 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
806 if( poGDS->m_osExternalFilename.empty() )
807 {
808 if( !poGDS->m_bIsLabelWritten )
809 poGDS->WriteLabel();
810 }
811 if( eRWFlag == GF_Write &&
812 poGDS->m_bHasSrcNoData && poGDS->m_dfSrcNoData != m_dfNoData )
813 {
814 const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
815 if( eBufType == eDataType && nPixelSpace == nDTSize &&
816 nLineSpace == nPixelSpace * nBufXSize )
817 {
818 RemapNoData( eDataType, pData, nBufXSize * nBufYSize,
819 poGDS->m_dfSrcNoData, m_dfNoData );
820 }
821 else
822 {
823 const GByte* pabySrc = reinterpret_cast<GByte*>(pData);
824 GByte* pabyTemp = reinterpret_cast<GByte*>(
825 VSI_MALLOC3_VERBOSE(nDTSize, nBufXSize, nBufYSize));
826 for( int i = 0; i < nBufYSize; i++ )
827 {
828 GDALCopyWords( pabySrc + i * nLineSpace, eBufType,
829 static_cast<int>(nPixelSpace),
830 pabyTemp + i * nBufXSize * nDTSize,
831 eDataType, nDTSize,
832 nBufXSize );
833 }
834 RemapNoData( eDataType, pabyTemp, nBufXSize * nBufYSize,
835 poGDS->m_dfSrcNoData, m_dfNoData );
836 CPLErr eErr = RawRasterBand::IRasterIO( eRWFlag,
837 nXOff, nYOff, nXSize, nYSize,
838 pabyTemp, nBufXSize, nBufYSize,
839 eDataType,
840 nDTSize, nDTSize*nBufXSize,
841 psExtraArg );
842 VSIFree(pabyTemp);
843 return eErr;
844 }
845 }
846 return RawRasterBand::IRasterIO( eRWFlag,
847 nXOff, nYOff, nXSize, nYSize,
848 pData, nBufXSize, nBufYSize,
849 eBufType,
850 nPixelSpace, nLineSpace,
851 psExtraArg );
852 }
853
854 /************************************************************************/
855 /* SetMaskBand() */
856 /************************************************************************/
857
SetMaskBand(GDALRasterBand * poMaskBand)858 void ISIS3RawRasterBand::SetMaskBand(GDALRasterBand* poMaskBand)
859 {
860 bOwnMask = true;
861 poMask = poMaskBand;
862 nMaskFlags = 0;
863 }
864
865 /************************************************************************/
866 /* GetOffset() */
867 /************************************************************************/
868
GetOffset(int * pbSuccess)869 double ISIS3RawRasterBand::GetOffset( int *pbSuccess )
870 {
871 if( pbSuccess )
872 *pbSuccess = m_bHasOffset;
873 return m_dfOffset;
874 }
875
876 /************************************************************************/
877 /* GetScale() */
878 /************************************************************************/
879
GetScale(int * pbSuccess)880 double ISIS3RawRasterBand::GetScale( int *pbSuccess )
881 {
882 if( pbSuccess )
883 *pbSuccess = m_bHasScale;
884 return m_dfScale;
885 }
886
887 /************************************************************************/
888 /* SetOffset() */
889 /************************************************************************/
890
SetOffset(double dfNewOffset)891 CPLErr ISIS3RawRasterBand::SetOffset( double dfNewOffset )
892 {
893 m_dfOffset = dfNewOffset;
894 m_bHasOffset = true;
895 return CE_None;
896 }
897
898 /************************************************************************/
899 /* SetScale() */
900 /************************************************************************/
901
SetScale(double dfNewScale)902 CPLErr ISIS3RawRasterBand::SetScale( double dfNewScale )
903 {
904 m_dfScale = dfNewScale;
905 m_bHasScale = true;
906 return CE_None;
907 }
908
909 /************************************************************************/
910 /* GetNoDataValue() */
911 /************************************************************************/
912
GetNoDataValue(int * pbSuccess)913 double ISIS3RawRasterBand::GetNoDataValue( int *pbSuccess )
914 {
915 if( pbSuccess )
916 *pbSuccess = true;
917 return m_dfNoData;
918 }
919
920 /************************************************************************/
921 /* SetNoDataValue() */
922 /************************************************************************/
923
SetNoDataValue(double dfNewNoData)924 CPLErr ISIS3RawRasterBand::SetNoDataValue( double dfNewNoData )
925 {
926 m_dfNoData = dfNewNoData;
927 return CE_None;
928 }
929
930 /************************************************************************/
931 /* ISIS3WrapperRasterBand() */
932 /************************************************************************/
933
ISIS3WrapperRasterBand(GDALRasterBand * poBaseBandIn)934 ISIS3WrapperRasterBand::ISIS3WrapperRasterBand( GDALRasterBand* poBaseBandIn ) :
935 m_poBaseBand(poBaseBandIn),
936 m_bHasOffset(false),
937 m_bHasScale(false),
938 m_dfOffset(0.0),
939 m_dfScale(1.0),
940 m_dfNoData(0.0)
941 {
942 eDataType = m_poBaseBand->GetRasterDataType();
943 m_poBaseBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
944 }
945
946 /************************************************************************/
947 /* SetMaskBand() */
948 /************************************************************************/
949
SetMaskBand(GDALRasterBand * poMaskBand)950 void ISIS3WrapperRasterBand::SetMaskBand(GDALRasterBand* poMaskBand)
951 {
952 bOwnMask = true;
953 poMask = poMaskBand;
954 nMaskFlags = 0;
955 }
956
957 /************************************************************************/
958 /* GetOffset() */
959 /************************************************************************/
960
GetOffset(int * pbSuccess)961 double ISIS3WrapperRasterBand::GetOffset( int *pbSuccess )
962 {
963 if( pbSuccess )
964 *pbSuccess = m_bHasOffset;
965 return m_dfOffset;
966 }
967
968 /************************************************************************/
969 /* GetScale() */
970 /************************************************************************/
971
GetScale(int * pbSuccess)972 double ISIS3WrapperRasterBand::GetScale( int *pbSuccess )
973 {
974 if( pbSuccess )
975 *pbSuccess = m_bHasScale;
976 return m_dfScale;
977 }
978
979 /************************************************************************/
980 /* SetOffset() */
981 /************************************************************************/
982
SetOffset(double dfNewOffset)983 CPLErr ISIS3WrapperRasterBand::SetOffset( double dfNewOffset )
984 {
985 m_dfOffset = dfNewOffset;
986 m_bHasOffset = true;
987
988 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
989 if( poGDS->m_poExternalDS && eAccess == GA_Update )
990 poGDS->m_poExternalDS->GetRasterBand(nBand)->SetOffset(dfNewOffset);
991
992 return CE_None;
993 }
994
995 /************************************************************************/
996 /* SetScale() */
997 /************************************************************************/
998
SetScale(double dfNewScale)999 CPLErr ISIS3WrapperRasterBand::SetScale( double dfNewScale )
1000 {
1001 m_dfScale = dfNewScale;
1002 m_bHasScale = true;
1003
1004 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
1005 if( poGDS->m_poExternalDS && eAccess == GA_Update )
1006 poGDS->m_poExternalDS->GetRasterBand(nBand)->SetScale(dfNewScale);
1007
1008 return CE_None;
1009 }
1010
1011 /************************************************************************/
1012 /* GetNoDataValue() */
1013 /************************************************************************/
1014
GetNoDataValue(int * pbSuccess)1015 double ISIS3WrapperRasterBand::GetNoDataValue( int *pbSuccess )
1016 {
1017 if( pbSuccess )
1018 *pbSuccess = true;
1019 return m_dfNoData;
1020 }
1021
1022 /************************************************************************/
1023 /* SetNoDataValue() */
1024 /************************************************************************/
1025
SetNoDataValue(double dfNewNoData)1026 CPLErr ISIS3WrapperRasterBand::SetNoDataValue( double dfNewNoData )
1027 {
1028 m_dfNoData = dfNewNoData;
1029
1030 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
1031 if( poGDS->m_poExternalDS && eAccess == GA_Update )
1032 poGDS->m_poExternalDS->GetRasterBand(nBand)->SetNoDataValue(dfNewNoData);
1033
1034 return CE_None;
1035 }
1036
1037 /************************************************************************/
1038 /* InitFile() */
1039 /************************************************************************/
1040
InitFile()1041 void ISIS3WrapperRasterBand::InitFile()
1042 {
1043 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
1044 if( poGDS->m_bGeoTIFFAsRegularExternal && !poGDS->m_bGeoTIFFInitDone )
1045 {
1046 poGDS->m_bGeoTIFFInitDone = true;
1047
1048 const int nBands = poGDS->GetRasterCount();
1049 // We need to make sure that blocks are written in the right order
1050 for( int i = 0; i < nBands; i++ )
1051 {
1052 poGDS->m_poExternalDS->GetRasterBand(i+1)->Fill(m_dfNoData);
1053 }
1054 poGDS->m_poExternalDS->FlushCache();
1055
1056 // Check that blocks are effectively written in expected order.
1057 const int nBlockSizeBytes = nBlockXSize * nBlockYSize *
1058 GDALGetDataTypeSizeBytes(eDataType);
1059
1060 GIntBig nLastOffset = 0;
1061 bool bGoOn = true;
1062 const int l_nBlocksPerRow = DIV_ROUND_UP(nRasterXSize, nBlockXSize);
1063 const int l_nBlocksPerColumn = DIV_ROUND_UP(nRasterYSize, nBlockYSize);
1064 for( int i = 0; i < nBands && bGoOn; i++ )
1065 {
1066 for( int y = 0; y < l_nBlocksPerColumn && bGoOn; y++ )
1067 {
1068 for( int x = 0; x < l_nBlocksPerRow && bGoOn; x++ )
1069 {
1070 const char* pszBlockOffset = poGDS->m_poExternalDS->
1071 GetRasterBand(i+1)->GetMetadataItem(
1072 CPLSPrintf("BLOCK_OFFSET_%d_%d", x, y), "TIFF");
1073 if( pszBlockOffset )
1074 {
1075 GIntBig nOffset = CPLAtoGIntBig(pszBlockOffset);
1076 if( i != 0 || x != 0 || y != 0 )
1077 {
1078 if( nOffset != nLastOffset + nBlockSizeBytes )
1079 {
1080 CPLError(CE_Warning, CPLE_AppDefined,
1081 "Block %d,%d band %d not at expected "
1082 "offset",
1083 x, y, i+1);
1084 bGoOn = false;
1085 poGDS->m_bGeoTIFFAsRegularExternal = false;
1086 }
1087 }
1088 nLastOffset = nOffset;
1089 }
1090 else
1091 {
1092 CPLError(CE_Warning, CPLE_AppDefined,
1093 "Block %d,%d band %d not at expected "
1094 "offset",
1095 x, y, i+1);
1096 bGoOn = false;
1097 poGDS->m_bGeoTIFFAsRegularExternal = false;
1098 }
1099 }
1100 }
1101 }
1102 }
1103 }
1104
1105 /************************************************************************/
1106 /* Fill() */
1107 /************************************************************************/
1108
Fill(double dfRealValue,double dfImaginaryValue)1109 CPLErr ISIS3WrapperRasterBand::Fill(double dfRealValue, double dfImaginaryValue)
1110 {
1111 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
1112 if( poGDS->m_bHasSrcNoData && poGDS->m_dfSrcNoData == dfRealValue )
1113 {
1114 dfRealValue = m_dfNoData;
1115 }
1116 if( poGDS->m_bGeoTIFFAsRegularExternal && !poGDS->m_bGeoTIFFInitDone )
1117 {
1118 InitFile();
1119 }
1120
1121 return GDALProxyRasterBand::Fill( dfRealValue, dfImaginaryValue );
1122 }
1123
1124 /************************************************************************/
1125 /* IWriteBlock() */
1126 /************************************************************************/
1127
IWriteBlock(int nXBlock,int nYBlock,void * pImage)1128 CPLErr ISIS3WrapperRasterBand::IWriteBlock( int nXBlock, int nYBlock,
1129 void *pImage )
1130
1131 {
1132 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
1133 if( poGDS->m_bHasSrcNoData && poGDS->m_dfSrcNoData != m_dfNoData )
1134 {
1135 RemapNoData( eDataType, pImage, nBlockXSize * nBlockYSize,
1136 poGDS->m_dfSrcNoData, m_dfNoData );
1137 }
1138 if( poGDS->m_bGeoTIFFAsRegularExternal && !poGDS->m_bGeoTIFFInitDone )
1139 {
1140 InitFile();
1141 }
1142
1143 return GDALProxyRasterBand::IWriteBlock( nXBlock, nYBlock, pImage );
1144 }
1145
1146 /************************************************************************/
1147 /* IRasterIO() */
1148 /************************************************************************/
1149
IRasterIO(GDALRWFlag eRWFlag,int nXOff,int nYOff,int nXSize,int nYSize,void * pData,int nBufXSize,int nBufYSize,GDALDataType eBufType,GSpacing nPixelSpace,GSpacing nLineSpace,GDALRasterIOExtraArg * psExtraArg)1150 CPLErr ISIS3WrapperRasterBand::IRasterIO( GDALRWFlag eRWFlag,
1151 int nXOff, int nYOff, int nXSize, int nYSize,
1152 void * pData, int nBufXSize, int nBufYSize,
1153 GDALDataType eBufType,
1154 GSpacing nPixelSpace, GSpacing nLineSpace,
1155 GDALRasterIOExtraArg* psExtraArg )
1156
1157 {
1158 ISIS3Dataset* poGDS = reinterpret_cast<ISIS3Dataset*>(poDS);
1159 if( eRWFlag == GF_Write && poGDS->m_bGeoTIFFAsRegularExternal &&
1160 !poGDS->m_bGeoTIFFInitDone )
1161 {
1162 InitFile();
1163 }
1164 if( eRWFlag == GF_Write &&
1165 poGDS->m_bHasSrcNoData && poGDS->m_dfSrcNoData != m_dfNoData )
1166 {
1167 const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
1168 if( eBufType == eDataType && nPixelSpace == nDTSize &&
1169 nLineSpace == nPixelSpace * nBufXSize )
1170 {
1171 RemapNoData( eDataType, pData, nBufXSize * nBufYSize,
1172 poGDS->m_dfSrcNoData, m_dfNoData );
1173 }
1174 else
1175 {
1176 const GByte* pabySrc = reinterpret_cast<GByte*>(pData);
1177 GByte* pabyTemp = reinterpret_cast<GByte*>(
1178 VSI_MALLOC3_VERBOSE(nDTSize, nBufXSize, nBufYSize));
1179 for( int i = 0; i < nBufYSize; i++ )
1180 {
1181 GDALCopyWords( pabySrc + i * nLineSpace, eBufType,
1182 static_cast<int>(nPixelSpace),
1183 pabyTemp + i * nBufXSize * nDTSize,
1184 eDataType, nDTSize,
1185 nBufXSize );
1186 }
1187 RemapNoData( eDataType, pabyTemp, nBufXSize * nBufYSize,
1188 poGDS->m_dfSrcNoData, m_dfNoData );
1189 CPLErr eErr = GDALProxyRasterBand::IRasterIO( eRWFlag,
1190 nXOff, nYOff, nXSize, nYSize,
1191 pabyTemp, nBufXSize, nBufYSize,
1192 eDataType,
1193 nDTSize, nDTSize*nBufXSize,
1194 psExtraArg );
1195 VSIFree(pabyTemp);
1196 return eErr;
1197 }
1198 }
1199 return GDALProxyRasterBand::IRasterIO( eRWFlag,
1200 nXOff, nYOff, nXSize, nYSize,
1201 pData, nBufXSize, nBufYSize,
1202 eBufType,
1203 nPixelSpace, nLineSpace,
1204 psExtraArg );
1205 }
1206
1207 /************************************************************************/
1208 /* ISISMaskBand() */
1209 /************************************************************************/
1210
ISISMaskBand(GDALRasterBand * poBaseBand)1211 ISISMaskBand::ISISMaskBand( GDALRasterBand* poBaseBand )
1212 : m_poBaseBand(poBaseBand)
1213 , m_pBuffer(nullptr)
1214 {
1215 eDataType = GDT_Byte;
1216 poBaseBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
1217 nRasterXSize = poBaseBand->GetXSize();
1218 nRasterYSize = poBaseBand->GetYSize();
1219 }
1220
1221 /************************************************************************/
1222 /* ~ISISMaskBand() */
1223 /************************************************************************/
1224
~ISISMaskBand()1225 ISISMaskBand::~ISISMaskBand()
1226 {
1227 VSIFree(m_pBuffer);
1228 }
1229
1230 /************************************************************************/
1231 /* FillMask() */
1232 /************************************************************************/
1233
1234 template<class T>
FillMask(void * pvBuffer,GByte * pabyDst,int nReqXSize,int nReqYSize,int nBlockXSize,T NULL_VAL,T LOW_REPR_SAT,T LOW_INSTR_SAT,T HIGH_INSTR_SAT,T HIGH_REPR_SAT)1235 static void FillMask (void* pvBuffer,
1236 GByte* pabyDst,
1237 int nReqXSize, int nReqYSize,
1238 int nBlockXSize,
1239 T NULL_VAL, T LOW_REPR_SAT, T LOW_INSTR_SAT,
1240 T HIGH_INSTR_SAT, T HIGH_REPR_SAT)
1241 {
1242 const T* pSrc = static_cast<T*>(pvBuffer);
1243 for( int y = 0; y < nReqYSize; y++ )
1244 {
1245 for( int x = 0; x < nReqXSize; x++ )
1246 {
1247 const T nSrc = pSrc[y * nBlockXSize + x];
1248 if( nSrc == NULL_VAL ||
1249 nSrc == LOW_REPR_SAT ||
1250 nSrc == LOW_INSTR_SAT ||
1251 nSrc == HIGH_INSTR_SAT ||
1252 nSrc == HIGH_REPR_SAT )
1253 {
1254 pabyDst[y * nBlockXSize + x] = 0;
1255 }
1256 else
1257 {
1258 pabyDst[y * nBlockXSize + x] = 255;
1259 }
1260 }
1261 }
1262 }
1263
1264 /************************************************************************/
1265 /* IReadBlock() */
1266 /************************************************************************/
1267
IReadBlock(int nXBlock,int nYBlock,void * pImage)1268 CPLErr ISISMaskBand::IReadBlock( int nXBlock, int nYBlock, void *pImage )
1269
1270 {
1271 const GDALDataType eSrcDT = m_poBaseBand->GetRasterDataType();
1272 const int nSrcDTSize = GDALGetDataTypeSizeBytes(eSrcDT);
1273 if( m_pBuffer == nullptr )
1274 {
1275 m_pBuffer = VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, nSrcDTSize);
1276 if( m_pBuffer == nullptr )
1277 return CE_Failure;
1278 }
1279
1280 int nXOff = nXBlock * nBlockXSize;
1281 int nReqXSize = nBlockXSize;
1282 if( nXOff + nReqXSize > nRasterXSize )
1283 nReqXSize = nRasterXSize - nXOff;
1284 int nYOff = nYBlock * nBlockYSize;
1285 int nReqYSize = nBlockYSize;
1286 if( nYOff + nReqYSize > nRasterYSize )
1287 nReqYSize = nRasterYSize - nYOff;
1288
1289 if( m_poBaseBand->RasterIO( GF_Read,
1290 nXOff, nYOff, nReqXSize, nReqYSize,
1291 m_pBuffer,
1292 nReqXSize, nReqYSize,
1293 eSrcDT,
1294 nSrcDTSize,
1295 nSrcDTSize * nBlockXSize,
1296 nullptr ) != CE_None )
1297 {
1298 return CE_Failure;
1299 }
1300
1301 GByte* pabyDst = static_cast<GByte*>(pImage);
1302 if( eSrcDT == GDT_Byte )
1303 {
1304 FillMask<GByte>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
1305 NULL1, LOW_REPR_SAT1, LOW_INSTR_SAT1,
1306 HIGH_INSTR_SAT1, HIGH_REPR_SAT1);
1307 }
1308 else if( eSrcDT == GDT_UInt16 )
1309 {
1310 FillMask<GUInt16>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
1311 NULLU2, LOW_REPR_SATU2, LOW_INSTR_SATU2,
1312 HIGH_INSTR_SATU2, HIGH_REPR_SATU2);
1313 }
1314 else if( eSrcDT == GDT_Int16 )
1315 {
1316 FillMask<GInt16>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
1317 NULL2, LOW_REPR_SAT2, LOW_INSTR_SAT2,
1318 HIGH_INSTR_SAT2, HIGH_REPR_SAT2);
1319 }
1320 else
1321 {
1322 CPLAssert( eSrcDT == GDT_Float32 );
1323 FillMask<float>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
1324 NULL4, LOW_REPR_SAT4, LOW_INSTR_SAT4,
1325 HIGH_INSTR_SAT4, HIGH_REPR_SAT4);
1326 }
1327
1328 return CE_None;
1329 }
1330
1331 /************************************************************************/
1332 /* ISIS3Dataset() */
1333 /************************************************************************/
1334
ISIS3Dataset()1335 ISIS3Dataset::ISIS3Dataset() :
1336 m_fpLabel(nullptr),
1337 m_fpImage(nullptr),
1338 m_poExternalDS(nullptr),
1339 m_bGeoTIFFAsRegularExternal(false),
1340 m_bGeoTIFFInitDone(true),
1341 m_bIsLabelWritten(true),
1342 m_bIsTiled(false),
1343 m_bInitToNodata(false),
1344 m_bGotTransform(false),
1345 m_bHasSrcNoData(false),
1346 m_dfSrcNoData(0.0),
1347 m_bForce360(false),
1348 m_bWriteBoundingDegrees(true),
1349 m_bUseSrcLabel(true),
1350 m_bUseSrcMapping(false),
1351 m_bUseSrcHistory(true),
1352 m_bAddGDALHistory(true)
1353 {
1354 m_oKeywords.SetStripSurroundingQuotes(true);
1355 m_adfGeoTransform[0] = 0.0;
1356 m_adfGeoTransform[1] = 1.0;
1357 m_adfGeoTransform[2] = 0.0;
1358 m_adfGeoTransform[3] = 0.0;
1359 m_adfGeoTransform[4] = 0.0;
1360 m_adfGeoTransform[5] = 1.0;
1361
1362 // Deinit JSON objects
1363 m_oJSonLabel.Deinit();
1364 m_oSrcJSonLabel.Deinit();
1365 }
1366
1367 /************************************************************************/
1368 /* ~ISIS3Dataset() */
1369 /************************************************************************/
1370
~ISIS3Dataset()1371 ISIS3Dataset::~ISIS3Dataset()
1372
1373 {
1374 if( !m_bIsLabelWritten )
1375 WriteLabel();
1376 if( m_poExternalDS && m_bGeoTIFFAsRegularExternal && !m_bGeoTIFFInitDone )
1377 {
1378 reinterpret_cast<ISIS3WrapperRasterBand*>(GetRasterBand(1))->
1379 InitFile();
1380 }
1381 ISIS3Dataset::FlushCache();
1382 if( m_fpLabel != nullptr )
1383 VSIFCloseL( m_fpLabel );
1384 if( m_fpImage != nullptr && m_fpImage != m_fpLabel )
1385 VSIFCloseL( m_fpImage );
1386
1387 ISIS3Dataset::CloseDependentDatasets();
1388 }
1389
1390 /************************************************************************/
1391 /* CloseDependentDatasets() */
1392 /************************************************************************/
1393
CloseDependentDatasets()1394 int ISIS3Dataset::CloseDependentDatasets()
1395 {
1396 int bHasDroppedRef = GDALPamDataset::CloseDependentDatasets();
1397
1398 if( m_poExternalDS )
1399 {
1400 bHasDroppedRef = FALSE;
1401 delete m_poExternalDS;
1402 m_poExternalDS = nullptr;
1403 }
1404
1405 for( int iBand = 0; iBand < nBands; iBand++ )
1406 {
1407 delete papoBands[iBand];
1408 }
1409 nBands = 0;
1410
1411 return bHasDroppedRef;
1412 }
1413
1414 /************************************************************************/
1415 /* GetFileList() */
1416 /************************************************************************/
1417
GetFileList()1418 char **ISIS3Dataset::GetFileList()
1419
1420 {
1421 char **papszFileList = GDALPamDataset::GetFileList();
1422
1423 if( !m_osExternalFilename.empty() )
1424 papszFileList = CSLAddString( papszFileList, m_osExternalFilename );
1425 for( int i = 0; i < m_aosAdditionalFiles.Count(); ++i )
1426 {
1427 if( CSLFindString(papszFileList, m_aosAdditionalFiles[i]) < 0 )
1428 {
1429 papszFileList = CSLAddString( papszFileList,
1430 m_aosAdditionalFiles[i] );
1431 }
1432 }
1433
1434 return papszFileList;
1435 }
1436
1437 /************************************************************************/
1438 /* GetSpatialRef() */
1439 /************************************************************************/
1440
GetSpatialRef() const1441 const OGRSpatialReference* ISIS3Dataset::GetSpatialRef() const
1442
1443 {
1444 if( !m_oSRS.IsEmpty() )
1445 return &m_oSRS;
1446
1447 return GDALPamDataset::GetSpatialRef();
1448 }
1449
1450 /************************************************************************/
1451 /* SetSpatialRef() */
1452 /************************************************************************/
1453
SetSpatialRef(const OGRSpatialReference * poSRS)1454 CPLErr ISIS3Dataset::SetSpatialRef( const OGRSpatialReference* poSRS )
1455 {
1456 if( eAccess == GA_ReadOnly )
1457 return GDALPamDataset::SetSpatialRef( poSRS );
1458 if( poSRS )
1459 m_oSRS = *poSRS;
1460 else
1461 m_oSRS.Clear();
1462 if( m_poExternalDS )
1463 m_poExternalDS->SetSpatialRef(poSRS);
1464 InvalidateLabel();
1465 return CE_None;
1466 }
1467
1468 /************************************************************************/
1469 /* GetGeoTransform() */
1470 /************************************************************************/
1471
GetGeoTransform(double * padfTransform)1472 CPLErr ISIS3Dataset::GetGeoTransform( double * padfTransform )
1473
1474 {
1475 if( m_bGotTransform )
1476 {
1477 memcpy( padfTransform, m_adfGeoTransform, sizeof(double) * 6 );
1478 return CE_None;
1479 }
1480
1481 return GDALPamDataset::GetGeoTransform( padfTransform );
1482 }
1483
1484 /************************************************************************/
1485 /* SetGeoTransform() */
1486 /************************************************************************/
1487
SetGeoTransform(double * padfTransform)1488 CPLErr ISIS3Dataset::SetGeoTransform( double * padfTransform )
1489
1490 {
1491 if( eAccess == GA_ReadOnly )
1492 return GDALPamDataset::SetGeoTransform( padfTransform );
1493 if( padfTransform[1] <= 0.0 || padfTransform[1] != -padfTransform[5] ||
1494 padfTransform[2] != 0.0 || padfTransform[4] != 0.0 )
1495 {
1496 CPLError(CE_Failure, CPLE_NotSupported,
1497 "Only north-up geotransform with square pixels supported");
1498 return CE_Failure;
1499 }
1500 m_bGotTransform = true;
1501 memcpy( m_adfGeoTransform, padfTransform, sizeof(double) * 6 );
1502 if( m_poExternalDS )
1503 m_poExternalDS->SetGeoTransform(padfTransform);
1504 InvalidateLabel();
1505 return CE_None;
1506 }
1507
1508 /************************************************************************/
1509 /* GetMetadataDomainList() */
1510 /************************************************************************/
1511
GetMetadataDomainList()1512 char **ISIS3Dataset::GetMetadataDomainList()
1513 {
1514 return BuildMetadataDomainList(
1515 nullptr, FALSE, "", "json:ISIS3", nullptr);
1516 }
1517
1518 /************************************************************************/
1519 /* GetMetadata() */
1520 /************************************************************************/
1521
GetMetadata(const char * pszDomain)1522 char **ISIS3Dataset::GetMetadata( const char* pszDomain )
1523 {
1524 if( pszDomain != nullptr && EQUAL( pszDomain, "json:ISIS3" ) )
1525 {
1526 if( m_aosISIS3MD.empty() )
1527 {
1528 if( eAccess == GA_Update && !m_oJSonLabel.IsValid() )
1529 {
1530 BuildLabel();
1531 }
1532 CPLAssert( m_oJSonLabel.IsValid() );
1533 const CPLString osJson = m_oJSonLabel.Format(CPLJSONObject::PrettyFormat::Pretty);
1534 m_aosISIS3MD.InsertString(0, osJson.c_str());
1535 }
1536 return m_aosISIS3MD.List();
1537 }
1538 return GDALPamDataset::GetMetadata(pszDomain);
1539 }
1540
1541 /************************************************************************/
1542 /* InvalidateLabel() */
1543 /************************************************************************/
1544
InvalidateLabel()1545 void ISIS3Dataset::InvalidateLabel()
1546 {
1547 m_oJSonLabel.Deinit();
1548 m_aosISIS3MD.Clear();
1549 }
1550
1551 /************************************************************************/
1552 /* SetMetadata() */
1553 /************************************************************************/
1554
SetMetadata(char ** papszMD,const char * pszDomain)1555 CPLErr ISIS3Dataset::SetMetadata( char** papszMD, const char* pszDomain )
1556 {
1557 if( m_bUseSrcLabel && eAccess == GA_Update && pszDomain != nullptr &&
1558 EQUAL( pszDomain, "json:ISIS3" ) )
1559 {
1560 m_oSrcJSonLabel.Deinit();
1561 InvalidateLabel();
1562 if( papszMD != nullptr && papszMD[0] != nullptr )
1563 {
1564 CPLJSONDocument oJSONDocument;
1565 const GByte *pabyData = reinterpret_cast<const GByte *>(papszMD[0]);
1566 if( !oJSONDocument.LoadMemory( pabyData ) )
1567 {
1568 return CE_Failure;
1569 }
1570
1571 m_oSrcJSonLabel = oJSONDocument.GetRoot();
1572 if( !m_oSrcJSonLabel.IsValid() )
1573 {
1574 return CE_Failure;
1575 }
1576 }
1577 return CE_None;
1578 }
1579 return GDALPamDataset::SetMetadata(papszMD, pszDomain);
1580 }
1581
1582 /************************************************************************/
1583 /* Identify() */
1584 /************************************************************************/
Identify(GDALOpenInfo * poOpenInfo)1585 int ISIS3Dataset::Identify( GDALOpenInfo * poOpenInfo )
1586
1587 {
1588 if( poOpenInfo->fpL != nullptr &&
1589 poOpenInfo->pabyHeader != nullptr &&
1590 strstr((const char *)poOpenInfo->pabyHeader,"IsisCube") != nullptr )
1591 return TRUE;
1592
1593 return FALSE;
1594 }
1595
1596 /************************************************************************/
1597 /* GetRawBinaryLayout() */
1598 /************************************************************************/
1599
GetRawBinaryLayout(GDALDataset::RawBinaryLayout & sLayout)1600 bool ISIS3Dataset::GetRawBinaryLayout(GDALDataset::RawBinaryLayout& sLayout)
1601 {
1602 if( m_sLayout.osRawFilename.empty() )
1603 return false;
1604 sLayout = m_sLayout;
1605 return true;
1606 }
1607
1608 /************************************************************************/
1609 /* GetValueAndUnits() */
1610 /************************************************************************/
1611
GetValueAndUnits(const CPLJSONObject & obj,std::vector<double> & adfValues,std::vector<std::string> & aosUnits,int nExpectedVals)1612 static void GetValueAndUnits(const CPLJSONObject& obj,
1613 std::vector<double>& adfValues,
1614 std::vector<std::string>& aosUnits,
1615 int nExpectedVals)
1616 {
1617 if( obj.GetType() == CPLJSONObject::Type::Integer ||
1618 obj.GetType() == CPLJSONObject::Type::Double )
1619 {
1620 adfValues.push_back(obj.ToDouble());
1621 }
1622 else if( obj.GetType() == CPLJSONObject::Type::Object )
1623 {
1624 auto oValue = obj.GetObj("value");
1625 auto oUnit = obj.GetObj("unit");
1626 if( oValue.IsValid() &&
1627 (oValue.GetType() == CPLJSONObject::Type::Integer ||
1628 oValue.GetType() == CPLJSONObject::Type::Double ||
1629 oValue.GetType() == CPLJSONObject::Type::Array) &&
1630 oUnit.IsValid() && oUnit.GetType() == CPLJSONObject::Type::String )
1631 {
1632 if( oValue.GetType() == CPLJSONObject::Type::Array )
1633 {
1634 GetValueAndUnits(oValue, adfValues, aosUnits, nExpectedVals);
1635 }
1636 else
1637 {
1638 adfValues.push_back(oValue.ToDouble());
1639 }
1640 aosUnits.push_back(oUnit.ToString());
1641 }
1642 }
1643 else if( obj.GetType() == CPLJSONObject::Type::Array )
1644 {
1645 auto oArray = obj.ToArray();
1646 if( oArray.Size() == nExpectedVals )
1647 {
1648 for( int i = 0; i < nExpectedVals; i++ )
1649 {
1650 if( oArray[i].GetType() == CPLJSONObject::Type::Integer ||
1651 oArray[i].GetType() == CPLJSONObject::Type::Double )
1652 {
1653 adfValues.push_back(oArray[i].ToDouble());
1654 }
1655 else
1656 {
1657 adfValues.clear();
1658 return;
1659 }
1660 }
1661 }
1662 }
1663 }
1664
1665 /************************************************************************/
1666 /* Open() */
1667 /************************************************************************/
1668
Open(GDALOpenInfo * poOpenInfo)1669 GDALDataset *ISIS3Dataset::Open( GDALOpenInfo * poOpenInfo )
1670
1671 {
1672 /* -------------------------------------------------------------------- */
1673 /* Does this look like a CUBE dataset? */
1674 /* -------------------------------------------------------------------- */
1675 if( !Identify( poOpenInfo ) )
1676 return nullptr;
1677
1678 /* -------------------------------------------------------------------- */
1679 /* Open the file using the large file API. */
1680 /* -------------------------------------------------------------------- */
1681 ISIS3Dataset *poDS = new ISIS3Dataset();
1682
1683 if( ! poDS->m_oKeywords.Ingest( poOpenInfo->fpL, 0 ) )
1684 {
1685 VSIFCloseL( poOpenInfo->fpL );
1686 poOpenInfo->fpL = nullptr;
1687 delete poDS;
1688 return nullptr;
1689 }
1690 poDS->m_oJSonLabel = poDS->m_oKeywords.GetJsonObject();
1691 poDS->m_oJSonLabel.Add( "_filename", poOpenInfo->pszFilename );
1692
1693 // Find additional files from the label
1694 for( const CPLJSONObject& oObj : poDS->m_oJSonLabel.GetChildren() )
1695 {
1696 if( oObj.GetType() == CPLJSONObject::Type::Object )
1697 {
1698 CPLString osContainerName = oObj.GetName();
1699 CPLJSONObject oContainerName = oObj.GetObj( "_container_name" );
1700 if( oContainerName.GetType() == CPLJSONObject::Type::String )
1701 {
1702 osContainerName = oContainerName.ToString();
1703 }
1704
1705 CPLJSONObject oFilename = oObj.GetObj( "^" + osContainerName );
1706 if( oFilename.GetType() == CPLJSONObject::Type::String )
1707 {
1708 VSIStatBufL sStat;
1709 CPLString osFilename( CPLFormFilename(
1710 CPLGetPath(poOpenInfo->pszFilename),
1711 oFilename.ToString().c_str(),
1712 nullptr ) );
1713 if( VSIStatL( osFilename, &sStat ) == 0 )
1714 {
1715 poDS->m_aosAdditionalFiles.AddString(osFilename);
1716 }
1717 else
1718 {
1719 CPLDebug("ISIS3", "File %s referenced but not foud",
1720 osFilename.c_str());
1721 }
1722 }
1723 }
1724 }
1725
1726 VSIFCloseL( poOpenInfo->fpL );
1727 poOpenInfo->fpL = nullptr;
1728
1729 /* -------------------------------------------------------------------- */
1730 /* Assume user is pointing to label (i.e. .lbl) file for detached option */
1731 /* -------------------------------------------------------------------- */
1732 // Image can be inline or detached and point to an image name
1733 // the Format can be Tiled or Raw
1734 // Object = Core
1735 // StartByte = 65537
1736 // Format = Tile
1737 // TileSamples = 128
1738 // TileLines = 128
1739 //OR-----
1740 // Object = Core
1741 // StartByte = 1
1742 // ^Core = r0200357_detatched.cub
1743 // Format = BandSequential
1744 //OR-----
1745 // Object = Core
1746 // StartByte = 1
1747 // ^Core = r0200357_detached_tiled.cub
1748 // Format = Tile
1749 // TileSamples = 128
1750 // TileLines = 128
1751 //OR-----
1752 // Object = Core
1753 // StartByte = 1
1754 // ^Core = some.tif
1755 // Format = GeoTIFF
1756
1757 /* -------------------------------------------------------------------- */
1758 /* What file contains the actual data? */
1759 /* -------------------------------------------------------------------- */
1760 const char *pszCore = poDS->GetKeyword( "IsisCube.Core.^Core" );
1761 CPLString osQubeFile;
1762
1763 if( EQUAL(pszCore,"") )
1764 osQubeFile = poOpenInfo->pszFilename;
1765 else
1766 {
1767 CPLString osPath = CPLGetPath( poOpenInfo->pszFilename );
1768 osQubeFile = CPLFormFilename( osPath, pszCore, nullptr );
1769 poDS->m_osExternalFilename = osQubeFile;
1770 }
1771
1772 /* -------------------------------------------------------------------- */
1773 /* Check if file an ISIS3 header file? Read a few lines of text */
1774 /* searching for something starting with nrows or ncols. */
1775 /* -------------------------------------------------------------------- */
1776
1777 /************* Skipbytes *****************************/
1778 int nSkipBytes = atoi(poDS->GetKeyword("IsisCube.Core.StartByte", "1"));
1779 if( nSkipBytes <= 1 )
1780 nSkipBytes = 0;
1781 else
1782 nSkipBytes -= 1;
1783
1784 /******* Grab format type (BandSequential, Tiled) *******/
1785 CPLString osFormat = poDS->GetKeyword( "IsisCube.Core.Format" );
1786
1787 int tileSizeX = 0;
1788 int tileSizeY = 0;
1789
1790 if (EQUAL(osFormat,"Tile") )
1791 {
1792 poDS->m_bIsTiled = true;
1793 /******* Get Tile Sizes *********/
1794 tileSizeX = atoi(poDS->GetKeyword("IsisCube.Core.TileSamples"));
1795 tileSizeY = atoi(poDS->GetKeyword("IsisCube.Core.TileLines"));
1796 if (tileSizeX <= 0 || tileSizeY <= 0)
1797 {
1798 CPLError( CE_Failure, CPLE_OpenFailed,
1799 "Wrong tile dimensions : %d x %d",
1800 tileSizeX, tileSizeY);
1801 delete poDS;
1802 return nullptr;
1803 }
1804 }
1805 else if (!EQUAL(osFormat,"BandSequential") &&
1806 !EQUAL(osFormat,"GeoTIFF") )
1807 {
1808 CPLError( CE_Failure, CPLE_OpenFailed,
1809 "%s format not supported.", osFormat.c_str());
1810 delete poDS;
1811 return nullptr;
1812 }
1813
1814 /*********** Grab samples lines band ************/
1815 const int nCols = atoi(poDS->GetKeyword("IsisCube.Core.Dimensions.Samples"));
1816 const int nRows = atoi(poDS->GetKeyword("IsisCube.Core.Dimensions.Lines"));
1817 const int nBands = atoi(poDS->GetKeyword("IsisCube.Core.Dimensions.Bands"));
1818
1819 /****** Grab format type - ISIS3 only supports 8,U16,S16,32 *****/
1820 GDALDataType eDataType = GDT_Byte;
1821 double dfNoData = 0.0;
1822
1823 const char *itype = poDS->GetKeyword( "IsisCube.Core.Pixels.Type" );
1824 if (EQUAL(itype,"UnsignedByte") ) {
1825 eDataType = GDT_Byte;
1826 dfNoData = NULL1;
1827 }
1828 else if (EQUAL(itype,"UnsignedWord") ) {
1829 eDataType = GDT_UInt16;
1830 dfNoData = NULLU2;
1831 }
1832 else if (EQUAL(itype,"SignedWord") ) {
1833 eDataType = GDT_Int16;
1834 dfNoData = NULL2;
1835 }
1836 else if (EQUAL(itype,"Real") || EQUAL(itype,"") ) {
1837 eDataType = GDT_Float32;
1838 dfNoData = NULL4;
1839 }
1840 else {
1841 CPLError( CE_Failure, CPLE_OpenFailed,
1842 "%s pixel type not supported.", itype);
1843 delete poDS;
1844 return nullptr;
1845 }
1846
1847 /*********** Grab samples lines band ************/
1848
1849 //default to MSB
1850 const bool bIsLSB = EQUAL(
1851 poDS->GetKeyword( "IsisCube.Core.Pixels.ByteOrder"),"Lsb");
1852
1853 /*********** Grab Cellsize ************/
1854 double dfXDim = 1.0;
1855 double dfYDim = 1.0;
1856
1857 const char* pszRes = poDS->GetKeyword("IsisCube.Mapping.PixelResolution");
1858 if (strlen(pszRes) > 0 ) {
1859 dfXDim = CPLAtof(pszRes); /* values are in meters */
1860 dfYDim = -CPLAtof(pszRes);
1861 }
1862
1863 /*********** Grab UpperLeftCornerY ************/
1864 double dfULYMap = 0.5;
1865
1866 const char* pszULY = poDS->GetKeyword("IsisCube.Mapping.UpperLeftCornerY");
1867 if (strlen(pszULY) > 0) {
1868 dfULYMap = CPLAtof(pszULY);
1869 }
1870
1871 /*********** Grab UpperLeftCornerX ************/
1872 double dfULXMap = 0.5;
1873
1874 const char* pszULX = poDS->GetKeyword("IsisCube.Mapping.UpperLeftCornerX");
1875 if( strlen(pszULX) > 0 ) {
1876 dfULXMap = CPLAtof(pszULX);
1877 }
1878
1879 /*********** Grab TARGET_NAME ************/
1880 /**** This is the planets name i.e. Mars ***/
1881 const char *target_name = poDS->GetKeyword("IsisCube.Mapping.TargetName");
1882
1883 #ifdef notdef
1884 const double dfLongitudeMulFactor =
1885 EQUAL(poDS->GetKeyword( "IsisCube.Mapping.LongitudeDirection", "PositiveEast"), "PositiveEast") ? 1 : -1;
1886 #else
1887 const double dfLongitudeMulFactor = 1;
1888 #endif
1889
1890 /*********** Grab MAP_PROJECTION_TYPE ************/
1891 const char *map_proj_name =
1892 poDS->GetKeyword( "IsisCube.Mapping.ProjectionName");
1893
1894 /*********** Grab SEMI-MAJOR ************/
1895 const double semi_major =
1896 CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.EquatorialRadius"));
1897
1898 /*********** Grab semi-minor ************/
1899 const double semi_minor =
1900 CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.PolarRadius"));
1901
1902 /*********** Grab CENTER_LAT ************/
1903 const double center_lat =
1904 CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.CenterLatitude"));
1905
1906 /*********** Grab CENTER_LON ************/
1907 const double center_lon =
1908 CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.CenterLongitude")) * dfLongitudeMulFactor;
1909
1910 /*********** Grab 1st std parallel ************/
1911 const double first_std_parallel =
1912 CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.FirstStandardParallel"));
1913
1914 /*********** Grab 2nd std parallel ************/
1915 const double second_std_parallel =
1916 CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.SecondStandardParallel"));
1917
1918 /*********** Grab scaleFactor ************/
1919 const double scaleFactor =
1920 CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.scaleFactor", "1.0"));
1921
1922 /*** grab LatitudeType = Planetographic ****/
1923 // Need to further study how ocentric/ographic will effect the gdal library
1924 // So far we will use this fact to define a sphere or ellipse for some
1925 // projections
1926
1927 // Frank - may need to talk this over
1928 bool bIsGeographic = true;
1929 if (EQUAL( poDS->GetKeyword("IsisCube.Mapping.LatitudeType"),
1930 "Planetocentric" ))
1931 bIsGeographic = false;
1932
1933 //Set oSRS projection and parameters
1934 //############################################################
1935 //ISIS3 Projection types
1936 // Equirectangular
1937 // LambertConformal
1938 // Mercator
1939 // ObliqueCylindrical
1940 // Orthographic
1941 // PolarStereographic
1942 // SimpleCylindrical
1943 // Sinusoidal
1944 // TransverseMercator
1945
1946 #ifdef DEBUG
1947 CPLDebug( "ISIS3", "using projection %s", map_proj_name);
1948 #endif
1949
1950 OGRSpatialReference oSRS;
1951 bool bProjectionSet = true;
1952
1953 if ((EQUAL( map_proj_name, "Equirectangular" )) ||
1954 (EQUAL( map_proj_name, "SimpleCylindrical" )) ) {
1955 oSRS.SetEquirectangular2 ( 0.0, center_lon, center_lat, 0, 0 );
1956 } else if (EQUAL( map_proj_name, "Orthographic" )) {
1957 oSRS.SetOrthographic ( center_lat, center_lon, 0, 0 );
1958 } else if (EQUAL( map_proj_name, "Sinusoidal" )) {
1959 oSRS.SetSinusoidal ( center_lon, 0, 0 );
1960 } else if (EQUAL( map_proj_name, "Mercator" )) {
1961 oSRS.SetMercator ( center_lat, center_lon, scaleFactor, 0, 0 );
1962 } else if (EQUAL( map_proj_name, "PolarStereographic" )) {
1963 oSRS.SetPS ( center_lat, center_lon, scaleFactor, 0, 0 );
1964 } else if (EQUAL( map_proj_name, "TransverseMercator" )) {
1965 oSRS.SetTM ( center_lat, center_lon, scaleFactor, 0, 0 );
1966 } else if (EQUAL( map_proj_name, "LambertConformal" )) {
1967 oSRS.SetLCC ( first_std_parallel, second_std_parallel, center_lat, center_lon, 0, 0 );
1968 } else if (EQUAL( map_proj_name, "PointPerspective" )) {
1969 // Distance parameter is the distance to the center of the body, and is given in km
1970 const double distance = CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.Distance")) * 1000.0;
1971 const double height_above_ground = distance - semi_major;
1972 oSRS.SetVerticalPerspective(center_lat, center_lon, 0, height_above_ground, 0, 0);
1973 } else if (EQUAL( map_proj_name, "ObliqueCylindrical" )) {
1974 const double poleLatitude = CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.PoleLatitude"));
1975 const double poleLongitude = CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.PoleLongitude")) * dfLongitudeMulFactor;
1976 const double poleRotation = CPLAtof(poDS->GetKeyword( "IsisCube.Mapping.PoleRotation"));
1977 CPLString oProj4String;
1978 // ISIS3 rotated pole doesn't use the same conventions than PROJ ob_tran
1979 // Compare the sign difference in https://github.com/USGS-Astrogeology/ISIS3/blob/3.8.0/isis/src/base/objs/ObliqueCylindrical/ObliqueCylindrical.cpp#L244
1980 // and https://github.com/OSGeo/PROJ/blob/6.2/src/projections/ob_tran.cpp#L34
1981 // They can be compensated by modifying the poleLatitude to 180-poleLatitude
1982 // There's also a sign difference for the poleRotation parameter
1983 // The existence of those different conventions is acknowledged in
1984 // https://pds-imaging.jpl.nasa.gov/documentation/Cassini_BIDRSIS.PDF in the middle of page 10
1985 oProj4String.Printf(
1986 "+proj=ob_tran +o_proj=eqc +o_lon_p=%.18g +o_lat_p=%.18g +lon_0=%.18g",
1987 -poleRotation,
1988 180-poleLatitude,
1989 poleLongitude);
1990 oSRS.SetFromUserInput(oProj4String);
1991 } else {
1992 CPLDebug( "ISIS3",
1993 "Dataset projection %s is not supported. Continuing...",
1994 map_proj_name );
1995 bProjectionSet = false;
1996 }
1997
1998 if (bProjectionSet) {
1999 //Create projection name, i.e. MERCATOR MARS and set as ProjCS keyword
2000 CPLString osProjTargetName(map_proj_name);
2001 osProjTargetName += " ";
2002 osProjTargetName += target_name;
2003 oSRS.SetProjCS(osProjTargetName); //set ProjCS keyword
2004
2005 //The geographic/geocentric name will be the same basic name as the body name
2006 //'GCS' = Geographic/Geocentric Coordinate System
2007 CPLString osGeogName("GCS_");
2008 osGeogName += target_name;
2009
2010 //The datum name will be the same basic name as the planet
2011 CPLString osDatumName("D_");
2012 osDatumName += target_name;
2013
2014 CPLString osSphereName(target_name);
2015 //strcat(osSphereName, "_IAU_IAG"); //Might not be IAU defined so don't add
2016
2017 //calculate inverse flattening from major and minor axis: 1/f = a/(a-b)
2018 double iflattening = 0.0;
2019 if ((semi_major - semi_minor) < 0.0000001)
2020 iflattening = 0;
2021 else
2022 iflattening = semi_major / (semi_major - semi_minor);
2023
2024 //Set the body size but take into consideration which proj is being used to help w/ proj4 compatibility
2025 //The use of a Sphere, polar radius or ellipse here is based on how ISIS does it internally
2026 if ( ( (EQUAL( map_proj_name, "Stereographic" ) && (fabs(center_lat) == 90)) ) ||
2027 (EQUAL( map_proj_name, "PolarStereographic" )) )
2028 {
2029 if (bIsGeographic) {
2030 //Geograpraphic, so set an ellipse
2031 oSRS.SetGeogCS( osGeogName, osDatumName, osSphereName,
2032 semi_major, iflattening,
2033 "Reference_Meridian", 0.0 );
2034 } else {
2035 //Geocentric, so force a sphere using the semi-minor axis. I hope...
2036 osSphereName += "_polarRadius";
2037 oSRS.SetGeogCS( osGeogName, osDatumName, osSphereName,
2038 semi_minor, 0.0,
2039 "Reference_Meridian", 0.0 );
2040 }
2041 }
2042 else if ( (EQUAL( map_proj_name, "SimpleCylindrical" )) ||
2043 (EQUAL( map_proj_name, "Orthographic" )) ||
2044 (EQUAL( map_proj_name, "Stereographic" )) ||
2045 (EQUAL( map_proj_name, "Sinusoidal" )) ||
2046 (EQUAL( map_proj_name, "PointPerspective" )) ) {
2047 // ISIS uses the spherical equation for these projections
2048 // so force a sphere.
2049 oSRS.SetGeogCS( osGeogName, osDatumName, osSphereName,
2050 semi_major, 0.0,
2051 "Reference_Meridian", 0.0 );
2052 }
2053 else if (EQUAL( map_proj_name, "Equirectangular" )) {
2054 //Calculate localRadius using ISIS3 simple elliptical method
2055 // not the more standard Radius of Curvature method
2056 //PI = 4 * atan(1);
2057 const double radLat = center_lat * M_PI / 180; // in radians
2058 const double meanRadius =
2059 sqrt( pow( semi_minor * cos( radLat ), 2)
2060 + pow( semi_major * sin( radLat ), 2) );
2061 const double localRadius = ( meanRadius == 0.0 ) ?
2062 0.0 : semi_major * semi_minor / meanRadius;
2063 osSphereName += "_localRadius";
2064 oSRS.SetGeogCS( osGeogName, osDatumName, osSphereName,
2065 localRadius, 0.0,
2066 "Reference_Meridian", 0.0 );
2067 }
2068 else {
2069 //All other projections: Mercator, Transverse Mercator, Lambert Conformal, etc.
2070 //Geographic, so set an ellipse
2071 if (bIsGeographic) {
2072 oSRS.SetGeogCS( osGeogName, osDatumName, osSphereName,
2073 semi_major, iflattening,
2074 "Reference_Meridian", 0.0 );
2075 } else {
2076 //Geocentric, so force a sphere. I hope...
2077 oSRS.SetGeogCS( osGeogName, osDatumName, osSphereName,
2078 semi_major, 0.0,
2079 "Reference_Meridian", 0.0 );
2080 }
2081 }
2082
2083 // translate back into a projection string.
2084 poDS->m_oSRS = oSRS;
2085 poDS->m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2086 }
2087
2088 /* END ISIS3 Label Read */
2089 /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
2090
2091 /* -------------------------------------------------------------------- */
2092 /* Did we get the required keywords? If not we return with */
2093 /* this never having been considered to be a match. This isn't */
2094 /* an error! */
2095 /* -------------------------------------------------------------------- */
2096 if( !GDALCheckDatasetDimensions(nCols, nRows) ||
2097 !GDALCheckBandCount(nBands, false) )
2098 {
2099 delete poDS;
2100 return nullptr;
2101 }
2102
2103 /* -------------------------------------------------------------------- */
2104 /* Capture some information from the file that is of interest. */
2105 /* -------------------------------------------------------------------- */
2106 poDS->nRasterXSize = nCols;
2107 poDS->nRasterYSize = nRows;
2108
2109 /* -------------------------------------------------------------------- */
2110 /* Open target binary file. */
2111 /* -------------------------------------------------------------------- */
2112 if( EQUAL(osFormat,"GeoTIFF") )
2113 {
2114 if( nSkipBytes != 0 )
2115 {
2116 CPLError(CE_Warning, CPLE_NotSupported,
2117 "Ignoring StartByte=%d for format=GeoTIFF",
2118 1+nSkipBytes);
2119 }
2120 if( osQubeFile == poOpenInfo->pszFilename )
2121 {
2122 CPLError( CE_Failure, CPLE_AppDefined,
2123 "A ^Core file must be set");
2124 delete poDS;
2125 return nullptr;
2126 }
2127 poDS->m_poExternalDS = reinterpret_cast<GDALDataset *>(
2128 GDALOpen( osQubeFile, poOpenInfo->eAccess ) );
2129 if( poDS->m_poExternalDS == nullptr )
2130 {
2131 delete poDS;
2132 return nullptr;
2133 }
2134 if( poDS->m_poExternalDS->GetRasterXSize() != poDS->nRasterXSize ||
2135 poDS->m_poExternalDS->GetRasterYSize() != poDS->nRasterYSize ||
2136 poDS->m_poExternalDS->GetRasterCount() != nBands ||
2137 poDS->m_poExternalDS->GetRasterBand(1)->GetRasterDataType() !=
2138 eDataType )
2139 {
2140 CPLError( CE_Failure, CPLE_AppDefined,
2141 "%s has incompatible characteristics with the ones "
2142 "declared in the label.",
2143 osQubeFile.c_str() );
2144 delete poDS;
2145 return nullptr;
2146 }
2147 }
2148 else
2149 {
2150 if( poOpenInfo->eAccess == GA_ReadOnly )
2151 poDS->m_fpImage = VSIFOpenL( osQubeFile, "r" );
2152 else
2153 poDS->m_fpImage = VSIFOpenL( osQubeFile, "r+" );
2154
2155 if( poDS->m_fpImage == nullptr )
2156 {
2157 CPLError( CE_Failure, CPLE_OpenFailed,
2158 "Failed to open %s: %s.",
2159 osQubeFile.c_str(),
2160 VSIStrerror( errno ) );
2161 delete poDS;
2162 return nullptr;
2163 }
2164
2165 // Sanity checks in case the external raw file appears to be a
2166 // TIFF file
2167 if( EQUAL(CPLGetExtension(osQubeFile), "tif") )
2168 {
2169 GDALDataset* poTIF_DS = reinterpret_cast<GDALDataset*>(
2170 GDALOpen(osQubeFile, GA_ReadOnly));
2171 if( poTIF_DS )
2172 {
2173 bool bWarned = false;
2174 if( poTIF_DS->GetRasterXSize() != poDS->nRasterXSize ||
2175 poTIF_DS->GetRasterYSize() != poDS->nRasterYSize ||
2176 poTIF_DS->GetRasterCount() != nBands ||
2177 poTIF_DS->GetRasterBand(1)->GetRasterDataType() !=
2178 eDataType ||
2179 poTIF_DS->GetMetadataItem("COMPRESSION",
2180 "IMAGE_STRUCTURE") != nullptr )
2181 {
2182 bWarned = true;
2183 CPLError( CE_Warning, CPLE_AppDefined,
2184 "%s has incompatible characteristics with the ones "
2185 "declared in the label.",
2186 osQubeFile.c_str() );
2187 }
2188 int nBlockXSize = 1, nBlockYSize = 1;
2189 poTIF_DS->GetRasterBand(1)->GetBlockSize(&nBlockXSize,
2190 &nBlockYSize);
2191 if( (poDS->m_bIsTiled && (nBlockXSize != tileSizeX ||
2192 nBlockYSize != tileSizeY) ) ||
2193 (!poDS->m_bIsTiled && (nBlockXSize != nCols ||
2194 (nBands > 1 && nBlockYSize != 1))) )
2195 {
2196 if( !bWarned )
2197 {
2198 bWarned = true;
2199 CPLError( CE_Warning, CPLE_AppDefined,
2200 "%s has incompatible characteristics with the ones "
2201 "declared in the label.",
2202 osQubeFile.c_str() );
2203 }
2204 }
2205 // to please Clang Static Analyzer
2206 nBlockXSize = std::max(1, nBlockXSize);
2207 nBlockYSize = std::max(1, nBlockYSize);
2208
2209 // Check that blocks are effectively written in expected order.
2210 const int nBlockSizeBytes = nBlockXSize * nBlockYSize *
2211 GDALGetDataTypeSizeBytes(eDataType);
2212 bool bGoOn = !bWarned;
2213 const int l_nBlocksPerRow =
2214 DIV_ROUND_UP(nCols, nBlockXSize);
2215 const int l_nBlocksPerColumn =
2216 DIV_ROUND_UP(nRows, nBlockYSize);
2217 int nBlockNo = 0;
2218 for( int i = 0; i < nBands && bGoOn; i++ )
2219 {
2220 for( int y = 0; y < l_nBlocksPerColumn && bGoOn; y++ )
2221 {
2222 for( int x = 0; x < l_nBlocksPerRow && bGoOn; x++ )
2223 {
2224 const char* pszBlockOffset = poTIF_DS->
2225 GetRasterBand(i+1)->GetMetadataItem(
2226 CPLSPrintf("BLOCK_OFFSET_%d_%d", x, y),
2227 "TIFF");
2228 if( pszBlockOffset )
2229 {
2230 GIntBig nOffset = CPLAtoGIntBig(pszBlockOffset);
2231 if( nOffset != nSkipBytes + nBlockNo *
2232 nBlockSizeBytes )
2233 {
2234 //bWarned = true;
2235 CPLError( CE_Warning, CPLE_AppDefined,
2236 "%s has incompatible "
2237 "characteristics with the ones "
2238 "declared in the label.",
2239 osQubeFile.c_str() );
2240 bGoOn = false;
2241 }
2242 }
2243 nBlockNo ++;
2244 }
2245 }
2246 }
2247
2248 delete poTIF_DS;
2249 }
2250 }
2251 }
2252
2253 poDS->eAccess = poOpenInfo->eAccess;
2254
2255 /* -------------------------------------------------------------------- */
2256 /* Compute the line offset. */
2257 /* -------------------------------------------------------------------- */
2258 int nLineOffset = 0;
2259 int nPixelOffset = 0;
2260 vsi_l_offset nBandOffset = 0;
2261
2262 if( EQUAL(osFormat,"BandSequential") )
2263 {
2264 const int nItemSize = GDALGetDataTypeSizeBytes(eDataType);
2265 nPixelOffset = nItemSize;
2266 try
2267 {
2268 nLineOffset = (CPLSM(nPixelOffset) * CPLSM(nCols)).v();
2269 }
2270 catch( const CPLSafeIntOverflow& )
2271 {
2272 delete poDS;
2273 return nullptr;
2274 }
2275 nBandOffset = static_cast<vsi_l_offset>(nLineOffset) * nRows;
2276
2277 poDS->m_sLayout.osRawFilename = osQubeFile;
2278 if( nBands > 1 )
2279 poDS->m_sLayout.eInterleaving = RawBinaryLayout::Interleaving::BSQ;
2280 poDS->m_sLayout.eDataType = eDataType;
2281 poDS->m_sLayout.bLittleEndianOrder = bIsLSB;
2282 poDS->m_sLayout.nImageOffset = nSkipBytes;
2283 poDS->m_sLayout.nPixelOffset = nPixelOffset;
2284 poDS->m_sLayout.nLineOffset = nLineOffset;
2285 poDS->m_sLayout.nBandOffset = static_cast<GIntBig>(nBandOffset);
2286 }
2287 /* else Tiled or external */
2288
2289 /* -------------------------------------------------------------------- */
2290 /* Extract BandBin info. */
2291 /* -------------------------------------------------------------------- */
2292 std::vector<std::string> aosBandNames;
2293 std::vector<std::string> aosBandUnits;
2294 std::vector<double> adfWavelengths;
2295 std::vector<std::string> aosWavelengthsUnit;
2296 std::vector<double> adfBandwidth;
2297 std::vector<std::string> aosBandwidthUnit;
2298 const auto oBandBin = poDS->m_oJSonLabel.GetObj( "IsisCube/BandBin" );
2299 if( oBandBin.IsValid() && oBandBin.GetType() == CPLJSONObject::Type::Object )
2300 {
2301 for( const auto& child: oBandBin.GetChildren() )
2302 {
2303 if( CPLString(child.GetName()).ifind("name") != std::string::npos )
2304 {
2305 // Use "name" in priority
2306 if( EQUAL(child.GetName().c_str(), "name") )
2307 {
2308 aosBandNames.clear();
2309 }
2310 else if( !aosBandNames.empty() )
2311 {
2312 continue;
2313 }
2314
2315 if( child.GetType() == CPLJSONObject::Type::String && nBands == 1 )
2316 {
2317 aosBandNames.push_back(child.ToString());
2318 }
2319 else if( child.GetType() == CPLJSONObject::Type::Array )
2320 {
2321 auto oArray = child.ToArray();
2322 if( oArray.Size() == nBands )
2323 {
2324 for( int i = 0; i < nBands; i++ )
2325 {
2326 if( oArray[i].GetType() == CPLJSONObject::Type::String )
2327 {
2328 aosBandNames.push_back(oArray[i].ToString());
2329 }
2330 else
2331 {
2332 aosBandNames.clear();
2333 break;
2334 }
2335 }
2336 }
2337 }
2338 }
2339 else if( EQUAL(child.GetName().c_str(), "BandSuffixUnit") &&
2340 child.GetType() == CPLJSONObject::Type::Array )
2341 {
2342 auto oArray = child.ToArray();
2343 if( oArray.Size() == nBands )
2344 {
2345 for( int i = 0; i < nBands; i++ )
2346 {
2347 if( oArray[i].GetType() == CPLJSONObject::Type::String )
2348 {
2349 aosBandUnits.push_back(oArray[i].ToString());
2350 }
2351 else
2352 {
2353 aosBandUnits.clear();
2354 break;
2355 }
2356 }
2357 }
2358 }
2359 else if( EQUAL(child.GetName().c_str(), "BandBinCenter") ||
2360 EQUAL(child.GetName().c_str(), "Center") )
2361 {
2362 GetValueAndUnits(child, adfWavelengths, aosWavelengthsUnit,
2363 nBands);
2364 }
2365 else if( EQUAL(child.GetName().c_str(), "BandBinUnit") &&
2366 child.GetType() == CPLJSONObject::Type::String )
2367 {
2368 CPLString unit(child.ToString());
2369 if( STARTS_WITH_CI(unit, "micromet") ||
2370 EQUAL(unit, "um") ||
2371 STARTS_WITH_CI(unit, "nanomet") ||
2372 EQUAL(unit, "nm") )
2373 {
2374 aosWavelengthsUnit.push_back(child.ToString());
2375 }
2376 }
2377 else if( EQUAL(child.GetName().c_str(), "Width") )
2378 {
2379 GetValueAndUnits(child, adfBandwidth, aosBandwidthUnit,
2380 nBands);
2381 }
2382 }
2383
2384 if( !adfWavelengths.empty() && aosWavelengthsUnit.size() == 1 )
2385 {
2386 for( int i = 1; i < nBands; i++ )
2387 {
2388 aosWavelengthsUnit.push_back(aosWavelengthsUnit[0]);
2389 }
2390 }
2391 if( !adfBandwidth.empty() && aosBandwidthUnit.size() == 1 )
2392 {
2393 for( int i = 1; i < nBands; i++ )
2394 {
2395 aosBandwidthUnit.push_back(aosBandwidthUnit[0]);
2396 }
2397 }
2398 }
2399
2400 /* -------------------------------------------------------------------- */
2401 /* Create band information objects. */
2402 /* -------------------------------------------------------------------- */
2403 #ifdef CPL_LSB
2404 const bool bNativeOrder = bIsLSB;
2405 #else
2406 const bool bNativeOrder = !bIsLSB;
2407 #endif
2408
2409 for( int i = 0; i < nBands; i++ )
2410 {
2411 GDALRasterBand *poBand = nullptr;
2412
2413 if( poDS->m_poExternalDS != nullptr )
2414 {
2415 ISIS3WrapperRasterBand* poISISBand =
2416 new ISIS3WrapperRasterBand(
2417 poDS->m_poExternalDS->GetRasterBand( i+1 ) );
2418 poBand = poISISBand;
2419 poDS->SetBand( i+1, poBand );
2420
2421 poISISBand->SetMaskBand( new ISISMaskBand(poISISBand) );
2422 }
2423 else if( poDS->m_bIsTiled )
2424 {
2425 CPLErrorReset();
2426 ISISTiledBand* poISISBand =
2427 new ISISTiledBand( poDS, poDS->m_fpImage, i+1, eDataType,
2428 tileSizeX, tileSizeY,
2429 nSkipBytes, 0, 0,
2430 bNativeOrder );
2431 if( CPLGetLastErrorType() != CE_None )
2432 {
2433 delete poISISBand;
2434 delete poDS;
2435 return nullptr;
2436 }
2437 poBand = poISISBand;
2438 poDS->SetBand( i+1, poBand );
2439
2440 poISISBand->SetMaskBand( new ISISMaskBand(poISISBand) );
2441 }
2442 else
2443 {
2444 ISIS3RawRasterBand* poISISBand =
2445 new ISIS3RawRasterBand( poDS, i+1, poDS->m_fpImage,
2446 nSkipBytes + nBandOffset * i,
2447 nPixelOffset, nLineOffset, eDataType,
2448 bNativeOrder );
2449
2450 poBand = poISISBand;
2451 poDS->SetBand( i+1, poBand );
2452
2453 poISISBand->SetMaskBand( new ISISMaskBand(poISISBand) );
2454 }
2455
2456 if( i < static_cast<int>(aosBandNames.size()) )
2457 {
2458 poBand->SetDescription(aosBandNames[i].c_str());
2459 }
2460 if( i < static_cast<int>(adfWavelengths.size()) &&
2461 i < static_cast<int>(aosWavelengthsUnit.size()) )
2462 {
2463 poBand->SetMetadataItem("WAVELENGTH", CPLSPrintf("%f", adfWavelengths[i]));
2464 poBand->SetMetadataItem("WAVELENGTH_UNIT", aosWavelengthsUnit[i].c_str());
2465 if( i < static_cast<int>(adfBandwidth.size()) &&
2466 i < static_cast<int>(aosBandwidthUnit.size()) )
2467 {
2468 poBand->SetMetadataItem("BANDWIDTH", CPLSPrintf("%f", adfBandwidth[i]));
2469 poBand->SetMetadataItem("BANDWIDTH_UNIT", aosBandwidthUnit[i].c_str());
2470 }
2471 }
2472 if( i < static_cast<int>(aosBandUnits.size()) )
2473 {
2474 poBand->SetUnitType(aosBandUnits[i].c_str());
2475 }
2476
2477 poBand->SetNoDataValue( dfNoData );
2478
2479 // Set offset/scale values.
2480 const double dfOffset =
2481 CPLAtofM(poDS->GetKeyword("IsisCube.Core.Pixels.Base","0.0"));
2482 const double dfScale =
2483 CPLAtofM(poDS->GetKeyword("IsisCube.Core.Pixels.Multiplier","1.0"));
2484 if( dfOffset != 0.0 || dfScale != 1.0 )
2485 {
2486 poBand->SetOffset(dfOffset);
2487 poBand->SetScale(dfScale);
2488 }
2489 }
2490
2491 /* -------------------------------------------------------------------- */
2492 /* Check for a .prj file. For ISIS3 I would like to keep this in */
2493 /* -------------------------------------------------------------------- */
2494 const CPLString osPath = CPLGetPath( poOpenInfo->pszFilename );
2495 const CPLString osName = CPLGetBasename(poOpenInfo->pszFilename);
2496 const char *pszPrjFile = CPLFormCIFilename( osPath, osName, "prj" );
2497
2498 VSILFILE *fp = VSIFOpenL( pszPrjFile, "r" );
2499 if( fp != nullptr )
2500 {
2501 VSIFCloseL( fp );
2502
2503 char **papszLines = CSLLoad( pszPrjFile );
2504
2505 OGRSpatialReference oSRS2;
2506 if( oSRS2.importFromESRI( papszLines ) == OGRERR_NONE )
2507 {
2508 poDS->m_aosAdditionalFiles.AddString( pszPrjFile );
2509 poDS->m_oSRS = oSRS2;
2510 poDS->m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2511 }
2512
2513 CSLDestroy( papszLines );
2514 }
2515
2516 if( dfULXMap != 0.5 || dfULYMap != 0.5 || dfXDim != 1.0 || dfYDim != 1.0 )
2517 {
2518 poDS->m_bGotTransform = true;
2519 poDS->m_adfGeoTransform[0] = dfULXMap;
2520 poDS->m_adfGeoTransform[1] = dfXDim;
2521 poDS->m_adfGeoTransform[2] = 0.0;
2522 poDS->m_adfGeoTransform[3] = dfULYMap;
2523 poDS->m_adfGeoTransform[4] = 0.0;
2524 poDS->m_adfGeoTransform[5] = dfYDim;
2525 }
2526
2527 if( !poDS->m_bGotTransform )
2528 {
2529 poDS->m_bGotTransform =
2530 CPL_TO_BOOL(GDALReadWorldFile( poOpenInfo->pszFilename, "cbw",
2531 poDS->m_adfGeoTransform ));
2532 if( poDS->m_bGotTransform )
2533 {
2534 poDS->m_aosAdditionalFiles.AddString(
2535 CPLResetExtension(poOpenInfo->pszFilename, "cbw") );
2536 }
2537 }
2538
2539 if( !poDS->m_bGotTransform )
2540 {
2541 poDS->m_bGotTransform =
2542 CPL_TO_BOOL(GDALReadWorldFile( poOpenInfo->pszFilename, "wld",
2543 poDS->m_adfGeoTransform ));
2544 if( poDS->m_bGotTransform )
2545 {
2546 poDS->m_aosAdditionalFiles.AddString(
2547 CPLResetExtension(poOpenInfo->pszFilename, "wld") );
2548 }
2549 }
2550
2551 /* -------------------------------------------------------------------- */
2552 /* Initialize any PAM information. */
2553 /* -------------------------------------------------------------------- */
2554 poDS->SetDescription( poOpenInfo->pszFilename );
2555 poDS->TryLoadXML();
2556
2557 /* -------------------------------------------------------------------- */
2558 /* Check for overviews. */
2559 /* -------------------------------------------------------------------- */
2560 poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
2561
2562 return poDS;
2563 }
2564
2565 /************************************************************************/
2566 /* GetKeyword() */
2567 /************************************************************************/
2568
GetKeyword(const char * pszPath,const char * pszDefault)2569 const char *ISIS3Dataset::GetKeyword( const char *pszPath,
2570 const char *pszDefault )
2571
2572 {
2573 return m_oKeywords.GetKeyword( pszPath, pszDefault );
2574 }
2575
2576 /************************************************************************/
2577 /* FixLong() */
2578 /************************************************************************/
2579
FixLong(double dfLong)2580 double ISIS3Dataset::FixLong( double dfLong )
2581 {
2582 if( m_osLongitudeDirection == "PositiveWest" )
2583 dfLong = -dfLong;
2584 if( m_bForce360 && dfLong < 0 )
2585 dfLong += 360.0;
2586 return dfLong;
2587 }
2588
2589 /************************************************************************/
2590 /* BuildLabel() */
2591 /************************************************************************/
2592
BuildLabel()2593 void ISIS3Dataset::BuildLabel()
2594 {
2595 CPLJSONObject oLabel = m_oSrcJSonLabel;
2596 if( !oLabel.IsValid() )
2597 {
2598 oLabel = CPLJSONObject();
2599 }
2600 // If we have a source label, then edit it directly
2601 CPLJSONObject oIsisCube = GetOrCreateJSONObject(oLabel, "IsisCube");
2602 oIsisCube.Set( "_type", "object");
2603
2604 if( !m_osComment.empty() )
2605 oIsisCube.Set( "_comment", m_osComment );
2606
2607 CPLJSONObject oCore = GetOrCreateJSONObject(oIsisCube, "Core");
2608 if( oCore.GetType() != CPLJSONObject::Type::Object )
2609 {
2610 oIsisCube.Delete( "Core" );
2611 oCore = CPLJSONObject();
2612 oIsisCube.Add("Core", oCore);
2613 }
2614 oCore.Set( "_type", "object" );
2615
2616 if( !m_osExternalFilename.empty() )
2617 {
2618 if( m_poExternalDS && m_bGeoTIFFAsRegularExternal )
2619 {
2620 if( !m_bGeoTIFFInitDone )
2621 {
2622 reinterpret_cast<ISIS3WrapperRasterBand*>(GetRasterBand(1))->
2623 InitFile();
2624 }
2625
2626 const char* pszOffset = m_poExternalDS->GetRasterBand(1)->
2627 GetMetadataItem("BLOCK_OFFSET_0_0", "TIFF");
2628 if( pszOffset )
2629 {
2630 oCore.Set( "StartByte", 1 + atoi(pszOffset) );
2631 }
2632 else
2633 {
2634 // Shouldn't happen normally
2635 CPLError(CE_Warning, CPLE_AppDefined,
2636 "Missing BLOCK_OFFSET_0_0");
2637 m_bGeoTIFFAsRegularExternal = false;
2638 oCore.Set( "StartByte", 1 );
2639 }
2640 }
2641 else
2642 {
2643 oCore.Set( "StartByte", 1 );
2644 }
2645 if( !m_osExternalFilename.empty() )
2646 {
2647 const CPLString osExternalFilename = CPLGetFilename(m_osExternalFilename);
2648 oCore.Set( "^Core", osExternalFilename );
2649 }
2650 }
2651 else
2652 {
2653 oCore.Set( "StartByte", pszSTARTBYTE_PLACEHOLDER );
2654 oCore.Delete( "^Core" );
2655 }
2656
2657 if( m_poExternalDS && !m_bGeoTIFFAsRegularExternal )
2658 {
2659 oCore.Set( "Format", "GeoTIFF" );
2660 oCore.Delete( "TileSamples" );
2661 oCore.Delete( "TileLines" );
2662 }
2663 else if( m_bIsTiled )
2664 {
2665 oCore.Set( "Format", "Tile");
2666 int nBlockXSize = 1, nBlockYSize = 1;
2667 GetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
2668 oCore.Set( "TileSamples", nBlockXSize );
2669 oCore.Set( "TileLines", nBlockYSize );
2670 }
2671 else
2672 {
2673 oCore.Set( "Format", "BandSequential" );
2674 oCore.Delete( "TileSamples" );
2675 oCore.Delete( "TileLines" );
2676 }
2677
2678 CPLJSONObject oDimensions = GetOrCreateJSONObject(oCore, "Dimensions");
2679 oDimensions.Set( "_type", "group" );
2680 oDimensions.Set( "Samples", nRasterXSize );
2681 oDimensions.Set( "Lines", nRasterYSize );
2682 oDimensions.Set( "Bands", nBands );
2683
2684 CPLJSONObject oPixels = GetOrCreateJSONObject(oCore, "Pixels");
2685 oPixels.Set( "_type", "group" );
2686 const GDALDataType eDT = GetRasterBand(1)->GetRasterDataType();
2687 oPixels.Set( "Type",
2688 (eDT == GDT_Byte) ? "UnsignedByte" :
2689 (eDT == GDT_UInt16) ? "UnsignedWord" :
2690 (eDT == GDT_Int16) ? "SignedWord" :
2691 "Real" );
2692
2693 oPixels.Set( "ByteOrder", "Lsb" );
2694 oPixels.Set( "Base", GetRasterBand(1)->GetOffset() );
2695 oPixels.Set( "Multiplier", GetRasterBand(1)->GetScale() );
2696
2697 const OGRSpatialReference& oSRS = m_oSRS;
2698
2699 if( !m_bUseSrcMapping )
2700 {
2701 oIsisCube.Delete( "Mapping" );
2702 }
2703
2704 CPLJSONObject oMapping = GetOrCreateJSONObject(oIsisCube, "Mapping");
2705 if( m_bUseSrcMapping && oMapping.IsValid() &&
2706 oMapping.GetType() == CPLJSONObject::Type::Object )
2707 {
2708 if( !m_osTargetName.empty() )
2709 oMapping.Set( "TargetName", m_osTargetName );
2710 if( !m_osLatitudeType.empty() )
2711 oMapping.Set( "LatitudeType", m_osLatitudeType );
2712 if( !m_osLongitudeDirection.empty() )
2713 oMapping.Set( "LongitudeDirection", m_osLongitudeDirection );
2714 }
2715 else if( !m_bUseSrcMapping && !m_oSRS.IsEmpty() )
2716 {
2717 oMapping.Add( "_type", "group" );
2718
2719 if( oSRS.IsProjected() || oSRS.IsGeographic() )
2720 {
2721 const char* pszDatum = oSRS.GetAttrValue("DATUM");
2722 CPLString osTargetName( m_osTargetName );
2723 if( osTargetName.empty() )
2724 {
2725 if( pszDatum && STARTS_WITH(pszDatum, "D_") )
2726 {
2727 osTargetName = pszDatum + 2;
2728 }
2729 else if( pszDatum )
2730 {
2731 osTargetName = pszDatum;
2732 }
2733 }
2734 if( !osTargetName.empty() )
2735 oMapping.Add( "TargetName", osTargetName );
2736
2737 oMapping.Add( "EquatorialRadius/value", oSRS.GetSemiMajor() );
2738 oMapping.Add( "EquatorialRadius/unit", "meters" );
2739 oMapping.Add( "PolarRadius/value", oSRS.GetSemiMinor() );
2740 oMapping.Add( "PolarRadius/unit", "meters" );
2741
2742 if( !m_osLatitudeType.empty() )
2743 oMapping.Add( "LatitudeType", m_osLatitudeType );
2744 else
2745 oMapping.Add( "LatitudeType", "Planetocentric" );
2746
2747 if( !m_osLongitudeDirection.empty() )
2748 oMapping.Add( "LongitudeDirection", m_osLongitudeDirection );
2749 else
2750 oMapping.Add( "LongitudeDirection", "PositiveEast" );
2751
2752 double adfX[4] = {0};
2753 double adfY[4] = {0};
2754 bool bLongLatCorners = false;
2755 if( m_bGotTransform )
2756 {
2757 for( int i = 0; i < 4; i++ )
2758 {
2759 adfX[i] = m_adfGeoTransform[0] + (i%2) *
2760 nRasterXSize * m_adfGeoTransform[1];
2761 adfY[i] = m_adfGeoTransform[3] +
2762 ( (i == 0 || i == 3) ? 0 : 1 ) *
2763 nRasterYSize * m_adfGeoTransform[5];
2764 }
2765 if( oSRS.IsGeographic() )
2766 {
2767 bLongLatCorners = true;
2768 }
2769 else
2770 {
2771 OGRSpatialReference* poSRSLongLat = oSRS.CloneGeogCS();
2772 if( poSRSLongLat )
2773 {
2774 poSRSLongLat->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2775 OGRCoordinateTransformation* poCT =
2776 OGRCreateCoordinateTransformation(&oSRS, poSRSLongLat);
2777 if( poCT )
2778 {
2779 if( poCT->Transform(4, adfX, adfY) )
2780 {
2781 bLongLatCorners = true;
2782 }
2783 delete poCT;
2784 }
2785 delete poSRSLongLat;
2786 }
2787 }
2788 }
2789 if( bLongLatCorners )
2790 {
2791 for( int i = 0; i < 4; i++ )
2792 {
2793 adfX[i] = FixLong(adfX[i]);
2794 }
2795 }
2796
2797 if( bLongLatCorners && (
2798 m_bForce360 || adfX[0] <- 180.0 || adfX[3] > 180.0) )
2799 {
2800 oMapping.Add( "LongitudeDomain", 360 );
2801 }
2802 else
2803 {
2804 oMapping.Add( "LongitudeDomain", 180 );
2805 }
2806
2807 if( m_bWriteBoundingDegrees && !m_osBoundingDegrees.empty() )
2808 {
2809 char** papszTokens =
2810 CSLTokenizeString2(m_osBoundingDegrees, ",", 0);
2811 if( CSLCount(papszTokens) == 4 )
2812 {
2813 oMapping.Add( "MinimumLatitude", CPLAtof(papszTokens[1]) );
2814 oMapping.Add( "MinimumLongitude", CPLAtof(papszTokens[0]) );
2815 oMapping.Add( "MaximumLatitude", CPLAtof(papszTokens[3]) );
2816 oMapping.Add( "MaximumLongitude", CPLAtof(papszTokens[2]) );
2817 }
2818 CSLDestroy(papszTokens);
2819 }
2820 else if( m_bWriteBoundingDegrees && bLongLatCorners )
2821 {
2822 oMapping.Add( "MinimumLatitude", std::min(
2823 std::min(adfY[0], adfY[1]), std::min(adfY[2],adfY[3])) );
2824 oMapping.Add( "MinimumLongitude", std::min(
2825 std::min(adfX[0], adfX[1]), std::min(adfX[2],adfX[3])) );
2826 oMapping.Add( "MaximumLatitude", std::max(
2827 std::max(adfY[0], adfY[1]), std::max(adfY[2],adfY[3])) );
2828 oMapping.Add( "MaximumLongitude", std::max(
2829 std::max(adfX[0], adfX[1]), std::max(adfX[2],adfX[3])) );
2830 }
2831
2832 const char* pszProjection = oSRS.GetAttrValue("PROJECTION");
2833 if( pszProjection == nullptr )
2834 {
2835 oMapping.Add( "ProjectionName", "SimpleCylindrical" );
2836 oMapping.Add( "CenterLongitude", 0.0 );
2837 oMapping.Add( "CenterLatitude", 0.0 );
2838 oMapping.Add( "CenterLatitudeRadius", oSRS.GetSemiMajor() );
2839 }
2840 else if( EQUAL(pszProjection, SRS_PT_EQUIRECTANGULAR) )
2841 {
2842 oMapping.Add( "ProjectionName", "Equirectangular" );
2843 if( oSRS.GetNormProjParm( SRS_PP_LATITUDE_OF_ORIGIN, 0.0 )
2844 != 0.0 )
2845 {
2846 CPLError(CE_Warning, CPLE_NotSupported,
2847 "Ignoring %s. Only 0 value supported",
2848 SRS_PP_LATITUDE_OF_ORIGIN);
2849 }
2850 oMapping.Add( "CenterLongitude",
2851 FixLong(oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)) );
2852 const double dfCenterLat =
2853 oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0);
2854 oMapping.Add( "CenterLatitude", dfCenterLat );
2855
2856 // in radians
2857 const double radLat = dfCenterLat * M_PI / 180;
2858 const double semi_major = oSRS.GetSemiMajor();
2859 const double semi_minor = oSRS.GetSemiMinor();
2860 const double localRadius
2861 = semi_major * semi_minor
2862 / sqrt( pow( semi_minor * cos( radLat ), 2)
2863 + pow( semi_major * sin( radLat ), 2) );
2864 oMapping.Add( "CenterLatitudeRadius", localRadius );
2865 }
2866
2867 else if( EQUAL(pszProjection, SRS_PT_ORTHOGRAPHIC) )
2868 {
2869 oMapping.Add( "ProjectionName", "Orthographic" );
2870 oMapping.Add( "CenterLongitude", FixLong(
2871 oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)) );
2872 oMapping.Add( "CenterLatitude",
2873 oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0) );
2874 }
2875
2876 else if( EQUAL(pszProjection, SRS_PT_SINUSOIDAL) )
2877 {
2878 oMapping.Add( "ProjectionName", "Sinusoidal" );
2879 oMapping.Add( "CenterLongitude", FixLong(
2880 oSRS.GetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER, 0.0)) );
2881 }
2882
2883 else if( EQUAL(pszProjection, SRS_PT_MERCATOR_1SP) )
2884 {
2885 oMapping.Add( "ProjectionName", "Mercator" );
2886 oMapping.Add( "CenterLongitude", FixLong(
2887 oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)) );
2888 oMapping.Add( "CenterLatitude",
2889 oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0) );
2890 oMapping.Add( "scaleFactor",
2891 oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0) );
2892 }
2893
2894 else if( EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC) )
2895 {
2896 oMapping.Add( "ProjectionName", "PolarStereographic" );
2897 oMapping.Add( "CenterLongitude", FixLong(
2898 oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)) );
2899 oMapping.Add( "CenterLatitude",
2900 oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0) );
2901 oMapping.Add( "scaleFactor",
2902 oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0) );
2903 }
2904
2905 else if( EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR) )
2906 {
2907 oMapping.Add( "ProjectionName", "TransverseMercator" );
2908 oMapping.Add( "CenterLongitude", FixLong(
2909 oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)) );
2910 oMapping.Add( "CenterLatitude",
2911 oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0) );
2912 oMapping.Add( "scaleFactor",
2913 oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0) );
2914 }
2915
2916 else if( EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) )
2917 {
2918 oMapping.Add( "ProjectionName", "LambertConformal" );
2919 oMapping.Add( "CenterLongitude", FixLong(
2920 oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)) );
2921 oMapping.Add( "CenterLatitude",
2922 oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0) );
2923 oMapping.Add( "FirstStandardParallel",
2924 oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0) );
2925 oMapping.Add( "SecondStandardParallel",
2926 oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0) );
2927 }
2928
2929 else if( EQUAL(pszProjection, "Vertical Perspective") ) // PROJ 7 required
2930 {
2931 oMapping.Add( "ProjectionName", "PointPerspective" );
2932 oMapping.Add( "CenterLongitude", FixLong(
2933 oSRS.GetNormProjParm("Longitude of topocentric origin", 0.0)) );
2934 oMapping.Add( "CenterLatitude",
2935 oSRS.GetNormProjParm("Latitude of topocentric origin", 0.0) );
2936 // ISIS3 value is the distance from center of ellipsoid, in km
2937 oMapping.Add( "Distance",
2938 (oSRS.GetNormProjParm("Viewpoint height", 0.0) + oSRS.GetSemiMajor()) / 1000.0 );
2939 }
2940
2941 else if( EQUAL(pszProjection, "custom_proj4") )
2942 {
2943 const char* pszProj4 = oSRS.GetExtension("PROJCS", "PROJ4", nullptr);
2944 if( pszProj4 && strstr(pszProj4, "+proj=ob_tran" ) &&
2945 strstr(pszProj4, "+o_proj=eqc") )
2946 {
2947 const auto FetchParam = [](const char* pszProj4Str, const char* pszKey)
2948 {
2949 CPLString needle;
2950 needle.Printf("+%s=", pszKey);
2951 const char* pszVal = strstr(pszProj4Str, needle.c_str());
2952 if( pszVal )
2953 return CPLAtof(pszVal+needle.size());
2954 return 0.0;
2955 };
2956
2957 double dfLonP = FetchParam(pszProj4, "o_lon_p");
2958 double dfLatP = FetchParam(pszProj4, "o_lat_p");
2959 double dfLon0 = FetchParam(pszProj4, "lon_0");
2960 double dfPoleRotation = -dfLonP;
2961 double dfPoleLatitude = 180 - dfLatP;
2962 double dfPoleLongitude = dfLon0;
2963 oMapping.Add( "ProjectionName", "ObliqueCylindrical" );
2964 oMapping.Add( "PoleLatitude", dfPoleLatitude );
2965 oMapping.Add( "PoleLongitude", FixLong(dfPoleLongitude) );
2966 oMapping.Add( "PoleRotation", dfPoleRotation );
2967 }
2968 else
2969 {
2970 CPLError(CE_Warning, CPLE_NotSupported,
2971 "Projection %s not supported",
2972 pszProjection);
2973 }
2974 }
2975 else
2976 {
2977 CPLError(CE_Warning, CPLE_NotSupported,
2978 "Projection %s not supported",
2979 pszProjection);
2980 }
2981
2982
2983 if( oMapping["ProjectionName"].IsValid() )
2984 {
2985 if( oSRS.GetNormProjParm( SRS_PP_FALSE_EASTING, 0.0 ) != 0.0 )
2986 {
2987 CPLError(CE_Warning, CPLE_NotSupported,
2988 "Ignoring %s. Only 0 value supported",
2989 SRS_PP_FALSE_EASTING);
2990 }
2991 if( oSRS.GetNormProjParm( SRS_PP_FALSE_NORTHING, 0.0 ) != 0.0 )
2992 {
2993 CPLError(CE_Warning, CPLE_AppDefined,
2994 "Ignoring %s. Only 0 value supported",
2995 SRS_PP_FALSE_NORTHING);
2996 }
2997 }
2998 }
2999 else
3000 {
3001 CPLError(CE_Warning, CPLE_NotSupported,
3002 "SRS not supported");
3003 }
3004 }
3005
3006 if( !m_bUseSrcMapping && m_bGotTransform )
3007 {
3008 oMapping.Add( "_type", "group" );
3009
3010 const double dfDegToMeter = oSRS.GetSemiMajor() * M_PI / 180.0;
3011 if( !m_oSRS.IsEmpty() && oSRS.IsProjected() )
3012 {
3013 const double dfLinearUnits = oSRS.GetLinearUnits();
3014 // Maybe we should deal differently with non meter units ?
3015 const double dfRes = m_adfGeoTransform[1] * dfLinearUnits;
3016 const double dfScale = dfDegToMeter / dfRes;
3017 oMapping.Add( "UpperLeftCornerX", m_adfGeoTransform[0] );
3018 oMapping.Add( "UpperLeftCornerY", m_adfGeoTransform[3] );
3019 oMapping.Add( "PixelResolution/value", dfRes );
3020 oMapping.Add( "PixelResolution/unit", "meters/pixel" );
3021 oMapping.Add( "Scale/value", dfScale );
3022 oMapping.Add( "Scale/unit", "pixels/degree" );
3023 }
3024 else if( !m_oSRS.IsEmpty() && oSRS.IsGeographic() )
3025 {
3026 const double dfScale = 1.0 / m_adfGeoTransform[1];
3027 const double dfRes = m_adfGeoTransform[1] * dfDegToMeter;
3028 oMapping.Add( "UpperLeftCornerX", m_adfGeoTransform[0] * dfDegToMeter );
3029 oMapping.Add( "UpperLeftCornerY", m_adfGeoTransform[3] * dfDegToMeter );
3030 oMapping.Add( "PixelResolution/value", dfRes );
3031 oMapping.Add( "PixelResolution/unit", "meters/pixel" );
3032 oMapping.Add( "Scale/value", dfScale );
3033 oMapping.Add( "Scale/unit", "pixels/degree" );
3034 }
3035 else
3036 {
3037 oMapping.Add( "UpperLeftCornerX", m_adfGeoTransform[0] );
3038 oMapping.Add( "UpperLeftCornerY", m_adfGeoTransform[3] );
3039 oMapping.Add( "PixelResolution", m_adfGeoTransform[1] );
3040 }
3041 }
3042
3043 CPLJSONObject oLabelLabel = GetOrCreateJSONObject(oLabel, "Label");
3044 oLabelLabel.Set( "_type", "object" );
3045 oLabelLabel.Set( "Bytes", pszLABEL_BYTES_PLACEHOLDER );
3046
3047 // Deal with History object
3048 BuildHistory();
3049
3050 oLabel.Delete( "History" );
3051 if( !m_osHistory.empty() )
3052 {
3053 CPLJSONObject oHistory;
3054 oHistory.Add( "_type", "object" );
3055 oHistory.Add( "Name", "IsisCube" );
3056 if( m_osExternalFilename.empty() )
3057 oHistory.Add( "StartByte", pszHISTORY_STARTBYTE_PLACEHOLDER );
3058 else
3059 oHistory.Add( "StartByte", 1 );
3060 oHistory.Add( "Bytes", static_cast<GIntBig>(m_osHistory.size()) );
3061 if( !m_osExternalFilename.empty() )
3062 {
3063 CPLString osFilename(CPLGetBasename(GetDescription()));
3064 osFilename += ".History.IsisCube";
3065 oHistory.Add( "^History", osFilename );
3066 }
3067 oLabel.Add( "History", oHistory );
3068 }
3069
3070 // Deal with other objects that have StartByte & Bytes
3071 m_aoNonPixelSections.clear();
3072 if( m_oSrcJSonLabel.IsValid() )
3073 {
3074 CPLString osLabelSrcFilename;
3075 CPLJSONObject oFilename = oLabel["_filename"];
3076 if( oFilename.GetType() == CPLJSONObject::Type::String )
3077 {
3078 osLabelSrcFilename = oFilename.ToString();
3079 }
3080
3081 for( CPLJSONObject& oObj : oLabel.GetChildren() )
3082 {
3083 CPLString osKey = oObj.GetName();
3084 if( osKey == "History" )
3085 {
3086 continue;
3087 }
3088
3089 CPLJSONObject oBytes = oObj.GetObj( "Bytes" );
3090 if( oBytes.GetType() != CPLJSONObject::Type::Integer ||
3091 oBytes.ToInteger() <= 0 )
3092 {
3093 continue;
3094 }
3095
3096 CPLJSONObject oStartByte = oObj.GetObj( "StartByte" );
3097 if( oStartByte.GetType() != CPLJSONObject::Type::Integer ||
3098 oStartByte.ToInteger() <= 0 )
3099 {
3100 continue;
3101 }
3102
3103 if( osLabelSrcFilename.empty() )
3104 {
3105 CPLError(CE_Warning, CPLE_AppDefined,
3106 "Cannot find _filename attribute in "
3107 "source ISIS3 metadata. Removing object "
3108 "%s from the label.",
3109 osKey.c_str());
3110 oLabel.Delete( osKey );
3111 continue;
3112 }
3113
3114 NonPixelSection oSection;
3115 oSection.osSrcFilename = osLabelSrcFilename;
3116 oSection.nSrcOffset = static_cast<vsi_l_offset>(
3117 oObj.GetInteger("StartByte")) - 1U;
3118 oSection.nSize = static_cast<vsi_l_offset>(
3119 oObj.GetInteger("Bytes"));
3120
3121 CPLString osName;
3122 CPLJSONObject oName = oObj.GetObj( "Name" );
3123 if( oName.GetType() == CPLJSONObject::Type::String )
3124 {
3125 osName = oName.ToString();
3126 }
3127
3128 CPLString osContainerName(osKey);
3129 CPLJSONObject oContainerName = oObj.GetObj( "_container_name" );
3130 if( oContainerName.GetType() == CPLJSONObject::Type::String )
3131 {
3132 osContainerName = oContainerName.ToString();
3133 }
3134
3135 const CPLString osKeyFilename( "^" + osContainerName );
3136 CPLJSONObject oFilenameCap = oObj.GetObj( osKeyFilename );
3137 if( oFilenameCap.GetType() == CPLJSONObject::Type::String )
3138 {
3139 VSIStatBufL sStat;
3140 const CPLString osSrcFilename( CPLFormFilename(
3141 CPLGetPath(osLabelSrcFilename),
3142 oFilenameCap.ToString().c_str(),
3143 nullptr ) );
3144 if( VSIStatL( osSrcFilename, &sStat ) == 0 )
3145 {
3146 oSection.osSrcFilename = osSrcFilename;
3147 }
3148 else
3149 {
3150 CPLError(CE_Warning, CPLE_AppDefined,
3151 "Object %s points to %s, which does "
3152 "not exist. Removing this section "
3153 "from the label",
3154 osKey.c_str(),
3155 osSrcFilename.c_str());
3156 oLabel.Delete( osKey );
3157 continue;
3158 }
3159 }
3160
3161 if( !m_osExternalFilename.empty() )
3162 {
3163 oObj.Set( "StartByte", 1 );
3164 }
3165 else
3166 {
3167 CPLString osPlaceHolder;
3168 osPlaceHolder.Printf(
3169 "!*^PLACEHOLDER_%d_STARTBYTE^*!",
3170 static_cast<int>(m_aoNonPixelSections.size()) + 1);
3171 oObj.Set( "StartByte", osPlaceHolder );
3172 oSection.osPlaceHolder = osPlaceHolder;
3173 }
3174
3175 if( !m_osExternalFilename.empty() )
3176 {
3177 CPLString osDstFilename( CPLGetBasename(GetDescription()) );
3178 osDstFilename += ".";
3179 osDstFilename += osContainerName;
3180 if( !osName.empty() )
3181 {
3182 osDstFilename += ".";
3183 osDstFilename += osName;
3184 }
3185
3186 oSection.osDstFilename = CPLFormFilename(
3187 CPLGetPath( GetDescription() ),
3188 osDstFilename,
3189 nullptr );
3190
3191 oObj.Set( osKeyFilename, osDstFilename );
3192 }
3193 else
3194 {
3195 oObj.Delete( osKeyFilename );
3196 }
3197
3198 m_aoNonPixelSections.push_back(oSection);
3199 }
3200 }
3201 m_oJSonLabel = oLabel;
3202 }
3203
3204 /************************************************************************/
3205 /* BuildHistory() */
3206 /************************************************************************/
3207
BuildHistory()3208 void ISIS3Dataset::BuildHistory()
3209 {
3210 CPLString osHistory;
3211
3212 if( m_oSrcJSonLabel.IsValid() && m_bUseSrcHistory )
3213 {
3214 vsi_l_offset nHistoryOffset = 0;
3215 int nHistorySize = 0;
3216 CPLString osSrcFilename;
3217
3218 CPLJSONObject oFilename = m_oSrcJSonLabel["_filename"];
3219 if( oFilename.GetType() == CPLJSONObject::Type::String )
3220 {
3221 osSrcFilename = oFilename.ToString();
3222 }
3223 CPLString osHistoryFilename(osSrcFilename);
3224 CPLJSONObject oHistory = m_oSrcJSonLabel["History"];
3225 if( oHistory.GetType() == CPLJSONObject::Type::Object )
3226 {
3227 CPLJSONObject oHistoryFilename = oHistory["^History"];
3228 if( oHistoryFilename.GetType() == CPLJSONObject::Type::String )
3229 {
3230 osHistoryFilename =
3231 CPLFormFilename( CPLGetPath(osSrcFilename),
3232 oHistoryFilename.ToString().c_str(),
3233 nullptr );
3234 }
3235
3236 CPLJSONObject oStartByte = oHistory["StartByte"];
3237 if( oStartByte.GetType() == CPLJSONObject::Type::Integer )
3238 {
3239 if( oStartByte.ToInteger() > 0 )
3240 {
3241 nHistoryOffset = static_cast<vsi_l_offset>(
3242 oStartByte.ToInteger()) - 1U;
3243 }
3244 }
3245
3246 CPLJSONObject oBytes = oHistory["Bytes"];
3247 if( oBytes.GetType() == CPLJSONObject::Type::Integer )
3248 {
3249 nHistorySize = static_cast<int>( oBytes.ToInteger() );
3250 }
3251 }
3252
3253 if( osHistoryFilename.empty() )
3254 {
3255 CPLDebug("ISIS3", "Cannot find filename for source history");
3256 }
3257 else if( nHistorySize <= 0 || nHistorySize > 1000000 )
3258 {
3259 CPLDebug("ISIS3",
3260 "Invalid or missing value for History.Bytes "
3261 "for source history");
3262 }
3263 else
3264 {
3265 VSILFILE* fpHistory = VSIFOpenL(osHistoryFilename, "rb");
3266 if( fpHistory != nullptr )
3267 {
3268 VSIFSeekL(fpHistory, nHistoryOffset, SEEK_SET);
3269 osHistory.resize( nHistorySize );
3270 if( VSIFReadL( &osHistory[0], nHistorySize, 1,
3271 fpHistory ) != 1 )
3272 {
3273 CPLError(CE_Warning, CPLE_FileIO,
3274 "Cannot read %d bytes at offset " CPL_FRMT_GUIB
3275 "of %s: history will not be preserved",
3276 nHistorySize, nHistoryOffset,
3277 osHistoryFilename.c_str());
3278 osHistory.clear();
3279 }
3280 VSIFCloseL(fpHistory);
3281 }
3282 else
3283 {
3284 CPLError(CE_Warning, CPLE_FileIO,
3285 "Cannot open %s: history will not be preserved",
3286 osHistoryFilename.c_str());
3287 }
3288 }
3289 }
3290
3291 if( m_bAddGDALHistory && !m_osGDALHistory.empty() )
3292 {
3293 if( !osHistory.empty() )
3294 osHistory += "\n";
3295 osHistory += m_osGDALHistory;
3296 }
3297 else if( m_bAddGDALHistory )
3298 {
3299 if( !osHistory.empty() )
3300 osHistory += "\n";
3301
3302 CPLJSONObject oHistoryObj;
3303 char szFullFilename[2048] = { 0 };
3304 if( !CPLGetExecPath(szFullFilename, sizeof(szFullFilename) - 1) )
3305 strcpy(szFullFilename, "unknown_program");
3306 const CPLString osProgram(CPLGetBasename(szFullFilename));
3307 const CPLString osPath(CPLGetPath(szFullFilename));
3308
3309 CPLJSONObject oObj;
3310 oHistoryObj.Add( osProgram, oObj );
3311
3312 oObj.Add( "_type", "object" );
3313 oObj.Add( "GdalVersion", GDALVersionInfo("RELEASE_NAME") );
3314 if( osPath != "." )
3315 oObj.Add( "ProgramPath", osPath );
3316 time_t nCurTime = time(nullptr);
3317 if( nCurTime != -1 )
3318 {
3319 struct tm mytm;
3320 CPLUnixTimeToYMDHMS(nCurTime, &mytm);
3321 oObj.Add( "ExecutionDateTime",
3322 CPLSPrintf("%04d-%02d-%02dT%02d:%02d:%02d",
3323 mytm.tm_year + 1900,
3324 mytm.tm_mon + 1,
3325 mytm.tm_mday,
3326 mytm.tm_hour,
3327 mytm.tm_min,
3328 mytm.tm_sec) );
3329 }
3330 char szHostname[256] = { 0 };
3331 if( gethostname(szHostname, sizeof(szHostname)-1) == 0 )
3332 {
3333 oObj.Add( "HostName", std::string(szHostname) );
3334 }
3335 const char* pszUsername = CPLGetConfigOption("USERNAME", nullptr);
3336 if( pszUsername == nullptr )
3337 pszUsername = CPLGetConfigOption("USER", nullptr);
3338 if( pszUsername != nullptr )
3339 {
3340 oObj.Add( "UserName", pszUsername );
3341 }
3342 oObj.Add( "Description", "GDAL conversion" );
3343
3344 CPLJSONObject oUserParameters;
3345 oObj.Add( "UserParameters", oUserParameters );
3346
3347 oUserParameters.Add( "_type", "group");
3348 if( !m_osFromFilename.empty() )
3349 {
3350 const CPLString osFromFilename = CPLGetFilename( m_osFromFilename );
3351 oUserParameters.Add( "FROM", osFromFilename );
3352 }
3353 if( nullptr != GetDescription() )
3354 {
3355 const CPLString osToFileName = CPLGetFilename( GetDescription() );
3356 oUserParameters.Add( "TO", osToFileName );
3357 }
3358 if( m_bForce360 )
3359 oUserParameters.Add( "Force_360", "true");
3360
3361 osHistory += SerializeAsPDL( oHistoryObj );
3362 }
3363
3364 m_osHistory = osHistory;
3365 }
3366
3367 /************************************************************************/
3368 /* WriteLabel() */
3369 /************************************************************************/
3370
WriteLabel()3371 void ISIS3Dataset::WriteLabel()
3372 {
3373 m_bIsLabelWritten = true;
3374
3375 if( !m_oJSonLabel.IsValid() )
3376 BuildLabel();
3377
3378 // Serialize label
3379 CPLString osLabel( SerializeAsPDL(m_oJSonLabel) );
3380 osLabel += "End\n";
3381 if( m_osExternalFilename.empty() && osLabel.size() < 65536 )
3382 {
3383 // In-line labels have conventionally a minimize size of 65536 bytes
3384 // See #2741
3385 osLabel.resize(65536);
3386 }
3387 char *pszLabel = &osLabel[0];
3388 const int nLabelSize = static_cast<int>(osLabel.size());
3389
3390 // Hack back StartByte value
3391 {
3392 char *pszStartByte = strstr(pszLabel, pszSTARTBYTE_PLACEHOLDER);
3393 if( pszStartByte != nullptr )
3394 {
3395 const char* pszOffset = CPLSPrintf("%d", 1 + nLabelSize);
3396 memcpy(pszStartByte, pszOffset, strlen(pszOffset));
3397 memset(pszStartByte + strlen(pszOffset), ' ',
3398 strlen(pszSTARTBYTE_PLACEHOLDER) - strlen(pszOffset));
3399 }
3400 }
3401
3402 // Hack back Label.Bytes value
3403 {
3404 char* pszLabelBytes = strstr(pszLabel, pszLABEL_BYTES_PLACEHOLDER);
3405 if( pszLabelBytes != nullptr )
3406 {
3407 const char* pszBytes = CPLSPrintf("%d", nLabelSize);
3408 memcpy(pszLabelBytes, pszBytes, strlen(pszBytes));
3409 memset(pszLabelBytes + strlen(pszBytes), ' ',
3410 strlen(pszLABEL_BYTES_PLACEHOLDER) - strlen(pszBytes));
3411 }
3412 }
3413
3414 const GDALDataType eType = GetRasterBand(1)->GetRasterDataType();
3415 const int nDTSize = GDALGetDataTypeSizeBytes(eType);
3416 vsi_l_offset nImagePixels = 0;
3417 if( m_poExternalDS == nullptr )
3418 {
3419 if( m_bIsTiled )
3420 {
3421 int nBlockXSize = 1, nBlockYSize = 1;
3422 GetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
3423 nImagePixels = static_cast<vsi_l_offset>(nBlockXSize) *
3424 nBlockYSize * nBands *
3425 DIV_ROUND_UP(nRasterXSize, nBlockXSize) *
3426 DIV_ROUND_UP(nRasterYSize, nBlockYSize);
3427 }
3428 else
3429 {
3430 nImagePixels = static_cast<vsi_l_offset>(nRasterXSize) *
3431 nRasterYSize * nBands;
3432 }
3433 }
3434
3435 // Hack back History.StartBytes value
3436 char* pszHistoryStartBytes = strstr(pszLabel,
3437 pszHISTORY_STARTBYTE_PLACEHOLDER);
3438
3439 vsi_l_offset nHistoryOffset = 0;
3440 vsi_l_offset nLastOffset = 0;
3441 if( pszHistoryStartBytes != nullptr )
3442 {
3443 CPLAssert( m_osExternalFilename.empty() );
3444 nHistoryOffset = nLabelSize + nImagePixels * nDTSize;
3445 nLastOffset = nHistoryOffset + m_osHistory.size();
3446 const char* pszStartByte = CPLSPrintf(CPL_FRMT_GUIB,
3447 nHistoryOffset + 1);
3448 CPLAssert(strlen(pszStartByte) <
3449 strlen(pszHISTORY_STARTBYTE_PLACEHOLDER));
3450 memcpy(pszHistoryStartBytes, pszStartByte, strlen(pszStartByte));
3451 memset(pszHistoryStartBytes + strlen(pszStartByte), ' ',
3452 strlen(pszHISTORY_STARTBYTE_PLACEHOLDER) - strlen(pszStartByte));
3453 }
3454
3455 // Replace placeholders in other sections
3456 for( size_t i = 0; i < m_aoNonPixelSections.size(); ++i )
3457 {
3458 if( !m_aoNonPixelSections[i].osPlaceHolder.empty() )
3459 {
3460 char* pszPlaceHolder = strstr(pszLabel,
3461 m_aoNonPixelSections[i].osPlaceHolder.c_str());
3462 CPLAssert( pszPlaceHolder != nullptr );
3463 const char* pszStartByte = CPLSPrintf(CPL_FRMT_GUIB,
3464 nLastOffset + 1);
3465 nLastOffset += m_aoNonPixelSections[i].nSize;
3466 CPLAssert(strlen(pszStartByte) <
3467 m_aoNonPixelSections[i].osPlaceHolder.size() );
3468
3469 memcpy(pszPlaceHolder, pszStartByte, strlen(pszStartByte));
3470 memset(pszPlaceHolder + strlen(pszStartByte), ' ',
3471 m_aoNonPixelSections[i].osPlaceHolder.size() -
3472 strlen(pszStartByte));
3473 }
3474 }
3475
3476 // Write to final file
3477 VSIFSeekL( m_fpLabel, 0, SEEK_SET );
3478 VSIFWriteL( pszLabel, 1, osLabel.size(), m_fpLabel);
3479
3480 if( m_osExternalFilename.empty() )
3481 {
3482 // Update image offset in bands
3483 if( m_bIsTiled )
3484 {
3485 for(int i=0;i<nBands;i++)
3486 {
3487 ISISTiledBand* poBand =
3488 reinterpret_cast<ISISTiledBand*>(GetRasterBand(i+1));
3489 poBand->m_nFirstTileOffset += nLabelSize;
3490 }
3491 }
3492 else
3493 {
3494 for(int i=0;i<nBands;i++)
3495 {
3496 ISIS3RawRasterBand* poBand =
3497 reinterpret_cast<ISIS3RawRasterBand*>(GetRasterBand(i+1));
3498 poBand->nImgOffset += nLabelSize;
3499 }
3500 }
3501 }
3502
3503 if( m_bInitToNodata )
3504 {
3505 // Initialize the image to nodata
3506 const double dfNoData = GetRasterBand(1)->GetNoDataValue();
3507 if( dfNoData == 0.0 )
3508 {
3509 VSIFTruncateL( m_fpImage, VSIFTellL(m_fpImage) +
3510 nImagePixels * nDTSize );
3511 }
3512 else if( nDTSize != 0 ) // to make Coverity not warn about div by 0
3513 {
3514 const int nPageSize = 4096; // Must be multiple of 4 since
3515 // Float32 is the largest type
3516 CPLAssert( (nPageSize % nDTSize) == 0 );
3517 const int nMaxPerPage = nPageSize / nDTSize;
3518 GByte* pabyTemp = static_cast<GByte*>(CPLMalloc(nPageSize));
3519 GDALCopyWords( &dfNoData, GDT_Float64, 0,
3520 pabyTemp, eType, nDTSize,
3521 nMaxPerPage );
3522 #ifdef CPL_MSB
3523 GDALSwapWords( pabyTemp, nDTSize, nMaxPerPage, nDTSize );
3524 #endif
3525 for( vsi_l_offset i = 0; i < nImagePixels; i += nMaxPerPage )
3526 {
3527 int n;
3528 if( i + nMaxPerPage <= nImagePixels )
3529 n = nMaxPerPage;
3530 else
3531 n = static_cast<int>(nImagePixels - i);
3532 if( VSIFWriteL( pabyTemp, n * nDTSize, 1, m_fpImage ) != 1 )
3533 {
3534 CPLError(CE_Failure, CPLE_FileIO,
3535 "Cannot initialize imagery to null");
3536 break;
3537 }
3538 }
3539
3540 CPLFree( pabyTemp );
3541 }
3542 }
3543
3544 // Write history
3545 if( !m_osHistory.empty() )
3546 {
3547 if( m_osExternalFilename.empty() )
3548 {
3549 VSIFSeekL( m_fpLabel, nHistoryOffset, SEEK_SET );
3550 VSIFWriteL( m_osHistory.c_str(), 1, m_osHistory.size(), m_fpLabel);
3551 }
3552 else
3553 {
3554 CPLString osFilename(CPLGetBasename(GetDescription()));
3555 osFilename += ".History.IsisCube";
3556 osFilename =
3557 CPLFormFilename(CPLGetPath(GetDescription()), osFilename, nullptr);
3558 VSILFILE* fp = VSIFOpenL(osFilename, "wb");
3559 if( fp )
3560 {
3561 m_aosAdditionalFiles.AddString(osFilename);
3562
3563 VSIFWriteL( m_osHistory.c_str(), 1,
3564 m_osHistory.size(), fp );
3565 VSIFCloseL(fp);
3566 }
3567 else
3568 {
3569 CPLError(CE_Warning, CPLE_FileIO,
3570 "Cannot write %s", osFilename.c_str());
3571 }
3572 }
3573 }
3574
3575 // Write other non pixel sections
3576 for( size_t i = 0; i < m_aoNonPixelSections.size(); ++i )
3577 {
3578 VSILFILE* fpSrc = VSIFOpenL(
3579 m_aoNonPixelSections[i].osSrcFilename, "rb");
3580 if( fpSrc == nullptr )
3581 {
3582 CPLError(CE_Warning, CPLE_FileIO,
3583 "Cannot open %s",
3584 m_aoNonPixelSections[i].osSrcFilename.c_str());
3585 continue;
3586 }
3587
3588 VSILFILE* fpDest = m_fpLabel;
3589 if( !m_aoNonPixelSections[i].osDstFilename.empty() )
3590 {
3591 fpDest = VSIFOpenL(m_aoNonPixelSections[i].osDstFilename, "wb");
3592 if( fpDest == nullptr )
3593 {
3594 CPLError(CE_Warning, CPLE_FileIO,
3595 "Cannot create %s",
3596 m_aoNonPixelSections[i].osDstFilename.c_str());
3597 VSIFCloseL(fpSrc);
3598 continue;
3599 }
3600
3601 m_aosAdditionalFiles.AddString(
3602 m_aoNonPixelSections[i].osDstFilename);
3603 }
3604
3605 VSIFSeekL(fpSrc, m_aoNonPixelSections[i].nSrcOffset, SEEK_SET);
3606 GByte abyBuffer[4096];
3607 vsi_l_offset nRemaining = m_aoNonPixelSections[i].nSize;
3608 while( nRemaining )
3609 {
3610 size_t nToRead = 4096;
3611 if( nRemaining < nToRead )
3612 nToRead = static_cast<size_t>(nRemaining);
3613 size_t nRead = VSIFReadL( abyBuffer, 1, nToRead, fpSrc );
3614 if( nRead != nToRead )
3615 {
3616 CPLError(CE_Warning, CPLE_FileIO,
3617 "Could not read " CPL_FRMT_GUIB " bytes from %s",
3618 m_aoNonPixelSections[i].nSize,
3619 m_aoNonPixelSections[i].osSrcFilename.c_str());
3620 break;
3621 }
3622 VSIFWriteL( abyBuffer, 1, nRead, fpDest );
3623 nRemaining -= nRead;
3624 }
3625
3626 VSIFCloseL( fpSrc );
3627 if( fpDest != m_fpLabel )
3628 VSIFCloseL(fpDest);
3629 }
3630 }
3631
3632 /************************************************************************/
3633 /* SerializeAsPDL() */
3634 /************************************************************************/
3635
SerializeAsPDL(const CPLJSONObject & oObj)3636 CPLString ISIS3Dataset::SerializeAsPDL( const CPLJSONObject &oObj )
3637 {
3638 CPLString osTmpFile( CPLSPrintf("/vsimem/isis3_%p", oObj.GetInternalHandle()) );
3639 VSILFILE* fpTmp = VSIFOpenL( osTmpFile, "wb+" );
3640 SerializeAsPDL( fpTmp, oObj );
3641 VSIFCloseL( fpTmp );
3642 CPLString osContent( reinterpret_cast<char*>(
3643 VSIGetMemFileBuffer( osTmpFile, nullptr, FALSE )) );
3644 VSIUnlink(osTmpFile);
3645 return osContent;
3646 }
3647
3648 /************************************************************************/
3649 /* SerializeAsPDL() */
3650 /************************************************************************/
3651
SerializeAsPDL(VSILFILE * fp,const CPLJSONObject & oObj,int nDepth)3652 void ISIS3Dataset::SerializeAsPDL( VSILFILE* fp, const CPLJSONObject &oObj,
3653 int nDepth )
3654 {
3655 CPLString osIndentation;
3656 for( int i = 0; i < nDepth; i++ )
3657 osIndentation += " ";
3658 const size_t WIDTH = 79;
3659
3660 std::vector<CPLJSONObject> aoChildren = oObj.GetChildren();
3661 size_t nMaxKeyLength = 0;
3662 for( const CPLJSONObject& oChild : aoChildren )
3663 {
3664 const CPLString osKey = oChild.GetName();
3665 if( EQUAL(osKey, "_type") ||
3666 EQUAL(osKey, "_container_name") ||
3667 EQUAL(osKey, "_filename") )
3668 {
3669 continue;
3670 }
3671
3672 const auto eType = oChild.GetType();
3673 if( eType == CPLJSONObject::Type::String ||
3674 eType == CPLJSONObject::Type::Integer ||
3675 eType == CPLJSONObject::Type::Double ||
3676 eType == CPLJSONObject::Type::Array )
3677 {
3678 if( osKey.size() > nMaxKeyLength )
3679 {
3680 nMaxKeyLength = osKey.size();
3681 }
3682 }
3683 else if( eType == CPLJSONObject::Type::Object )
3684 {
3685 CPLJSONObject oValue = oChild.GetObj( "value" );
3686 CPLJSONObject oUnit = oChild.GetObj( "unit" );
3687 if( oValue.IsValid() && oUnit.GetType() == CPLJSONObject::Type::String )
3688 {
3689 if( osKey.size() > nMaxKeyLength )
3690 {
3691 nMaxKeyLength = osKey.size();
3692 }
3693 }
3694 }
3695 }
3696
3697 for( const CPLJSONObject& oChild : aoChildren )
3698 {
3699 const CPLString osKey = oChild.GetName();
3700 if( EQUAL(osKey, "_type") ||
3701 EQUAL(osKey, "_container_name") ||
3702 EQUAL(osKey, "_filename") )
3703 {
3704 continue;
3705 }
3706 if( STARTS_WITH(osKey, "_comment") )
3707 {
3708 if( oChild.GetType() == CPLJSONObject::Type::String )
3709 {
3710 VSIFPrintfL(fp, "#%s\n", oChild.ToString().c_str() );
3711 }
3712 continue;
3713 }
3714 CPLString osPadding;
3715 size_t nLen = osKey.size();
3716 if( nLen < nMaxKeyLength )
3717 {
3718 osPadding.append( nMaxKeyLength - nLen, ' ' );
3719 }
3720
3721 const auto eType = oChild.GetType();
3722 if( eType == CPLJSONObject::Type::Object )
3723 {
3724 CPLJSONObject oType = oChild.GetObj( "_type" );
3725 CPLJSONObject oContainerName = oChild.GetObj( "_container_name" );
3726 CPLString osContainerName = osKey;
3727 if( oContainerName.GetType() == CPLJSONObject::Type::String )
3728 {
3729 osContainerName = oContainerName.ToString();
3730 }
3731 if( oType.GetType() == CPLJSONObject::Type::String )
3732 {
3733 const CPLString osType = oType.ToString();
3734 if( EQUAL(osType, "Object") )
3735 {
3736 if( nDepth == 0 && VSIFTellL(fp) != 0 )
3737 VSIFPrintfL(fp, "\n");
3738 VSIFPrintfL(fp, "%sObject = %s\n",
3739 osIndentation.c_str(), osContainerName.c_str());
3740 SerializeAsPDL( fp, oChild, nDepth + 1 );
3741 VSIFPrintfL(fp, "%sEnd_Object\n", osIndentation.c_str());
3742 }
3743 else if( EQUAL(osType, "Group") )
3744 {
3745 VSIFPrintfL(fp, "\n");
3746 VSIFPrintfL(fp, "%sGroup = %s\n",
3747 osIndentation.c_str(), osContainerName.c_str());
3748 SerializeAsPDL( fp, oChild, nDepth + 1 );
3749 VSIFPrintfL(fp, "%sEnd_Group\n", osIndentation.c_str());
3750 }
3751 }
3752 else
3753 {
3754 CPLJSONObject oValue = oChild.GetObj( "value" );
3755 CPLJSONObject oUnit = oChild.GetObj( "unit" );
3756 if( oValue.IsValid() &&
3757 oUnit.GetType() == CPLJSONObject::Type::String )
3758 {
3759 const CPLString osUnit = oUnit.ToString();
3760 const auto eValueType = oValue.GetType();
3761 if( eValueType == CPLJSONObject::Type::Integer )
3762 {
3763 VSIFPrintfL(fp, "%s%s%s = %d <%s>\n",
3764 osIndentation.c_str(), osKey.c_str(),
3765 osPadding.c_str(),
3766 oValue.ToInteger(), osUnit.c_str());
3767 }
3768 else if( eValueType == CPLJSONObject::Type::Double )
3769 {
3770 const double dfVal = oValue.ToDouble();
3771 if( dfVal >= INT_MIN && dfVal <= INT_MAX &&
3772 static_cast<int>(dfVal) == dfVal )
3773 {
3774 VSIFPrintfL(fp, "%s%s%s = %d.0 <%s>\n",
3775 osIndentation.c_str(), osKey.c_str(),
3776 osPadding.c_str(),
3777 static_cast<int>(dfVal), osUnit.c_str());
3778 }
3779 else
3780 {
3781 VSIFPrintfL(fp, "%s%s%s = %.18g <%s>\n",
3782 osIndentation.c_str(), osKey.c_str(),
3783 osPadding.c_str(),
3784 dfVal, osUnit.c_str());
3785 }
3786 }
3787 }
3788 }
3789 }
3790 else if( eType == CPLJSONObject::Type::String )
3791 {
3792 CPLString osVal = oChild.ToString();
3793 const char* pszVal = osVal.c_str();
3794 if( pszVal[0] == '\0' ||
3795 strchr(pszVal, ' ') || strstr(pszVal, "\\n") ||
3796 strstr(pszVal, "\\r") )
3797 {
3798 osVal.replaceAll("\\n", "\n");
3799 osVal.replaceAll("\\r", "\r");
3800 VSIFPrintfL(fp, "%s%s%s = \"%s\"\n",
3801 osIndentation.c_str(), osKey.c_str(),
3802 osPadding.c_str(), osVal.c_str());
3803 }
3804 else
3805 {
3806 if( osIndentation.size() + osKey.size() + osPadding.size() +
3807 strlen(" = ") + strlen(pszVal) > WIDTH &&
3808 osIndentation.size() + osKey.size() + osPadding.size() +
3809 strlen(" = ") < WIDTH )
3810 {
3811 size_t nFirstPos = osIndentation.size() + osKey.size() +
3812 osPadding.size() + strlen(" = ");
3813 VSIFPrintfL(fp, "%s%s%s = ",
3814 osIndentation.c_str(), osKey.c_str(),
3815 osPadding.c_str());
3816 size_t nCurPos = nFirstPos;
3817 for( int j = 0; pszVal[j] != '\0'; j++ )
3818 {
3819 nCurPos ++;
3820 if( nCurPos == WIDTH && pszVal[j+1] != '\0' )
3821 {
3822 VSIFPrintfL( fp, "-\n" );
3823 for( size_t k=0;k<nFirstPos;k++ )
3824 {
3825 const char chSpace = ' ';
3826 VSIFWriteL(&chSpace, 1, 1, fp);
3827 }
3828 nCurPos = nFirstPos + 1;
3829 }
3830 VSIFWriteL( &pszVal[j], 1, 1, fp );
3831 }
3832 VSIFPrintfL(fp, "\n");
3833 }
3834 else
3835 {
3836 VSIFPrintfL(fp, "%s%s%s = %s\n",
3837 osIndentation.c_str(), osKey.c_str(),
3838 osPadding.c_str(), pszVal);
3839 }
3840 }
3841 }
3842 else if( eType == CPLJSONObject::Type::Integer )
3843 {
3844 const int nVal = oChild.ToInteger();
3845 VSIFPrintfL(fp, "%s%s%s = %d\n",
3846 osIndentation.c_str(), osKey.c_str(),
3847 osPadding.c_str(), nVal);
3848 }
3849 else if( eType == CPLJSONObject::Type::Double )
3850 {
3851 const double dfVal = oChild.ToDouble();
3852 if( dfVal >= INT_MIN && dfVal <= INT_MAX &&
3853 static_cast<int>(dfVal) == dfVal )
3854 {
3855 VSIFPrintfL(fp, "%s%s%s = %d.0\n",
3856 osIndentation.c_str(), osKey.c_str(),
3857 osPadding.c_str(), static_cast<int>(dfVal));
3858 }
3859 else
3860 {
3861 VSIFPrintfL(fp, "%s%s%s = %.18g\n",
3862 osIndentation.c_str(), osKey.c_str(),
3863 osPadding.c_str(), dfVal);
3864 }
3865 }
3866 else if( eType == CPLJSONObject::Type::Array )
3867 {
3868 CPLJSONArray oArrayItem(oChild);
3869 const int nLength = oArrayItem.Size();
3870 size_t nFirstPos = osIndentation.size() + osKey.size() +
3871 osPadding.size() + strlen(" = (");
3872 VSIFPrintfL(fp, "%s%s%s = (",
3873 osIndentation.c_str(), osKey.c_str(),
3874 osPadding.c_str());
3875 size_t nCurPos = nFirstPos;
3876 for( int idx = 0; idx < nLength; idx++ )
3877 {
3878 CPLJSONObject oItem = oArrayItem[idx];
3879 const auto eArrayItemType = oItem.GetType();
3880 if( eArrayItemType == CPLJSONObject::Type::String )
3881 {
3882 CPLString osVal = oItem.ToString();
3883 const char* pszVal = osVal.c_str();
3884 if( pszVal[0] == '\0' ||
3885 strchr(pszVal, ' ') || strstr(pszVal, "\\n") ||
3886 strstr(pszVal, "\\r") )
3887 {
3888 osVal.replaceAll("\\n", "\n");
3889 osVal.replaceAll("\\r", "\r");
3890 VSIFPrintfL(fp, "\"%s\"", osVal.c_str());
3891 }
3892 else if( nFirstPos < WIDTH && nCurPos + strlen(pszVal) > WIDTH )
3893 {
3894 if( idx > 0 )
3895 {
3896 VSIFPrintfL( fp, "\n" );
3897 for( size_t j=0;j<nFirstPos;j++ )
3898 {
3899 const char chSpace = ' ';
3900 VSIFWriteL(&chSpace, 1, 1, fp);
3901 }
3902 nCurPos = nFirstPos;
3903 }
3904
3905 for( int j = 0; pszVal[j] != '\0'; j++ )
3906 {
3907 nCurPos ++;
3908 if( nCurPos == WIDTH && pszVal[j+1] != '\0' )
3909 {
3910 VSIFPrintfL( fp, "-\n" );
3911 for( size_t k=0;k<nFirstPos;k++ )
3912 {
3913 const char chSpace = ' ';
3914 VSIFWriteL(&chSpace, 1, 1, fp);
3915 }
3916 nCurPos = nFirstPos + 1;
3917 }
3918 VSIFWriteL( &pszVal[j], 1, 1, fp );
3919 }
3920 }
3921 else
3922 {
3923 VSIFPrintfL( fp, "%s", pszVal );
3924 nCurPos += strlen(pszVal);
3925 }
3926 }
3927 else if( eArrayItemType == CPLJSONObject::Type::Integer )
3928 {
3929 const int nVal = oItem.ToInteger();
3930 const char* pszVal = CPLSPrintf("%d", nVal);
3931 const size_t nValLen = strlen(pszVal);
3932 if( nFirstPos < WIDTH && idx > 0 &&
3933 nCurPos + nValLen > WIDTH )
3934 {
3935 VSIFPrintfL( fp, "\n" );
3936 for( size_t j=0;j<nFirstPos;j++ )
3937 {
3938 const char chSpace = ' ';
3939 VSIFWriteL(&chSpace, 1, 1, fp);
3940 }
3941 nCurPos = nFirstPos;
3942 }
3943 VSIFPrintfL( fp, "%d", nVal );
3944 nCurPos += nValLen;
3945 }
3946 else if( eArrayItemType == CPLJSONObject::Type::Double )
3947 {
3948 const double dfVal = oItem.ToDouble();
3949 CPLString osVal;
3950 if( dfVal >= INT_MIN && dfVal <= INT_MAX &&
3951 static_cast<int>(dfVal) == dfVal )
3952 {
3953 osVal = CPLSPrintf("%d.0", static_cast<int>(dfVal));
3954 }
3955 else
3956 {
3957 osVal = CPLSPrintf("%.18g", dfVal);
3958 }
3959 const size_t nValLen = osVal.size();
3960 if( nFirstPos < WIDTH && idx > 0 &&
3961 nCurPos + nValLen > WIDTH )
3962 {
3963 VSIFPrintfL( fp, "\n" );
3964 for( size_t j=0;j<nFirstPos;j++ )
3965 {
3966 const char chSpace = ' ';
3967 VSIFWriteL(&chSpace, 1, 1, fp);
3968 }
3969 nCurPos = nFirstPos;
3970 }
3971 VSIFPrintfL( fp, "%s", osVal.c_str() );
3972 nCurPos += nValLen;
3973 }
3974 if( idx < nLength - 1 )
3975 {
3976 VSIFPrintfL( fp, ", " );
3977 nCurPos += 2;
3978 }
3979 }
3980 VSIFPrintfL(fp, ")\n" );
3981 }
3982 }
3983 }
3984
3985 /************************************************************************/
3986 /* Create() */
3987 /************************************************************************/
3988
Create(const char * pszFilename,int nXSize,int nYSize,int nBands,GDALDataType eType,char ** papszOptions)3989 GDALDataset *ISIS3Dataset::Create(const char* pszFilename,
3990 int nXSize, int nYSize, int nBands,
3991 GDALDataType eType,
3992 char** papszOptions)
3993 {
3994 if( eType != GDT_Byte && eType != GDT_UInt16 && eType != GDT_Int16 &&
3995 eType != GDT_Float32 )
3996 {
3997 CPLError(CE_Failure, CPLE_NotSupported, "Unsupported data type");
3998 return nullptr;
3999 }
4000 if( nBands == 0 || nBands > 32767 )
4001 {
4002 CPLError(CE_Failure, CPLE_NotSupported, "Unsupported band count");
4003 return nullptr;
4004 }
4005
4006 const char* pszDataLocation = CSLFetchNameValueDef(papszOptions,
4007 "DATA_LOCATION",
4008 "LABEL");
4009 const bool bIsTiled = CPLFetchBool(papszOptions, "TILED", false);
4010 const int nBlockXSize = std::max(1,
4011 atoi(CSLFetchNameValueDef(papszOptions, "BLOCKXSIZE", "256")));
4012 const int nBlockYSize = std::max(1,
4013 atoi(CSLFetchNameValueDef(papszOptions, "BLOCKYSIZE", "256")));
4014 if( !EQUAL(pszDataLocation, "LABEL") &&
4015 !EQUAL( CPLGetExtension(pszFilename), "LBL") )
4016 {
4017 CPLError(CE_Failure, CPLE_NotSupported,
4018 "For DATA_LOCATION=%s, "
4019 "the main filename should have a .lbl extension",
4020 pszDataLocation);
4021 return nullptr;
4022 }
4023
4024 VSILFILE* fp = VSIFOpenExL(pszFilename, "wb", true);
4025 if( fp == nullptr )
4026 {
4027 CPLError( CE_Failure, CPLE_FileIO,
4028 "Cannot create %s: %s",
4029 pszFilename, VSIGetLastErrorMsg() );
4030 return nullptr;
4031 }
4032 VSILFILE* fpImage = nullptr;
4033 CPLString osExternalFilename;
4034 GDALDataset* poExternalDS = nullptr;
4035 bool bGeoTIFFAsRegularExternal = false;
4036 if( EQUAL(pszDataLocation, "EXTERNAL") )
4037 {
4038 osExternalFilename = CSLFetchNameValueDef(papszOptions,
4039 "EXTERNAL_FILENAME",
4040 CPLResetExtension(pszFilename, "cub"));
4041 fpImage = VSIFOpenExL(osExternalFilename, "wb", true);
4042 if( fpImage == nullptr )
4043 {
4044 CPLError( CE_Failure, CPLE_FileIO,
4045 "Cannot create %s: %s",
4046 osExternalFilename.c_str(), VSIGetLastErrorMsg() );
4047 VSIFCloseL(fp);
4048 return nullptr;
4049 }
4050 }
4051 else if( EQUAL(pszDataLocation, "GEOTIFF") )
4052 {
4053 osExternalFilename = CSLFetchNameValueDef(papszOptions,
4054 "EXTERNAL_FILENAME",
4055 CPLResetExtension(pszFilename, "tif"));
4056 GDALDriver* poDrv = static_cast<GDALDriver*>(
4057 GDALGetDriverByName("GTiff"));
4058 if( poDrv == nullptr )
4059 {
4060 CPLError( CE_Failure, CPLE_AppDefined,
4061 "Cannot find GTiff driver" );
4062 VSIFCloseL(fp);
4063 return nullptr;
4064 }
4065 char** papszGTiffOptions = nullptr;
4066 papszGTiffOptions = CSLSetNameValue(papszGTiffOptions,
4067 "ENDIANNESS", "LITTLE");
4068 if( bIsTiled )
4069 {
4070 papszGTiffOptions = CSLSetNameValue(papszGTiffOptions,
4071 "TILED", "YES");
4072 papszGTiffOptions = CSLSetNameValue(papszGTiffOptions,
4073 "BLOCKXSIZE",
4074 CPLSPrintf("%d", nBlockXSize));
4075 papszGTiffOptions = CSLSetNameValue(papszGTiffOptions,
4076 "BLOCKYSIZE",
4077 CPLSPrintf("%d", nBlockYSize));
4078 }
4079 const char* pszGTiffOptions = CSLFetchNameValueDef(papszOptions,
4080 "GEOTIFF_OPTIONS", "");
4081 char** papszTokens = CSLTokenizeString2( pszGTiffOptions, ",", 0 );
4082 for( int i = 0; papszTokens[i] != nullptr; i++ )
4083 {
4084 papszGTiffOptions = CSLAddString(papszGTiffOptions,
4085 papszTokens[i]);
4086 }
4087 CSLDestroy(papszTokens);
4088
4089 // If the user didn't specify any compression and
4090 // GEOTIFF_AS_REGULAR_EXTERNAL is set (or unspecified), then the
4091 // GeoTIFF file can be seen as a regular external raw file, provided
4092 // we make some provision on its organization.
4093 if( CSLFetchNameValue(papszGTiffOptions, "COMPRESS") == nullptr &&
4094 CPLFetchBool(papszOptions, "GEOTIFF_AS_REGULAR_EXTERNAL", true) )
4095 {
4096 bGeoTIFFAsRegularExternal = true;
4097 papszGTiffOptions = CSLSetNameValue(papszGTiffOptions,
4098 "INTERLEAVE", "BAND");
4099 // Will make sure that our blocks at nodata are not optimized
4100 // away but indeed well written
4101 papszGTiffOptions = CSLSetNameValue(papszGTiffOptions,
4102 "@WRITE_EMPTY_TILES_SYNCHRONOUSLY", "YES");
4103 if( !bIsTiled && nBands > 1 )
4104 {
4105 papszGTiffOptions = CSLSetNameValue(papszGTiffOptions,
4106 "BLOCKYSIZE", "1");
4107 }
4108 }
4109
4110 poExternalDS = poDrv->Create( osExternalFilename, nXSize, nYSize,
4111 nBands,
4112 eType, papszGTiffOptions );
4113 CSLDestroy(papszGTiffOptions);
4114 if( poExternalDS == nullptr )
4115 {
4116 CPLError( CE_Failure, CPLE_FileIO,
4117 "Cannot create %s",
4118 osExternalFilename.c_str() );
4119 VSIFCloseL(fp);
4120 return nullptr;
4121 }
4122 }
4123
4124 ISIS3Dataset* poDS = new ISIS3Dataset();
4125 poDS->SetDescription( pszFilename );
4126 poDS->eAccess = GA_Update;
4127 poDS->nRasterXSize = nXSize;
4128 poDS->nRasterYSize = nYSize;
4129 poDS->m_osExternalFilename = osExternalFilename;
4130 poDS->m_poExternalDS = poExternalDS;
4131 poDS->m_bGeoTIFFAsRegularExternal = bGeoTIFFAsRegularExternal;
4132 if( bGeoTIFFAsRegularExternal )
4133 poDS->m_bGeoTIFFInitDone = false;
4134 poDS->m_fpLabel = fp;
4135 poDS->m_fpImage = fpImage ? fpImage: fp;
4136 poDS->m_bIsLabelWritten = false;
4137 poDS->m_bIsTiled = bIsTiled;
4138 poDS->m_bInitToNodata = (poDS->m_poExternalDS == nullptr);
4139 poDS->m_osComment = CSLFetchNameValueDef(papszOptions, "COMMENT", "");
4140 poDS->m_osLatitudeType = CSLFetchNameValueDef(papszOptions,
4141 "LATITUDE_TYPE", "");
4142 poDS->m_osLongitudeDirection = CSLFetchNameValueDef(papszOptions,
4143 "LONGITUDE_DIRECTION", "");
4144 poDS->m_osTargetName = CSLFetchNameValueDef(papszOptions,
4145 "TARGET_NAME", "");
4146 poDS->m_bForce360 = CPLFetchBool(papszOptions, "FORCE_360", false);
4147 poDS->m_bWriteBoundingDegrees = CPLFetchBool(papszOptions,
4148 "WRITE_BOUNDING_DEGREES",
4149 true);
4150 poDS->m_osBoundingDegrees = CSLFetchNameValueDef(papszOptions,
4151 "BOUNDING_DEGREES", "");
4152 poDS->m_bUseSrcLabel = CPLFetchBool(papszOptions, "USE_SRC_LABEL", true);
4153 poDS->m_bUseSrcMapping =
4154 CPLFetchBool(papszOptions, "USE_SRC_MAPPING", false);
4155 poDS->m_bUseSrcHistory =
4156 CPLFetchBool(papszOptions, "USE_SRC_HISTORY", true);
4157 poDS->m_bAddGDALHistory =
4158 CPLFetchBool(papszOptions, "ADD_GDAL_HISTORY", true);
4159 if( poDS->m_bAddGDALHistory )
4160 {
4161 poDS->m_osGDALHistory = CSLFetchNameValueDef(papszOptions,
4162 "GDAL_HISTORY", "");
4163 }
4164 const double dfNoData = (eType == GDT_Byte) ? NULL1:
4165 (eType == GDT_UInt16) ? NULLU2:
4166 (eType == GDT_Int16) ? NULL2:
4167 /*(eType == GDT_Float32) ?*/ NULL4;
4168
4169 for( int i = 0; i < nBands; i++ )
4170 {
4171 GDALRasterBand *poBand = nullptr;
4172
4173 if( poDS->m_poExternalDS != nullptr )
4174 {
4175 ISIS3WrapperRasterBand* poISISBand =
4176 new ISIS3WrapperRasterBand(
4177 poDS->m_poExternalDS->GetRasterBand( i+1 ) );
4178 poBand = poISISBand;
4179 }
4180 else if( bIsTiled )
4181 {
4182 ISISTiledBand* poISISBand =
4183 new ISISTiledBand( poDS, poDS->m_fpImage, i+1, eType,
4184 nBlockXSize, nBlockYSize,
4185 0, //nSkipBytes, to be hacked
4186 // afterwards for in-label imagery
4187 0, 0,
4188 CPL_IS_LSB );
4189
4190 poBand = poISISBand;
4191 }
4192 else
4193 {
4194 const int nPixelOffset = GDALGetDataTypeSizeBytes(eType);
4195 const int nLineOffset = nPixelOffset * nXSize;
4196 const vsi_l_offset nBandOffset =
4197 static_cast<vsi_l_offset>(nLineOffset) * nYSize;
4198 ISIS3RawRasterBand* poISISBand =
4199 new ISIS3RawRasterBand( poDS, i+1, poDS->m_fpImage,
4200 nBandOffset * i, // nImgOffset, to be
4201 //hacked afterwards for in-label imagery
4202 nPixelOffset, nLineOffset, eType,
4203 CPL_IS_LSB );
4204
4205 poBand = poISISBand;
4206 }
4207 poDS->SetBand( i+1, poBand );
4208 poBand->SetNoDataValue(dfNoData);
4209 }
4210
4211 return poDS;
4212 }
4213
4214 /************************************************************************/
4215 /* GetUnderlyingDataset() */
4216 /************************************************************************/
4217
GetUnderlyingDataset(GDALDataset * poSrcDS)4218 static GDALDataset* GetUnderlyingDataset( GDALDataset* poSrcDS )
4219 {
4220 if( poSrcDS->GetDriver() != nullptr &&
4221 poSrcDS->GetDriver() == GDALGetDriverByName("VRT") )
4222 {
4223 VRTDataset* poVRTDS = reinterpret_cast<VRTDataset* >(poSrcDS);
4224 poSrcDS = poVRTDS->GetSingleSimpleSource();
4225 }
4226
4227 return poSrcDS;
4228 }
4229
4230 /************************************************************************/
4231 /* CreateCopy() */
4232 /************************************************************************/
4233
CreateCopy(const char * pszFilename,GDALDataset * poSrcDS,int,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressData)4234 GDALDataset* ISIS3Dataset::CreateCopy( const char *pszFilename,
4235 GDALDataset *poSrcDS,
4236 int /*bStrict*/,
4237 char ** papszOptions,
4238 GDALProgressFunc pfnProgress,
4239 void * pProgressData )
4240 {
4241 const char* pszDataLocation = CSLFetchNameValueDef(papszOptions,
4242 "DATA_LOCATION",
4243 "LABEL");
4244 GDALDataset* poSrcUnderlyingDS = GetUnderlyingDataset(poSrcDS);
4245 if( poSrcUnderlyingDS == nullptr )
4246 poSrcUnderlyingDS = poSrcDS;
4247 if( EQUAL(pszDataLocation, "GEOTIFF") &&
4248 strcmp(poSrcUnderlyingDS->GetDescription(),
4249 CSLFetchNameValueDef(papszOptions, "EXTERNAL_FILENAME",
4250 CPLResetExtension(pszFilename, "tif"))
4251 ) == 0 )
4252 {
4253 CPLError(CE_Failure, CPLE_NotSupported,
4254 "Output file has same name as input file");
4255 return nullptr;
4256 }
4257 if( poSrcDS->GetRasterCount() == 0 )
4258 {
4259 CPLError(CE_Failure, CPLE_NotSupported, "Unsupported band count");
4260 return nullptr;
4261 }
4262
4263 const int nXSize = poSrcDS->GetRasterXSize();
4264 const int nYSize = poSrcDS->GetRasterYSize();
4265 const int nBands = poSrcDS->GetRasterCount();
4266 GDALDataType eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
4267 ISIS3Dataset *poDS = reinterpret_cast<ISIS3Dataset*>(
4268 Create( pszFilename, nXSize, nYSize, nBands, eType, papszOptions ));
4269 if( poDS == nullptr )
4270 return nullptr;
4271 poDS->m_osFromFilename = poSrcUnderlyingDS->GetDescription();
4272
4273 double adfGeoTransform[6] = { 0.0 };
4274 if( poSrcDS->GetGeoTransform( adfGeoTransform ) == CE_None
4275 && (adfGeoTransform[0] != 0.0
4276 || adfGeoTransform[1] != 1.0
4277 || adfGeoTransform[2] != 0.0
4278 || adfGeoTransform[3] != 0.0
4279 || adfGeoTransform[4] != 0.0
4280 || adfGeoTransform[5] != 1.0) )
4281 {
4282 poDS->SetGeoTransform( adfGeoTransform );
4283 }
4284
4285 auto poSrcSRS = poSrcDS->GetSpatialRef();
4286 if( poSrcSRS )
4287 {
4288 poDS->SetSpatialRef( poSrcSRS );
4289 }
4290
4291 for(int i=1;i<=nBands;i++)
4292 {
4293 const double dfOffset = poSrcDS->GetRasterBand(i)->GetOffset();
4294 if( dfOffset != 0.0 )
4295 poDS->GetRasterBand(i)->SetOffset(dfOffset);
4296
4297 const double dfScale = poSrcDS->GetRasterBand(i)->GetScale();
4298 if( dfScale != 1.0 )
4299 poDS->GetRasterBand(i)->SetScale(dfScale);
4300 }
4301
4302 // Do we need to remap nodata ?
4303 int bHasNoData = FALSE;
4304 poDS->m_dfSrcNoData =
4305 poSrcDS->GetRasterBand(1)->GetNoDataValue(&bHasNoData);
4306 poDS->m_bHasSrcNoData = CPL_TO_BOOL(bHasNoData);
4307
4308 if( poDS->m_bUseSrcLabel )
4309 {
4310 char** papszMD_ISIS3 = poSrcDS->GetMetadata("json:ISIS3");
4311 if( papszMD_ISIS3 != nullptr )
4312 {
4313 poDS->SetMetadata( papszMD_ISIS3, "json:ISIS3" );
4314 }
4315 }
4316
4317 // We don't need to initialize the imagery as we are going to copy it
4318 // completely
4319 poDS->m_bInitToNodata = false;
4320 CPLErr eErr = GDALDatasetCopyWholeRaster( poSrcDS, poDS,
4321 nullptr, pfnProgress, pProgressData );
4322 poDS->FlushCache();
4323 poDS->m_bHasSrcNoData = false;
4324 if( eErr != CE_None )
4325 {
4326 delete poDS;
4327 return nullptr;
4328 }
4329
4330 return poDS;
4331 }
4332
4333 /************************************************************************/
4334 /* GDALRegister_ISIS3() */
4335 /************************************************************************/
4336
GDALRegister_ISIS3()4337 void GDALRegister_ISIS3()
4338
4339 {
4340 if( GDALGetDriverByName( "ISIS3" ) != nullptr )
4341 return;
4342
4343 GDALDriver *poDriver = new GDALDriver();
4344
4345 poDriver->SetDescription( "ISIS3" );
4346 poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
4347 poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
4348 "USGS Astrogeology ISIS cube (Version 3)" );
4349 poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
4350 "drivers/raster/isis3.html" );
4351 poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
4352 poDriver->SetMetadataItem( GDAL_DMD_EXTENSIONS, "lbl cub" );
4353 poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
4354 "Byte UInt16 Int16 Float32" );
4355 poDriver->SetMetadataItem( GDAL_DMD_OPENOPTIONLIST, "<OpenOptionList/>");
4356 poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
4357 "<CreationOptionList>"
4358 " <Option name='DATA_LOCATION' type='string-select' "
4359 "description='Location of pixel data' default='LABEL'>"
4360 " <Value>LABEL</Value>"
4361 " <Value>EXTERNAL</Value>"
4362 " <Value>GEOTIFF</Value>"
4363 " </Option>"
4364 " <Option name='GEOTIFF_AS_REGULAR_EXTERNAL' type='boolean' "
4365 "description='Whether the GeoTIFF file, if uncompressed, should be "
4366 "registered as a regular raw file' default='YES'/>"
4367 " <Option name='GEOTIFF_OPTIONS' type='string' "
4368 "description='Comma separated list of KEY=VALUE tuples to forward "
4369 "to the GeoTIFF driver'/>"
4370 " <Option name='EXTERNAL_FILENAME' type='string' "
4371 "description='Override default external filename. "
4372 "Only for DATA_LOCATION=EXTERNAL or GEOTIFF'/>"
4373 " <Option name='TILED' type='boolean' "
4374 "description='Whether the pixel data should be tiled' default='NO'/>"
4375 " <Option name='BLOCKXSIZE' type='int' "
4376 "description='Tile width' default='256'/>"
4377 " <Option name='BLOCKYSIZE' type='int' "
4378 "description='Tile height' default='256'/>"
4379 " <Option name='COMMENT' type='string' "
4380 "description='Comment to add into the label'/>"
4381 " <Option name='LATITUDE_TYPE' type='string-select' "
4382 "description='Value of Mapping.LatitudeType' default='Planetocentric'>"
4383 " <Value>Planetocentric</Value>"
4384 " <Value>Planetographic</Value>"
4385 " </Option>"
4386 " <Option name='LONGITUDE_DIRECTION' type='string-select' "
4387 "description='Value of Mapping.LongitudeDirection' "
4388 "default='PositiveEast'>"
4389 " <Value>PositiveEast</Value>"
4390 " <Value>PositiveWest</Value>"
4391 " </Option>"
4392 " <Option name='TARGET_NAME' type='string' description='Value of "
4393 "Mapping.TargetName'/>"
4394 " <Option name='FORCE_360' type='boolean' "
4395 "description='Whether to force longitudes in [0,360] range' default='NO'/>"
4396 " <Option name='WRITE_BOUNDING_DEGREES' type='boolean' "
4397 "description='Whether to write Min/MaximumLong/Latitude values' "
4398 "default='YES'/>"
4399 " <Option name='BOUNDING_DEGREES' type='string' "
4400 "description='Manually set bounding box with the syntax "
4401 "min_long,min_lat,max_long,max_lat'/>"
4402 " <Option name='USE_SRC_LABEL' type='boolean' "
4403 "description='Whether to use source label in ISIS3 to ISIS3 conversions' "
4404 "default='YES'/>"
4405 " <Option name='USE_SRC_MAPPING' type='boolean' "
4406 "description='Whether to use Mapping group from source label in "
4407 "ISIS3 to ISIS3 conversions' "
4408 "default='NO'/>"
4409 " <Option name='USE_SRC_HISTORY' type='boolean' "
4410 "description='Whether to use content pointed by the History object in "
4411 "ISIS3 to ISIS3 conversions' "
4412 "default='YES'/>"
4413 " <Option name='ADD_GDAL_HISTORY' type='boolean' "
4414 "description='Whether to add GDAL specific history in the content pointed "
4415 "by the History object in "
4416 "ISIS3 to ISIS3 conversions' "
4417 "default='YES'/>"
4418 " <Option name='GDAL_HISTORY' type='string' "
4419 "description='Manually defined GDAL history. Must be formatted as ISIS3 "
4420 "PDL. If not specified, it is automatically composed.'/>"
4421 "</CreationOptionList>"
4422 );
4423
4424 poDriver->pfnOpen = ISIS3Dataset::Open;
4425 poDriver->pfnIdentify = ISIS3Dataset::Identify;
4426 poDriver->pfnCreate = ISIS3Dataset::Create;
4427 poDriver->pfnCreateCopy = ISIS3Dataset::CreateCopy;
4428
4429 GetGDALDriverManager()->RegisterDriver( poDriver );
4430 }
4431