1 /******************************************************************************
2  *
3  * Project:  ISCE Raster Reader
4  * Purpose:  Implementation of the ISCE raster reader
5  * Author:   Matthieu Volat (ISTerre), matthieu.volat@ujf-grenoble.fr
6  *
7  ******************************************************************************
8  * Copyright (c) 2015, Matthieu Volat <matthieu.volat@ujf-grenoble.fr>
9  *
10  * Permission is hereby granted, free of charge, to any person obtaining a
11  * copy of this software and associated documentation files (the "Software"),
12  * to deal in the Software without restriction, including without limitation
13  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14  * and/or sell copies of the Software, and to permit persons to whom the
15  * Software is furnished to do so, subject to the following conditions:
16  *
17  * The above copyright notice and this permission notice shall be included
18  * in all copies or substantial portions of the Software.
19  *
20  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26  * DEALINGS IN THE SOFTWARE.
27  ****************************************************************************/
28 
29 #include "gdal_frmts.h"
30 #include "ogr_spatialref.h"
31 #include "rawdataset.h"
32 
33 CPL_CVSID("$Id: iscedataset.cpp fa752ad6eabafaf630a704e1892a9d837d683cb3 2021-03-06 17:04:38 +0100 Even Rouault $")
34 
35 static const char * const apszISCE2GDALDatatypes[] = {
36     "BYTE:Byte",
37     "CHAR:Byte",
38     "SHORT:Int16",
39     "INT:Int32",
40     "LONG:Int64",
41     "FLOAT:Float32",
42     "DOUBLE:Float64",
43     "CBYTE:Unknown",
44     "CCHAR:Unknown",
45     "CSHORT:CInt16",
46     "CINT:CInt32",
47     "CLONG:CInt64",
48     "CFLOAT:CFloat32",
49     "CDOUBLE:CFloat64",
50     nullptr };
51 
52 static const char * const apszGDAL2ISCEDatatypes[] = {
53     "Byte:BYTE",
54     "Int16:SHORT",
55     "Int32:INT",
56     "Int64:LONG",
57     "Float32:FLOAT",
58     "Float64:DOUBLE",
59     "CInt16:CSHORT",
60     "CInt32:CINT",
61     "CInt64:CLONG",
62     "CFloat32:CFLOAT",
63     "CFloat64:CDOUBLE",
64     nullptr };
65 
66 enum Scheme { BIL = 0, BIP = 1, BSQ = 2 };
67 static const char * const apszSchemeNames[] = { "BIL", "BIP", "BSQ", nullptr };
68 
69 /************************************************************************/
70 /* ==================================================================== */
71 /*                              ISCEDataset                             */
72 /* ==================================================================== */
73 /************************************************************************/
74 
75 class ISCERasterBand;
76 
77 class ISCEDataset final: public RawDataset
78 {
79     friend class ISCERasterBand;
80 
81     VSILFILE    *fpImage;
82 
83     char        *pszXMLFilename;
84 
85     enum Scheme eScheme;
86 
87     CPL_DISALLOW_COPY_ASSIGN(ISCEDataset)
88 
89   public:
90     ISCEDataset();
91     ~ISCEDataset() override;
92 
93     void FlushCache() override;
94     char **GetFileList() override;
95 
96     static int          Identify( GDALOpenInfo *poOpenInfo );
97     static GDALDataset *Open( GDALOpenInfo *poOpenInfo );
98     static GDALDataset *Open( GDALOpenInfo *poOpenInfo, bool bFileSizeCheck );
99     static GDALDataset *Create( const char *pszFilename,
100                                 int nXSize, int nYSize, int nBands,
101                                 GDALDataType eType, char **papszOptions );
102 };
103 
104 /************************************************************************/
105 /* ==================================================================== */
106 /*                            ISCERasterBand                            */
107 /* ==================================================================== */
108 /************************************************************************/
109 
110 class ISCERasterBand final: public RawRasterBand
111 {
112         CPL_DISALLOW_COPY_ASSIGN(ISCERasterBand)
113 
114     public:
115                 ISCERasterBand( GDALDataset *poDS, int nBand, VSILFILE *fpRaw,
116                                   vsi_l_offset nImgOffset, int nPixelOffset,
117                                   int nLineOffset,
118                                   GDALDataType eDataType, int bNativeOrder );
119 };
120 
121 /************************************************************************/
122 /*                           getXMLFilename()                           */
123 /************************************************************************/
124 
getXMLFilename(GDALOpenInfo * poOpenInfo)125 static CPLString getXMLFilename( GDALOpenInfo *poOpenInfo )
126 {
127     CPLString osXMLFilename;
128 
129     if( poOpenInfo->fpL == nullptr )
130         return CPLString();
131 
132     char **papszSiblingFiles = poOpenInfo->GetSiblingFiles();
133     if ( papszSiblingFiles == nullptr )
134     {
135         osXMLFilename = CPLFormFilename( nullptr, poOpenInfo->pszFilename,
136                                          "xml" );
137         VSIStatBufL psXMLStatBuf;
138         if ( VSIStatL( osXMLFilename, &psXMLStatBuf ) != 0 )
139         {
140             osXMLFilename = "";
141         }
142     }
143     else
144     {
145         /* ------------------------------------------------------------ */
146         /*      We need to tear apart the filename to form a .xml       */
147         /*      filename.                                               */
148         /* ------------------------------------------------------------ */
149         const CPLString osPath = CPLGetPath( poOpenInfo->pszFilename );
150         const CPLString osName = CPLGetFilename( poOpenInfo->pszFilename );
151 
152         const int iFile =
153             CSLFindString( papszSiblingFiles,
154                            CPLFormFilename( nullptr, osName, "xml" ) );
155         if( iFile >= 0 )
156         {
157             osXMLFilename = CPLFormFilename( osPath,
158                                              papszSiblingFiles[iFile],
159                                              nullptr );
160         }
161     }
162 
163     return osXMLFilename;
164 }
165 
166 /************************************************************************/
167 /*                             ISCEDataset()                            */
168 /************************************************************************/
169 
ISCEDataset()170 ISCEDataset::ISCEDataset() :
171     fpImage(nullptr),
172     pszXMLFilename(nullptr),
173     eScheme(BIL)
174 {}
175 
176 /************************************************************************/
177 /*                            ~ISCEDataset()                          */
178 /************************************************************************/
179 
~ISCEDataset(void)180 ISCEDataset::~ISCEDataset( void )
181 {
182     ISCEDataset::FlushCache();
183     if ( fpImage != nullptr )
184     {
185         if( VSIFCloseL( fpImage ) != 0 )
186         {
187             CPLError(CE_Failure, CPLE_FileIO, "I/O error");
188         }
189     }
190     CPLFree( pszXMLFilename );
191 }
192 
193 /************************************************************************/
194 /*                            FlushCache()                              */
195 /************************************************************************/
196 
FlushCache(void)197 void ISCEDataset::FlushCache( void )
198 {
199     RawDataset::FlushCache();
200 
201     GDALRasterBand *band = (GetRasterCount() > 0) ? GetRasterBand(1) : nullptr;
202 
203     if ( eAccess == GA_ReadOnly || band == nullptr )
204         return;
205 
206 /* -------------------------------------------------------------------- */
207 /*      Recreate a XML doc with the dataset information.                */
208 /* -------------------------------------------------------------------- */
209     char sBuf[64] = { '\0' };
210     CPLXMLNode *psDocNode = CPLCreateXMLNode( nullptr, CXT_Element, "imageFile" );
211 
212     CPLXMLNode *psTmpNode
213         = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
214     CPLAddXMLAttributeAndValue( psTmpNode, "name", "WIDTH" );
215     CPLsnprintf(sBuf, sizeof(sBuf), "%d", nRasterXSize);
216     CPLCreateXMLElementAndValue( psTmpNode, "value", sBuf );
217 
218     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
219     CPLAddXMLAttributeAndValue( psTmpNode, "name", "LENGTH" );
220     CPLsnprintf(sBuf, sizeof(sBuf), "%d", nRasterYSize);
221     CPLCreateXMLElementAndValue( psTmpNode, "value", sBuf );
222 
223     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
224     CPLAddXMLAttributeAndValue( psTmpNode, "name", "NUMBER_BANDS" );
225     CPLsnprintf(sBuf, sizeof(sBuf), "%d", nBands);
226     CPLCreateXMLElementAndValue( psTmpNode, "value", sBuf );
227 
228     const char *sType = GDALGetDataTypeName( band->GetRasterDataType() );
229     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
230     CPLAddXMLAttributeAndValue( psTmpNode, "name", "DATA_TYPE" );
231     CPLCreateXMLElementAndValue(
232         psTmpNode, "value",
233         CSLFetchNameValue(
234             const_cast<char **>(apszGDAL2ISCEDatatypes),
235             sType ) );
236 
237     const char *pszScheme = apszSchemeNames[eScheme];
238     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
239     CPLAddXMLAttributeAndValue( psTmpNode, "name", "SCHEME" );
240     CPLCreateXMLElementAndValue( psTmpNode, "value", pszScheme );
241 
242     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
243     CPLAddXMLAttributeAndValue( psTmpNode, "name", "BYTE_ORDER" );
244 #ifdef CPL_LSB
245     CPLCreateXMLElementAndValue( psTmpNode, "value", "l" );
246 #else
247     CPLCreateXMLElementAndValue( psTmpNode, "value", "b" );
248 #endif
249 
250     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
251     CPLAddXMLAttributeAndValue( psTmpNode, "name", "ACCESS_MODE" );
252     CPLCreateXMLElementAndValue( psTmpNode, "value", "read" );
253 
254     const char *pszFilename = CPLGetBasename( pszXMLFilename );
255     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
256     CPLAddXMLAttributeAndValue( psTmpNode, "name", "FILE_NAME" );
257     CPLCreateXMLElementAndValue( psTmpNode, "value", pszFilename );
258 
259 /* -------------------------------------------------------------------- */
260 /*      Then, add the ISCE domain metadata.                             */
261 /* -------------------------------------------------------------------- */
262     char **papszISCEMetadata = GetMetadata( "ISCE" );
263     for (int i = 0; i < CSLCount( papszISCEMetadata ); i++)
264     {
265         /* Get the tokens from the metadata item */
266         char **papszTokens = CSLTokenizeString2( papszISCEMetadata[i],
267                                                  "=",
268                                                  CSLT_STRIPLEADSPACES
269                                                  | CSLT_STRIPENDSPACES);
270         if ( CSLCount( papszTokens ) != 2 )
271         {
272             CPLDebug( "ISCE",
273                       "Line of header file could not be split at = into two"
274                       " elements: %s",
275                       papszISCEMetadata[i] );
276             CSLDestroy( papszTokens );
277             continue;
278         }
279 
280         /* Don't write it out if it is one of the bits of metadata that is
281          * written out elsewhere in this routine */
282         if ( EQUAL( papszTokens[0], "WIDTH" )
283               || EQUAL( papszTokens[0], "LENGTH" )
284               || EQUAL( papszTokens[0], "NUMBER_BANDS" )
285               || EQUAL( papszTokens[0], "DATA_TYPE" )
286               || EQUAL( papszTokens[0], "SCHEME" )
287               || EQUAL( papszTokens[0], "BYTE_ORDER" ) )
288         {
289             CSLDestroy( papszTokens );
290             continue;
291         }
292 
293         psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
294         CPLAddXMLAttributeAndValue( psTmpNode, "name", papszTokens[0] );
295         CPLCreateXMLElementAndValue( psTmpNode, "value", papszTokens[1] );
296 
297         CSLDestroy( papszTokens );
298     }
299 
300 /* -------------------------------------------------------------------- */
301 /*      Create the "Coordinate" component elements, possibly with       */
302 /*      georeferencing.                                                 */
303 /* -------------------------------------------------------------------- */
304     CPLXMLNode *psCoordinate1Node, *psCoordinate2Node;
305     double adfGeoTransform[6];
306 
307     /* Coordinate 1 */
308     psCoordinate1Node = CPLCreateXMLNode( psDocNode,
309                                           CXT_Element,
310                                           "component" );
311     CPLAddXMLAttributeAndValue( psCoordinate1Node, "name", "Coordinate1" );
312     CPLCreateXMLElementAndValue( psCoordinate1Node,
313                                  "factorymodule",
314                                  "isceobj.Image" );
315     CPLCreateXMLElementAndValue( psCoordinate1Node,
316                                  "factoryname",
317                                  "createCoordinate" );
318     CPLCreateXMLElementAndValue( psCoordinate1Node,
319                                  "doc",
320                                  "First coordinate of a 2D image (width)." );
321     /* Property name */
322     psTmpNode = CPLCreateXMLNode( psCoordinate1Node,
323                                   CXT_Element,
324                                   "property" );
325     CPLAddXMLAttributeAndValue( psTmpNode, "name", "name" );
326     CPLCreateXMLElementAndValue( psTmpNode,
327                                  "value",
328                                  "ImageCoordinate_name" );
329     /* Property family */
330     psTmpNode = CPLCreateXMLNode( psCoordinate1Node,
331                                   CXT_Element,
332                                   "property" );
333     CPLAddXMLAttributeAndValue( psTmpNode, "name", "family" );
334     CPLCreateXMLElementAndValue( psTmpNode,
335                                  "value",
336                                  "ImageCoordinate" );
337     /* Property size */
338     CPLsnprintf(sBuf, sizeof(sBuf), "%d", nRasterXSize);
339     psTmpNode = CPLCreateXMLNode( psCoordinate1Node,
340                                   CXT_Element,
341                                   "property" );
342     CPLAddXMLAttributeAndValue( psTmpNode, "name", "size" );
343     CPLCreateXMLElementAndValue( psTmpNode,
344                                  "value",
345                                  sBuf );
346 
347     /* Coordinate 2 */
348     psCoordinate2Node = CPLCreateXMLNode( psDocNode,
349                                           CXT_Element,
350                                           "component" );
351     CPLAddXMLAttributeAndValue( psCoordinate2Node, "name", "Coordinate2" );
352     CPLCreateXMLElementAndValue( psCoordinate2Node,
353                                  "factorymodule",
354                                  "isceobj.Image" );
355     CPLCreateXMLElementAndValue( psCoordinate2Node,
356                                  "factoryname",
357                                  "createCoordinate" );
358     /* Property name */
359     psTmpNode = CPLCreateXMLNode( psCoordinate2Node,
360                                   CXT_Element,
361                                   "property" );
362     CPLAddXMLAttributeAndValue( psTmpNode, "name", "name" );
363     CPLCreateXMLElementAndValue( psTmpNode,
364                                  "value",
365                                  "ImageCoordinate_name" );
366     /* Property family */
367     psTmpNode = CPLCreateXMLNode( psCoordinate2Node,
368                                   CXT_Element,
369                                   "property" );
370     CPLAddXMLAttributeAndValue( psTmpNode, "name", "family" );
371     CPLCreateXMLElementAndValue( psTmpNode,
372                                  "value",
373                                  "ImageCoordinate" );
374     /* Property size */
375     CPLsnprintf(sBuf, sizeof(sBuf), "%d", nRasterYSize);
376     psTmpNode = CPLCreateXMLNode( psCoordinate2Node,
377                                   CXT_Element,
378                                   "property" );
379     CPLAddXMLAttributeAndValue( psTmpNode, "name", "size" );
380     CPLCreateXMLElementAndValue( psTmpNode,
381                                  "value",
382                                  sBuf );
383 
384     if ( GetGeoTransform( adfGeoTransform ) == CE_None )
385     {
386         if ( adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0 )
387         {
388             CPLError( CE_Warning, CPLE_AppDefined,
389                       "ISCE format do not support geotransform with "
390                           "rotation, discarding info.");
391         }
392         else {
393             CPLsnprintf( sBuf, sizeof(sBuf), "%g", adfGeoTransform[0] );
394             psTmpNode = CPLCreateXMLNode( psCoordinate1Node,
395                                           CXT_Element,
396                                           "property" );
397             CPLAddXMLAttributeAndValue( psTmpNode, "name", "startingValue" );
398             CPLCreateXMLElementAndValue( psTmpNode,
399                                          "value",
400                                          sBuf );
401 
402             CPLsnprintf( sBuf, sizeof(sBuf), "%g", adfGeoTransform[1] );
403             psTmpNode = CPLCreateXMLNode( psCoordinate1Node,
404                                           CXT_Element,
405                                           "property" );
406             CPLAddXMLAttributeAndValue( psTmpNode, "name", "delta" );
407             CPLCreateXMLElementAndValue( psTmpNode,
408                                          "value",
409                                          sBuf );
410 
411             CPLsnprintf( sBuf, sizeof(sBuf), "%g", adfGeoTransform[3] );
412             psTmpNode = CPLCreateXMLNode( psCoordinate2Node,
413                                           CXT_Element,
414                                           "property" );
415             CPLAddXMLAttributeAndValue( psTmpNode, "name", "startingValue" );
416             CPLCreateXMLElementAndValue( psTmpNode,
417                                          "value",
418                                          sBuf );
419 
420             CPLsnprintf( sBuf, sizeof(sBuf), "%g", adfGeoTransform[5] );
421             psTmpNode = CPLCreateXMLNode( psCoordinate2Node,
422                                           CXT_Element,
423                                           "property" );
424             CPLAddXMLAttributeAndValue( psTmpNode, "name", "delta" );
425             CPLCreateXMLElementAndValue( psTmpNode,
426                                          "value",
427                                          sBuf );
428         }
429     }
430 
431 /* -------------------------------------------------------------------- */
432 /*      Write the XML file.                                             */
433 /* -------------------------------------------------------------------- */
434     CPLSerializeXMLTreeToFile( psDocNode, pszXMLFilename );
435 
436 /* -------------------------------------------------------------------- */
437 /*      Free the XML Doc.                                               */
438 /* -------------------------------------------------------------------- */
439     CPLDestroyXMLNode( psDocNode );
440 }
441 
442 /************************************************************************/
443 /*                            GetFileList()                             */
444 /************************************************************************/
445 
GetFileList()446 char **ISCEDataset::GetFileList()
447 {
448     /* Main data file, etc. */
449     char **papszFileList = RawDataset::GetFileList();
450 
451     /* XML file. */
452     papszFileList = CSLAddString( papszFileList, pszXMLFilename );
453 
454     return papszFileList;
455 }
456 
457 /************************************************************************/
458 /*                             Identify()                               */
459 /************************************************************************/
460 
Identify(GDALOpenInfo * poOpenInfo)461 int ISCEDataset::Identify( GDALOpenInfo *poOpenInfo )
462 {
463 /* -------------------------------------------------------------------- */
464 /*      TODO: This function is unusable now:                            */
465 /*          * we can't just check for the presence of a XML file        */
466 /*          * we cannot parse it to check basic tree (Identify() is     */
467 /*            supposed to be faster than this                           */
468 /*          * we could read only a few bytes and strstr() for           */
469 /*            "imageData", but what if a file is padded with comments   */
470 /*            and/or whitespaces? it would still be legit, but the      */
471 /*            driver would fail...                                      */
472 /* -------------------------------------------------------------------- */
473 /* -------------------------------------------------------------------- */
474 /*      Check if there is a .xml file                                   */
475 /* -------------------------------------------------------------------- */
476     CPLString osXMLFilename = getXMLFilename( poOpenInfo );
477     if ( osXMLFilename.empty() )
478     {
479         return false;
480     }
481 
482     return true;
483 }
484 
485 /************************************************************************/
486 /*                                Open()                                */
487 /************************************************************************/
488 
Open(GDALOpenInfo * poOpenInfo)489 GDALDataset *ISCEDataset::Open( GDALOpenInfo *poOpenInfo )
490 {
491     return Open(poOpenInfo, true);
492 }
493 
Open(GDALOpenInfo * poOpenInfo,bool bFileSizeCheck)494 GDALDataset *ISCEDataset::Open( GDALOpenInfo *poOpenInfo, bool bFileSizeCheck )
495 {
496 /* -------------------------------------------------------------------- */
497 /*      Confirm that the header is compatible with a ISCE dataset.    */
498 /* -------------------------------------------------------------------- */
499     if ( !Identify(poOpenInfo) || poOpenInfo->fpL == nullptr )
500     {
501         return nullptr;
502     }
503 
504 /* -------------------------------------------------------------------- */
505 /*      Open and parse the .xml file                                    */
506 /* -------------------------------------------------------------------- */
507     const CPLString osXMLFilename = getXMLFilename( poOpenInfo );
508     CPLXMLNode *psNode = CPLParseXMLFile( osXMLFilename );
509     if ( psNode == nullptr || CPLGetXMLNode( psNode, "=imageFile" ) == nullptr )
510     {
511         CPLDestroyXMLNode( psNode );
512         return nullptr;
513     }
514     CPLXMLNode *psCur  = CPLGetXMLNode( psNode, "=imageFile" )->psChild;
515     char **papszXmlProps = nullptr;
516     while ( psCur != nullptr )
517     {
518         if ( EQUAL( psCur->pszValue, "property" ) )
519         {
520             /* Top-level property */
521             const char *pszName = CPLGetXMLValue( psCur, "name", nullptr );
522             const char *pszValue = CPLGetXMLValue( psCur, "value", nullptr );
523             if ( pszName != nullptr && pszValue != nullptr)
524             {
525                 papszXmlProps = CSLSetNameValue( papszXmlProps,
526                                                  pszName, pszValue );
527             }
528         }
529         else if ( EQUAL( psCur->pszValue, "component" ) )
530         {
531             /* "components" elements in ISCE store set of properties.   */
532             /* For now, they are avoided as I am not sure the full      */
533             /* scope of these. An exception is made for the ones named  */
534             /* Coordinate1 and Coordinate2, because they may have the   */
535             /* georeferencing information.                              */
536             const char *pszCurName = CPLGetXMLValue( psCur, "name", nullptr );
537             if ( pszCurName != nullptr
538                 && ( EQUAL( pszCurName, "Coordinate1" )
539                     || EQUAL( pszCurName, "Coordinate2" ) ) )
540             {
541                 /* We need two subproperties: startingValue and delta.  */
542                 /* To simplify parsing code, we will store them in      */
543                 /* papszXmlProps with the coordinate name prefixed to   */
544                 /* the property name.                                   */
545                 CPLXMLNode *psCur2 = psCur->psChild;
546                 while ( psCur2 != nullptr )
547                 {
548                     if ( ! EQUAL( psCur2->pszValue, "property" ) )
549                     {
550                         psCur2 = psCur2->psNext;
551                         continue; /* Skip non property elements */
552                     }
553 
554                     const char
555                        *pszCur2Name = CPLGetXMLValue( psCur2, "name", nullptr ),
556                        *pszCur2Value = CPLGetXMLValue( psCur2, "value", nullptr );
557 
558                     if ( pszCur2Name == nullptr || pszCur2Value == nullptr )
559                     {
560                         psCur2 = psCur2->psNext;
561                         continue; /* Skip malformatted elements */
562                     }
563 
564                     if ( EQUAL( pszCur2Name, "startingValue" )
565                         || EQUAL( pszCur2Name, "delta" ) )
566                     {
567                         char szPropName[32];
568                         snprintf(szPropName, sizeof(szPropName), "%s%s",
569                                  pszCurName, pszCur2Name);
570 
571                         papszXmlProps =
572                             CSLSetNameValue( papszXmlProps,
573                                              szPropName,
574                                              pszCur2Value );
575                     }
576                     psCur2 = psCur2->psNext;
577                 }
578             }
579         }
580         psCur = psCur->psNext;
581     }
582 
583     CPLDestroyXMLNode( psNode );
584 
585 /* -------------------------------------------------------------------- */
586 /*      Fetch required fields.                                          */
587 /* -------------------------------------------------------------------- */
588     if ( CSLFetchNameValue( papszXmlProps, "WIDTH" ) == nullptr
589         || CSLFetchNameValue( papszXmlProps, "LENGTH" ) == nullptr
590         || CSLFetchNameValue( papszXmlProps, "NUMBER_BANDS" ) == nullptr
591         || CSLFetchNameValue( papszXmlProps, "DATA_TYPE" ) == nullptr
592         || CSLFetchNameValue( papszXmlProps, "SCHEME" ) == nullptr )
593     {
594         CSLDestroy( papszXmlProps );
595         return nullptr;
596     }
597     const int nWidth = atoi( CSLFetchNameValue( papszXmlProps, "WIDTH" ) );
598     const int nHeight = atoi( CSLFetchNameValue( papszXmlProps, "LENGTH" ) );
599     const int nBands = atoi( CSLFetchNameValue( papszXmlProps, "NUMBER_BANDS" ) );
600 
601     if (!GDALCheckDatasetDimensions(nWidth, nHeight) ||
602         !GDALCheckBandCount(nBands, FALSE))
603     {
604         CSLDestroy( papszXmlProps );
605         return nullptr;
606     }
607 
608 /* -------------------------------------------------------------------- */
609 /*      Update byte order info if image specify something.              */
610 /* -------------------------------------------------------------------- */
611     bool bNativeOrder = true;
612 
613     const char *pszByteOrder = CSLFetchNameValue( papszXmlProps,
614                                                     "BYTE_ORDER" );
615     if ( pszByteOrder != nullptr )
616     {
617 #ifdef CPL_LSB
618         if ( EQUAL( pszByteOrder, "b" ) )
619 #else
620         if ( EQUAL( pszByteOrder, "l" ) )
621 #endif
622             bNativeOrder = false;
623     }
624 
625 /* -------------------------------------------------------------------- */
626 /*      Create a corresponding GDALDataset.                             */
627 /* -------------------------------------------------------------------- */
628     ISCEDataset *poDS = new ISCEDataset();
629     poDS->nRasterXSize = nWidth;
630     poDS->nRasterYSize = nHeight;
631     poDS->eAccess = poOpenInfo->eAccess;
632     poDS->pszXMLFilename = CPLStrdup( osXMLFilename.c_str() );
633     poDS->fpImage = poOpenInfo->fpL;
634     poOpenInfo->fpL = nullptr;
635 
636 /* -------------------------------------------------------------------- */
637 /*      Create band information objects.                                */
638 /* -------------------------------------------------------------------- */
639     const char *pszDataType =
640         CSLFetchNameValue( const_cast<char **>(apszISCE2GDALDatatypes),
641                            CSLFetchNameValue( papszXmlProps, "DATA_TYPE" ) );
642     if( pszDataType == nullptr )
643     {
644         delete poDS;
645         CSLDestroy( papszXmlProps );
646         return nullptr;
647     }
648     const GDALDataType eDataType = GDALGetDataTypeByName( pszDataType );
649     const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
650     if( nDTSize == 0 )
651     {
652         delete poDS;
653         CSLDestroy( papszXmlProps );
654         return nullptr;
655     }
656     const char *pszScheme = CSLFetchNameValue( papszXmlProps, "SCHEME" );
657     int nPixelOffset = 0;
658     int nLineOffset = 0;
659     vsi_l_offset nBandOffset = 0;
660     bool bIntOverflow = false;
661     if( EQUAL( pszScheme, "BIL" ) )
662     {
663         poDS->eScheme = BIL;
664         nPixelOffset = nDTSize;
665         if( nWidth > INT_MAX / (nPixelOffset * nBands) )
666             bIntOverflow = true;
667         else
668         {
669             nLineOffset = nPixelOffset * nWidth * nBands;
670             nBandOffset = nDTSize * static_cast<vsi_l_offset>(nWidth);
671         }
672     }
673     else if( EQUAL( pszScheme, "BIP" ) )
674     {
675         poDS->eScheme = BIP;
676         nPixelOffset = nDTSize * nBands;
677         if( nWidth > INT_MAX / nPixelOffset )
678             bIntOverflow = true;
679         else
680         {
681             nLineOffset = nPixelOffset * nWidth;
682             if( nBands > 1 && nLineOffset < INT_MAX / nBands )
683             {
684                 // GDAL 2.1.0 had a value of nLineOffset that was equal to the theoretical
685                 // nLineOffset multiplied by nBands...
686                 VSIFSeekL( poDS->fpImage, 0, SEEK_END );
687                 const GUIntBig nWrongFileSize = nDTSize *
688                     nWidth * (static_cast<GUIntBig>(nHeight - 1) *
689                     nBands * nBands + nBands);
690                 if( VSIFTellL( poDS->fpImage ) == nWrongFileSize )
691                 {
692                     CPLError(CE_Warning, CPLE_AppDefined,
693                             "This file has been incorrectly generated by an older "
694                             "GDAL version whose line offset computation was erroneous. "
695                             "Taking that into account, but the file should be re-encoded ideally");
696                     nLineOffset = nLineOffset * nBands;
697                 }
698             }
699             nBandOffset = nDTSize;
700         }
701     }
702     else if ( EQUAL( pszScheme, "BSQ" ) )
703     {
704         poDS->eScheme = BSQ;
705         nPixelOffset = nDTSize;
706         if( nWidth > INT_MAX / nPixelOffset )
707             bIntOverflow = true;
708         else
709         {
710             nLineOffset = nPixelOffset * nWidth;
711             nBandOffset = nLineOffset * static_cast<vsi_l_offset>(nHeight);
712         }
713     }
714     else
715     {
716         CPLError( CE_Failure, CPLE_OpenFailed,
717                   "Unknown scheme \"%s\" within ISCE raster.",
718                   pszScheme );
719         CSLDestroy( papszXmlProps );
720         delete poDS;
721         return nullptr;
722     }
723 
724     if (bIntOverflow)
725     {
726         delete poDS;
727         CPLError( CE_Failure, CPLE_AppDefined,
728                   "Int overflow occurred." );
729         CSLDestroy( papszXmlProps );
730         return nullptr;
731     }
732 
733     if( bFileSizeCheck &&
734         !RAWDatasetCheckMemoryUsage(
735                         poDS->nRasterXSize, poDS->nRasterYSize, nBands,
736                         nDTSize,
737                         nPixelOffset, nLineOffset, 0, nBandOffset,
738                         poDS->fpImage) )
739     {
740         delete poDS;
741         CSLDestroy( papszXmlProps );
742         return nullptr;
743     }
744 
745     poDS->nBands = nBands;
746     for (int b = 0; b < nBands; b++)
747     {
748         poDS->SetBand( b + 1,
749                        new ISCERasterBand( poDS, b + 1, poDS->fpImage,
750                                            nBandOffset * b,
751                                            nPixelOffset, nLineOffset,
752                                            eDataType, bNativeOrder ) );
753     }
754 
755 /* -------------------------------------------------------------------- */
756 /*      Interpret georeferencing, if present.                           */
757 /* -------------------------------------------------------------------- */
758     if ( CSLFetchNameValue( papszXmlProps, "Coordinate1startingValue" ) != nullptr
759          && CSLFetchNameValue( papszXmlProps, "Coordinate1delta" ) != nullptr
760          && CSLFetchNameValue( papszXmlProps, "Coordinate2startingValue" ) != nullptr
761          && CSLFetchNameValue( papszXmlProps, "Coordinate2delta" ) != nullptr )
762     {
763         double adfGeoTransform[6];
764         adfGeoTransform[0] = CPLAtof( CSLFetchNameValue( papszXmlProps,
765                                                          "Coordinate1startingValue" ) );
766         adfGeoTransform[1] = CPLAtof( CSLFetchNameValue( papszXmlProps,
767                                                          "Coordinate1delta" ) );
768         adfGeoTransform[2] = 0.0;
769         adfGeoTransform[3] = CPLAtof( CSLFetchNameValue( papszXmlProps,
770                                                          "Coordinate2startingValue" ) );
771         adfGeoTransform[4] = 0.0;
772         adfGeoTransform[5] = CPLAtof( CSLFetchNameValue( papszXmlProps,
773                                                                "Coordinate2delta" ) );
774         poDS->SetGeoTransform( adfGeoTransform );
775 
776         /* ISCE format seems not to have a projection field, but uses   */
777         /* WGS84.                                                       */
778         poDS->SetProjection( SRS_WKT_WGS84_LAT_LONG );
779     }
780 
781 /* -------------------------------------------------------------------- */
782 /*      Set all the other header metadata into the ISCE domain          */
783 /* -------------------------------------------------------------------- */
784     for( int i = 0; papszXmlProps != nullptr && papszXmlProps[i] != nullptr; i++ )
785     {
786         char **papszTokens =
787             CSLTokenizeString2( papszXmlProps[i],
788                                 "=",
789                                 CSLT_STRIPLEADSPACES
790                                 | CSLT_STRIPENDSPACES);
791         if ( CSLCount(papszTokens) < 2
792               || EQUAL( papszTokens[0], "WIDTH" )
793               || EQUAL( papszTokens[0], "LENGTH" )
794               || EQUAL( papszTokens[0], "NUMBER_BANDS" )
795               || EQUAL( papszTokens[0], "DATA_TYPE" )
796               || EQUAL( papszTokens[0], "SCHEME" )
797               || EQUAL( papszTokens[0], "BYTE_ORDER" )
798               || EQUAL( papszTokens[0], "Coordinate1startingValue" )
799               || EQUAL( papszTokens[0], "Coordinate1delta" )
800               || EQUAL( papszTokens[0], "Coordinate2startingValue" )
801               || EQUAL( papszTokens[0], "Coordinate2delta" ) )
802         {
803             CSLDestroy( papszTokens );
804             continue;
805         }
806         poDS->SetMetadataItem(papszTokens[0], papszTokens[1], "ISCE");
807         CSLDestroy( papszTokens );
808     }
809 
810 /* -------------------------------------------------------------------- */
811 /*      Free papszXmlProps                                              */
812 /* -------------------------------------------------------------------- */
813     CSLDestroy( papszXmlProps );
814 
815 /* -------------------------------------------------------------------- */
816 /*      Initialize any PAM information.                                 */
817 /* -------------------------------------------------------------------- */
818     poDS->SetDescription( poOpenInfo->pszFilename );
819     poDS->TryLoadXML();
820 
821 /* -------------------------------------------------------------------- */
822 /*      Check for overviews.                                            */
823 /* -------------------------------------------------------------------- */
824     poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
825 
826     return poDS;
827 }
828 
829 /************************************************************************/
830 /*                              Create()                                */
831 /************************************************************************/
832 
Create(const char * pszFilename,int nXSize,int nYSize,int nBands,GDALDataType eType,char ** papszOptions)833 GDALDataset *ISCEDataset::Create( const char *pszFilename,
834                                   int nXSize, int nYSize, int nBands,
835                                   GDALDataType eType,
836                                   char **papszOptions )
837 {
838     const char *sType = GDALGetDataTypeName( eType );
839     const char *pszScheme = CSLFetchNameValueDef( papszOptions,
840                                                 "SCHEME",
841                                                 "BIP" );
842 
843 /* -------------------------------------------------------------------- */
844 /*      Try to create the file.                                         */
845 /* -------------------------------------------------------------------- */
846     VSILFILE *fp = VSIFOpenL( pszFilename, "wb" );
847     if( fp == nullptr )
848     {
849         CPLError( CE_Failure, CPLE_OpenFailed,
850                   "Attempt to create file `%s' failed.",
851                   pszFilename );
852         return nullptr;
853     }
854 
855 /* -------------------------------------------------------------------- */
856 /*      Just write out a couple of bytes to establish the binary        */
857 /*      file, and then close it.                                        */
858 /* -------------------------------------------------------------------- */
859     CPL_IGNORE_RET_VAL(VSIFWriteL( "\0\0", 2, 1, fp ));
860     CPL_IGNORE_RET_VAL(VSIFCloseL( fp ));
861 
862 /* -------------------------------------------------------------------- */
863 /*      Create a minimal XML document.                                  */
864 /* -------------------------------------------------------------------- */
865     CPLXMLNode *psDocNode = CPLCreateXMLNode( nullptr, CXT_Element, "imageFile" );
866 
867     CPLXMLNode *psTmpNode
868         = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
869     CPLAddXMLAttributeAndValue( psTmpNode, "name", "WIDTH" );
870     char sBuf[64] = { '\0' };
871     CPLsnprintf(sBuf, sizeof(sBuf), "%d", nXSize);
872     CPLCreateXMLElementAndValue( psTmpNode, "value", sBuf );
873 
874     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
875     CPLAddXMLAttributeAndValue( psTmpNode, "name", "LENGTH" );
876     CPLsnprintf(sBuf, sizeof(sBuf), "%d", nYSize);
877     CPLCreateXMLElementAndValue( psTmpNode, "value", sBuf );
878 
879     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
880     CPLAddXMLAttributeAndValue( psTmpNode, "name", "NUMBER_BANDS" );
881     CPLsnprintf(sBuf, sizeof(sBuf), "%d", nBands);
882     CPLCreateXMLElementAndValue( psTmpNode, "value", sBuf );
883 
884     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
885     CPLAddXMLAttributeAndValue( psTmpNode, "name", "DATA_TYPE" );
886     CPLCreateXMLElementAndValue(
887         psTmpNode, "value",
888         CSLFetchNameValue(
889             const_cast<char **>(apszGDAL2ISCEDatatypes),
890             sType ));
891 
892     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
893     CPLAddXMLAttributeAndValue( psTmpNode, "name", "SCHEME" );
894     CPLCreateXMLElementAndValue( psTmpNode, "value", pszScheme );
895 
896     psTmpNode = CPLCreateXMLNode( psDocNode, CXT_Element, "property" );
897     CPLAddXMLAttributeAndValue( psTmpNode, "name", "BYTE_ORDER" );
898 #ifdef CPL_LSB
899     CPLCreateXMLElementAndValue( psTmpNode, "value", "l" );
900 #else
901     CPLCreateXMLElementAndValue( psTmpNode, "value", "b" );
902 #endif
903 
904 /* -------------------------------------------------------------------- */
905 /*      Write the XML file.                                             */
906 /* -------------------------------------------------------------------- */
907     const char  *pszXMLFilename = CPLFormFilename( nullptr, pszFilename, "xml" );
908     CPLSerializeXMLTreeToFile( psDocNode, pszXMLFilename );
909 
910 /* -------------------------------------------------------------------- */
911 /*      Free the XML Doc.                                               */
912 /* -------------------------------------------------------------------- */
913     CPLDestroyXMLNode( psDocNode );
914 
915     GDALOpenInfo oOpenInfo( pszFilename, GA_Update );
916     return Open(&oOpenInfo, false);
917 }
918 
919 /************************************************************************/
920 /*                          ISCERasterBand()                            */
921 /************************************************************************/
922 
ISCERasterBand(GDALDataset * poDSIn,int nBandIn,VSILFILE * fpRawIn,vsi_l_offset nImgOffsetIn,int nPixelOffsetIn,int nLineOffsetIn,GDALDataType eDataTypeIn,int bNativeOrderIn)923 ISCERasterBand::ISCERasterBand( GDALDataset *poDSIn, int nBandIn, VSILFILE *fpRawIn,
924                                 vsi_l_offset nImgOffsetIn, int nPixelOffsetIn,
925                                 int nLineOffsetIn,
926                                 GDALDataType eDataTypeIn, int bNativeOrderIn ) :
927     RawRasterBand( poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
928                    nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
929                    RawRasterBand::OwnFP::NO )
930 {}
931 
932 /************************************************************************/
933 /*                         GDALRegister_ISCE()                          */
934 /************************************************************************/
935 
GDALRegister_ISCE()936 void GDALRegister_ISCE()
937 {
938     if( GDALGetDriverByName( "ISCE" ) != nullptr )
939         return;
940 
941     GDALDriver *poDriver = new GDALDriver();
942 
943     poDriver->SetDescription( "ISCE" );
944     poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "ISCE raster" );
945     poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "drivers/raster/isce.html" );
946     poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
947                                "Byte Int16 Int32 Int64 Float32"
948                                " Float64 CInt16 CInt64 CFloat32 "
949                                " CFloat64" );
950     poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
951                                "<CreationOptionList>"
952                                "   <Option name='SCHEME' type='string-select'>"
953                                "       <Value>BIP</Value>"
954                                "       <Value>BIL</Value>"
955                                "       <Value>BSQ</Value>"
956                                "   </Option>"
957                                "</CreationOptionList>" );
958     poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
959     poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
960 
961     poDriver->pfnOpen = ISCEDataset::Open;
962     poDriver->pfnCreate = ISCEDataset::Create;
963 
964     GetGDALDriverManager()->RegisterDriver( poDriver );
965 }
966