Coverage Report

Created: 2025-07-23 09:13

/src/gdal/frmts/pds/isis2dataset.cpp
Line
Count
Source (jump to first uncovered line)
1
/******************************************************************************
2
 *
3
 * Project:  ISIS Version 2 Driver
4
 * Purpose:  Implementation of ISIS2Dataset
5
 * Author:   Trent Hare (thare at usgs.gov),
6
 *           Robert Soricone (rsoricone at usgs.gov)
7
 *           Ludovic Mercier (ludovic.mercier at gmail.com)
8
 *           Frank Warmerdam (warmerdam at pobox.com)
9
 *
10
 * NOTE: Original code authored by Trent and Robert and placed in the public
11
 * domain as per US government policy.  I have (within my rights) appropriated
12
 * it and placed it under the following license.  This is not intended to
13
 * diminish Trent and Roberts contribution.
14
 ******************************************************************************
15
 * Copyright (c) 2006, Frank Warmerdam <warmerdam at pobox.com>
16
 * Copyright (c) 2008-2011, Even Rouault <even dot rouault at spatialys.com>
17
 *
18
 * SPDX-License-Identifier: MIT
19
 ****************************************************************************/
20
21
constexpr int NULL1 = 0;
22
constexpr int NULL2 = -32768;
23
constexpr double NULL3 = -3.4028226550889044521e+38;
24
25
#include "cpl_string.h"
26
#include "gdal_frmts.h"
27
#include "nasakeywordhandler.h"
28
#include "ogr_spatialref.h"
29
#include "rawdataset.h"
30
#include "pdsdrivercore.h"
31
32
/************************************************************************/
33
/* ==================================================================== */
34
/*                      ISISDataset     version2                        */
35
/* ==================================================================== */
36
/************************************************************************/
37
38
class ISIS2Dataset final : public RawDataset
39
{
40
    VSILFILE *fpImage{};  // image data file.
41
    CPLString osExternalCube{};
42
43
    NASAKeywordHandler oKeywords{};
44
45
    bool bGotTransform{};
46
    GDALGeoTransform m_gt{};
47
48
    OGRSpatialReference m_oSRS{};
49
50
    int parse_label(const char *file, char *keyword, char *value);
51
    int strstrip(char instr[], char outstr[], int position);
52
53
    CPLString oTempResult{};
54
55
    static void CleanString(CPLString &osInput);
56
57
    const char *GetKeyword(const char *pszPath, const char *pszDefault = "");
58
    const char *GetKeywordSub(const char *pszPath, int iSubscript,
59
                              const char *pszDefault = "");
60
61
    CPLErr Close() override;
62
63
    CPL_DISALLOW_COPY_ASSIGN(ISIS2Dataset)
64
65
  public:
66
    ISIS2Dataset();
67
    virtual ~ISIS2Dataset();
68
69
    virtual CPLErr GetGeoTransform(GDALGeoTransform &gt) const override;
70
    const OGRSpatialReference *GetSpatialRef() const override;
71
72
    virtual char **GetFileList() override;
73
74
    static GDALDataset *Open(GDALOpenInfo *);
75
};
76
77
/************************************************************************/
78
/*                            ISIS2Dataset()                            */
79
/************************************************************************/
80
81
ISIS2Dataset::ISIS2Dataset()
82
3.69k
{
83
3.69k
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
84
3.69k
}
85
86
/************************************************************************/
87
/*                            ~ISIS2Dataset()                            */
88
/************************************************************************/
89
90
ISIS2Dataset::~ISIS2Dataset()
91
92
3.69k
{
93
3.69k
    ISIS2Dataset::Close();
94
3.69k
}
95
96
/************************************************************************/
97
/*                              Close()                                 */
98
/************************************************************************/
99
100
CPLErr ISIS2Dataset::Close()
101
3.69k
{
102
3.69k
    CPLErr eErr = CE_None;
103
3.69k
    if (nOpenFlags != OPEN_FLAGS_CLOSED)
104
3.69k
    {
105
3.69k
        if (ISIS2Dataset::FlushCache(true) != CE_None)
106
0
            eErr = CE_Failure;
107
3.69k
        if (fpImage != nullptr)
108
4
        {
109
4
            if (VSIFCloseL(fpImage) != 0)
110
0
                eErr = CE_Failure;
111
4
        }
112
3.69k
        if (GDALPamDataset::Close() != CE_None)
113
0
            eErr = CE_Failure;
114
3.69k
    }
115
3.69k
    return eErr;
116
3.69k
}
117
118
/************************************************************************/
119
/*                            GetFileList()                             */
120
/************************************************************************/
121
122
char **ISIS2Dataset::GetFileList()
123
124
8
{
125
8
    char **papszFileList = GDALPamDataset::GetFileList();
126
127
8
    if (!osExternalCube.empty())
128
0
        papszFileList = CSLAddString(papszFileList, osExternalCube);
129
130
8
    return papszFileList;
131
8
}
132
133
/************************************************************************/
134
/*                         GetSpatialRef()                              */
135
/************************************************************************/
136
137
const OGRSpatialReference *ISIS2Dataset::GetSpatialRef() const
138
4
{
139
4
    if (!m_oSRS.IsEmpty())
140
2
        return &m_oSRS;
141
2
    return GDALPamDataset::GetSpatialRef();
142
4
}
143
144
/************************************************************************/
145
/*                          GetGeoTransform()                           */
146
/************************************************************************/
147
148
CPLErr ISIS2Dataset::GetGeoTransform(GDALGeoTransform &gt) const
149
150
4
{
151
4
    if (bGotTransform)
152
3
    {
153
3
        gt = m_gt;
154
3
        return CE_None;
155
3
    }
156
157
1
    return GDALPamDataset::GetGeoTransform(gt);
158
4
}
159
160
/************************************************************************/
161
/*                                Open()                                */
162
/************************************************************************/
163
164
GDALDataset *ISIS2Dataset::Open(GDALOpenInfo *poOpenInfo)
165
3.69k
{
166
    /* -------------------------------------------------------------------- */
167
    /*      Does this look like a CUBE or an IMAGE Primary Data Object?     */
168
    /* -------------------------------------------------------------------- */
169
3.69k
    if (!ISIS2DriverIdentify(poOpenInfo) || poOpenInfo->fpL == nullptr)
170
0
        return nullptr;
171
172
3.69k
    VSILFILE *fpQube = poOpenInfo->fpL;
173
3.69k
    poOpenInfo->fpL = nullptr;
174
175
3.69k
    auto poDS = std::make_unique<ISIS2Dataset>();
176
177
3.69k
    if (!poDS->oKeywords.Ingest(fpQube, 0))
178
3.27k
    {
179
3.27k
        VSIFCloseL(fpQube);
180
3.27k
        return nullptr;
181
3.27k
    }
182
183
421
    VSIFCloseL(fpQube);
184
185
    /* -------------------------------------------------------------------- */
186
    /*      We assume the user is pointing to the label (i.e. .lab) file.   */
187
    /* -------------------------------------------------------------------- */
188
    // QUBE can be inline or detached and point to an image name
189
    // ^QUBE = 76
190
    // ^QUBE = ("ui31s015.img",6441<BYTES>) - has another label on the image
191
    // ^QUBE = "ui31s015.img" - which implies no label or skip value
192
193
421
    const char *pszQube = poDS->GetKeyword("^QUBE");
194
421
    int nQube = 0;
195
421
    int bByteLocation = FALSE;
196
421
    CPLString osTargetFile = poOpenInfo->pszFilename;
197
198
421
    if (pszQube[0] == '"')
199
31
    {
200
31
        const CPLString osTPath = CPLGetPathSafe(poOpenInfo->pszFilename);
201
31
        CPLString osFilename = pszQube;
202
31
        poDS->CleanString(osFilename);
203
31
        osTargetFile = CPLFormCIFilenameSafe(osTPath, osFilename, nullptr);
204
31
        poDS->osExternalCube = osTargetFile;
205
31
    }
206
390
    else if (pszQube[0] == '(')
207
0
    {
208
0
        const CPLString osTPath = CPLGetPathSafe(poOpenInfo->pszFilename);
209
0
        CPLString osFilename = poDS->GetKeywordSub("^QUBE", 1, "");
210
0
        poDS->CleanString(osFilename);
211
0
        osTargetFile = CPLFormCIFilenameSafe(osTPath, osFilename, nullptr);
212
0
        poDS->osExternalCube = osTargetFile;
213
214
0
        nQube = atoi(poDS->GetKeywordSub("^QUBE", 2, "1"));
215
0
        if (strstr(poDS->GetKeywordSub("^QUBE", 2, "1"), "<BYTES>") != nullptr)
216
0
            bByteLocation = true;
217
0
    }
218
390
    else
219
390
    {
220
390
        nQube = atoi(pszQube);
221
390
        if (strstr(pszQube, "<BYTES>") != nullptr)
222
0
            bByteLocation = true;
223
390
    }
224
225
    /* -------------------------------------------------------------------- */
226
    /*      Check if file an ISIS2 header file?  Read a few lines of text   */
227
    /*      searching for something starting with nrows or ncols.           */
228
    /* -------------------------------------------------------------------- */
229
230
    /* -------------------------------------------------------------------- */
231
    /*      Checks to see if this is valid ISIS2 cube                       */
232
    /*      SUFFIX_ITEM tag in .cub file should be (0,0,0); no side-planes  */
233
    /* -------------------------------------------------------------------- */
234
421
    const int s_ix = atoi(poDS->GetKeywordSub("QUBE.SUFFIX_ITEMS", 1));
235
421
    const int s_iy = atoi(poDS->GetKeywordSub("QUBE.SUFFIX_ITEMS", 2));
236
421
    const int s_iz = atoi(poDS->GetKeywordSub("QUBE.SUFFIX_ITEMS", 3));
237
238
421
    if (s_ix != 0 || s_iy != 0 || s_iz != 0)
239
0
    {
240
0
        CPLError(CE_Failure, CPLE_OpenFailed,
241
0
                 "*** ISIS 2 cube file has invalid SUFFIX_ITEMS parameters:\n"
242
0
                 "*** gdal isis2 driver requires (0, 0, 0), thus no sideplanes "
243
0
                 "or backplanes\n"
244
0
                 "found: (%i, %i, %i)\n\n",
245
0
                 s_ix, s_iy, s_iz);
246
0
        return nullptr;
247
0
    }
248
249
    /**************** end SUFFIX_ITEM check ***********************/
250
251
    /***********   Grab layout type (BSQ, BIP, BIL) ************/
252
    //  AXIS_NAME = (SAMPLE,LINE,BAND)
253
    /***********************************************************/
254
255
421
    char szLayout[10] = "BSQ";  // default to band seq.
256
421
    const char *value = poDS->GetKeyword("QUBE.AXIS_NAME", "");
257
421
    if (EQUAL(value, "(SAMPLE,LINE,BAND)"))
258
5
        strcpy(szLayout, "BSQ");
259
416
    else if (EQUAL(value, "(BAND,LINE,SAMPLE)"))
260
0
        strcpy(szLayout, "BIP");
261
416
    else if (EQUAL(value, "(SAMPLE,BAND,LINE)") || EQUAL(value, ""))
262
416
        strcpy(szLayout, "BSQ");
263
0
    else
264
0
    {
265
0
        CPLError(CE_Failure, CPLE_OpenFailed,
266
0
                 "%s layout not supported. Abort\n\n", value);
267
0
        return nullptr;
268
0
    }
269
270
    /***********   Grab samples lines band ************/
271
421
    const int nCols = atoi(poDS->GetKeywordSub("QUBE.CORE_ITEMS", 1));
272
421
    const int nRows = atoi(poDS->GetKeywordSub("QUBE.CORE_ITEMS", 2));
273
421
    const int nBands = atoi(poDS->GetKeywordSub("QUBE.CORE_ITEMS", 3));
274
275
    /***********   Grab Qube record bytes  **********/
276
421
    const int record_bytes = atoi(poDS->GetKeyword("RECORD_BYTES"));
277
421
    if (record_bytes < 0)
278
0
    {
279
0
        return nullptr;
280
0
    }
281
282
421
    GUIntBig nSkipBytes = 0;
283
421
    if (nQube > 0 && bByteLocation)
284
0
        nSkipBytes = (nQube - 1);
285
421
    else if (nQube > 0)
286
4
        nSkipBytes = static_cast<GUIntBig>(nQube - 1) * record_bytes;
287
417
    else
288
417
        nSkipBytes = 0;
289
290
    /***********   Grab samples lines band ************/
291
421
    char chByteOrder = 'M';  // default to MSB
292
421
    CPLString osCoreItemType = poDS->GetKeyword("QUBE.CORE_ITEM_TYPE");
293
421
    if ((EQUAL(osCoreItemType, "PC_INTEGER")) ||
294
421
        (EQUAL(osCoreItemType, "PC_UNSIGNED_INTEGER")) ||
295
421
        (EQUAL(osCoreItemType, "PC_REAL")))
296
0
    {
297
0
        chByteOrder = 'I';
298
0
    }
299
300
    /********   Grab format type - isis2 only supports 8,16,32 *******/
301
421
    GDALDataType eDataType = GDT_Byte;
302
421
    bool bNoDataSet = false;
303
421
    double dfNoData = 0.0;
304
305
421
    int itype = atoi(poDS->GetKeyword("QUBE.CORE_ITEM_BYTES", ""));
306
421
    switch (itype)
307
421
    {
308
0
        case 1:
309
0
            eDataType = GDT_Byte;
310
0
            dfNoData = NULL1;
311
0
            bNoDataSet = true;
312
0
            break;
313
0
        case 2:
314
0
            if (strstr(osCoreItemType, "UNSIGNED") != nullptr)
315
0
            {
316
0
                dfNoData = 0;
317
0
                eDataType = GDT_UInt16;
318
0
            }
319
0
            else
320
0
            {
321
0
                dfNoData = NULL2;
322
0
                eDataType = GDT_Int16;
323
0
            }
324
0
            bNoDataSet = true;
325
0
            break;
326
5
        case 4:
327
5
            eDataType = GDT_Float32;
328
5
            dfNoData = NULL3;
329
5
            bNoDataSet = true;
330
5
            break;
331
0
        case 8:
332
0
            eDataType = GDT_Float64;
333
0
            dfNoData = NULL3;
334
0
            bNoDataSet = true;
335
0
            break;
336
416
        default:
337
416
            CPLError(CE_Failure, CPLE_AppDefined,
338
416
                     "Itype of %d is not supported in ISIS 2.", itype);
339
416
            return nullptr;
340
421
    }
341
342
    /***********   Grab Cellsize ************/
343
5
    double dfXDim = 1.0;
344
5
    double dfYDim = 1.0;
345
346
5
    value = poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.MAP_SCALE");
347
5
    if (strlen(value) > 0)
348
3
    {
349
        // Convert km to m
350
3
        dfXDim = static_cast<float>(CPLAtof(value) * 1000.0);
351
3
        dfYDim = static_cast<float>(CPLAtof(value) * 1000.0 * -1);
352
3
    }
353
354
    /***********   Grab LINE_PROJECTION_OFFSET ************/
355
5
    double dfULYMap = 0.5;
356
357
5
    value =
358
5
        poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.LINE_PROJECTION_OFFSET");
359
5
    if (strlen(value) > 0)
360
4
    {
361
4
        const double yulcenter = static_cast<float>(CPLAtof(value)) * dfYDim;
362
4
        dfULYMap = yulcenter - (dfYDim / 2);
363
4
    }
364
365
    /***********   Grab SAMPLE_PROJECTION_OFFSET ************/
366
5
    double dfULXMap = 0.5;
367
368
5
    value =
369
5
        poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.SAMPLE_PROJECTION_OFFSET");
370
5
    if (strlen(value) > 0)
371
3
    {
372
3
        const double xulcenter = static_cast<float>(CPLAtof(value)) * dfXDim;
373
3
        dfULXMap = xulcenter - (dfXDim / 2);
374
3
    }
375
376
    /***********  Grab TARGET_NAME  ************/
377
    /**** This is the planets name i.e. MARS ***/
378
5
    CPLString target_name = poDS->GetKeyword("QUBE.TARGET_NAME");
379
380
    /***********   Grab MAP_PROJECTION_TYPE ************/
381
5
    CPLString map_proj_name =
382
5
        poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.MAP_PROJECTION_TYPE");
383
5
    poDS->CleanString(map_proj_name);
384
385
    /***********   Grab SEMI-MAJOR ************/
386
5
    const double semi_major =
387
5
        CPLAtof(poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.A_AXIS_RADIUS")) *
388
5
        1000.0;
389
390
    /***********   Grab semi-minor ************/
391
5
    const double semi_minor =
392
5
        CPLAtof(poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.C_AXIS_RADIUS")) *
393
5
        1000.0;
394
395
    /***********   Grab CENTER_LAT ************/
396
5
    const double center_lat =
397
5
        CPLAtof(poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.CENTER_LATITUDE"));
398
399
    /***********   Grab CENTER_LON ************/
400
5
    const double center_lon =
401
5
        CPLAtof(poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.CENTER_LONGITUDE"));
402
403
    /***********   Grab 1st std parallel ************/
404
5
    const double first_std_parallel = CPLAtof(
405
5
        poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.FIRST_STANDARD_PARALLEL"));
406
407
    /***********   Grab 2nd std parallel ************/
408
5
    const double second_std_parallel = CPLAtof(
409
5
        poDS->GetKeyword("QUBE.IMAGE_MAP_PROJECTION.SECOND_STANDARD_PARALLEL"));
410
411
    /*** grab  PROJECTION_LATITUDE_TYPE = "PLANETOCENTRIC" ****/
412
    // Need to further study how ocentric/ographic will effect the gdal library.
413
    // So far we will use this fact to define a sphere or ellipse for some
414
    // projections Frank - may need to talk this over
415
5
    bool bIsGeographic = true;
416
5
    value =
417
5
        poDS->GetKeyword("CUBE.IMAGE_MAP_PROJECTION.PROJECTION_LATITUDE_TYPE");
418
5
    if (EQUAL(value, "\"PLANETOCENTRIC\""))
419
0
        bIsGeographic = false;
420
421
5
    CPLDebug("ISIS2", "using projection %s", map_proj_name.c_str());
422
423
5
    OGRSpatialReference oSRS;
424
5
    oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
425
5
    bool bProjectionSet = true;
426
427
    // Set oSRS projection and parameters
428
5
    if ((EQUAL(map_proj_name, "EQUIRECTANGULAR_CYLINDRICAL")) ||
429
5
        (EQUAL(map_proj_name, "EQUIRECTANGULAR")) ||
430
5
        (EQUAL(map_proj_name, "SIMPLE_CYLINDRICAL")))
431
3
    {
432
3
        oSRS.OGRSpatialReference::SetEquirectangular2(0.0, center_lon,
433
3
                                                      center_lat, 0, 0);
434
3
    }
435
2
    else if (EQUAL(map_proj_name, "ORTHOGRAPHIC"))
436
0
    {
437
0
        oSRS.OGRSpatialReference::SetOrthographic(center_lat, center_lon, 0, 0);
438
0
    }
439
2
    else if ((EQUAL(map_proj_name, "SINUSOIDAL")) ||
440
2
             (EQUAL(map_proj_name, "SINUSOIDAL_EQUAL-AREA")))
441
0
    {
442
0
        oSRS.OGRSpatialReference::SetSinusoidal(center_lon, 0, 0);
443
0
    }
444
2
    else if (EQUAL(map_proj_name, "MERCATOR"))
445
0
    {
446
0
        oSRS.OGRSpatialReference::SetMercator(center_lat, center_lon, 1, 0, 0);
447
0
    }
448
2
    else if (EQUAL(map_proj_name, "POLAR_STEREOGRAPHIC"))
449
0
    {
450
0
        oSRS.OGRSpatialReference::SetPS(center_lat, center_lon, 1, 0, 0);
451
0
    }
452
2
    else if (EQUAL(map_proj_name, "TRANSVERSE_MERCATOR"))
453
0
    {
454
0
        oSRS.OGRSpatialReference::SetTM(center_lat, center_lon, 1, 0, 0);
455
0
    }
456
2
    else if (EQUAL(map_proj_name, "LAMBERT_CONFORMAL_CONIC"))
457
0
    {
458
0
        oSRS.OGRSpatialReference::SetLCC(first_std_parallel,
459
0
                                         second_std_parallel, center_lat,
460
0
                                         center_lon, 0, 0);
461
0
    }
462
2
    else if (EQUAL(map_proj_name, ""))
463
1
    {
464
        /* no projection */
465
1
        bProjectionSet = false;
466
1
    }
467
1
    else
468
1
    {
469
1
        CPLDebug("ISIS2",
470
1
                 "Dataset projection %s is not supported. Continuing...",
471
1
                 map_proj_name.c_str());
472
1
        bProjectionSet = false;
473
1
    }
474
475
5
    if (bProjectionSet)
476
3
    {
477
        // Create projection name, i.e. MERCATOR MARS and set as ProjCS keyword
478
3
        const CPLString proj_target_name = map_proj_name + " " + target_name;
479
3
        oSRS.SetProjCS(proj_target_name);  // set ProjCS keyword
480
481
        // The geographic/geocentric name will be the same basic name as the
482
        // body name 'GCS' = Geographic/Geocentric Coordinate System
483
3
        const CPLString geog_name = "GCS_" + target_name;
484
485
        // The datum and sphere names will be the same basic name aas the planet
486
3
        const CPLString datum_name = "D_" + target_name;
487
488
3
        CPLString sphere_name = std::move(target_name);
489
490
        // calculate inverse flattening from major and minor axis: 1/f = a/(a-b)
491
3
        double iflattening = 0.0;
492
3
        if ((semi_major - semi_minor) < 0.0000001)
493
3
            iflattening = 0;
494
0
        else
495
0
            iflattening = semi_major / (semi_major - semi_minor);
496
497
        // Set the body size but take into consideration which proj is being
498
        // used to help w/ proj4 compatibility The use of a Sphere, polar radius
499
        // or ellipse here is based on how ISIS does it internally
500
3
        if (((EQUAL(map_proj_name, "STEREOGRAPHIC") &&
501
3
              (fabs(center_lat) == 90))) ||
502
3
            (EQUAL(map_proj_name, "POLAR_STEREOGRAPHIC")))
503
0
        {
504
0
            if (bIsGeographic)
505
0
            {
506
                // Geograpraphic, so set an ellipse
507
0
                oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_major,
508
0
                               iflattening, "Reference_Meridian", 0.0);
509
0
            }
510
0
            else
511
0
            {
512
                // Geocentric, so force a sphere using the semi-minor axis. I
513
                // hope...
514
0
                sphere_name += "_polarRadius";
515
0
                oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_minor,
516
0
                               0.0, "Reference_Meridian", 0.0);
517
0
            }
518
0
        }
519
3
        else if ((EQUAL(map_proj_name, "SIMPLE_CYLINDRICAL")) ||
520
3
                 (EQUAL(map_proj_name, "ORTHOGRAPHIC")) ||
521
3
                 (EQUAL(map_proj_name, "STEREOGRAPHIC")) ||
522
3
                 (EQUAL(map_proj_name, "SINUSOIDAL_EQUAL-AREA")) ||
523
3
                 (EQUAL(map_proj_name, "SINUSOIDAL")))
524
3
        {
525
            // ISIS uses the spherical equation for these projections so force
526
            // a sphere.
527
3
            oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_major, 0.0,
528
3
                           "Reference_Meridian", 0.0);
529
3
        }
530
0
        else if ((EQUAL(map_proj_name, "EQUIRECTANGULAR_CYLINDRICAL")) ||
531
0
                 (EQUAL(map_proj_name, "EQUIRECTANGULAR")))
532
0
        {
533
            // Calculate localRadius using ISIS3 simple elliptical method
534
            //   not the more standard Radius of Curvature method
535
            // PI = 4 * atan(1);
536
0
            if (center_lon == 0)
537
0
            {  // No need to calculate local radius
538
0
                oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_major,
539
0
                               0.0, "Reference_Meridian", 0.0);
540
0
            }
541
0
            else
542
0
            {
543
0
                const double radLat = center_lat * M_PI / 180;  // in radians
544
0
                const double localRadius =
545
0
                    semi_major * semi_minor /
546
0
                    sqrt(pow(semi_minor * cos(radLat), 2) +
547
0
                         pow(semi_major * sin(radLat), 2));
548
0
                sphere_name += "_localRadius";
549
0
                oSRS.SetGeogCS(geog_name, datum_name, sphere_name, localRadius,
550
0
                               0.0, "Reference_Meridian", 0.0);
551
0
                CPLDebug("ISIS2", "local radius: %f", localRadius);
552
0
            }
553
0
        }
554
0
        else
555
0
        {
556
            // All other projections: Mercator, Transverse Mercator, Lambert
557
            // Conformal, etc. Geographic, so set an ellipse
558
0
            if (bIsGeographic)
559
0
            {
560
0
                oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_major,
561
0
                               iflattening, "Reference_Meridian", 0.0);
562
0
            }
563
0
            else
564
0
            {
565
                // Geocentric, so force a sphere. I hope...
566
0
                oSRS.SetGeogCS(geog_name, datum_name, sphere_name, semi_major,
567
0
                               0.0, "Reference_Meridian", 0.0);
568
0
            }
569
0
        }
570
571
        // translate back into a projection string.
572
3
        poDS->m_oSRS = std::move(oSRS);
573
3
    }
574
575
    /* END ISIS2 Label Read */
576
    /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
577
578
    /* -------------------------------------------------------------------- */
579
    /*      Did we get the required keywords?  If not we return with        */
580
    /*      this never having been considered to be a match. This isn't     */
581
    /*      an error!                                                       */
582
    /* -------------------------------------------------------------------- */
583
5
    if (!GDALCheckDatasetDimensions(nCols, nRows) ||
584
5
        !GDALCheckBandCount(nBands, false))
585
1
    {
586
1
        return nullptr;
587
1
    }
588
589
    /* -------------------------------------------------------------------- */
590
    /*      Capture some information from the file that is of interest.     */
591
    /* -------------------------------------------------------------------- */
592
4
    poDS->nRasterXSize = nCols;
593
4
    poDS->nRasterYSize = nRows;
594
595
    /* -------------------------------------------------------------------- */
596
    /*      Open target binary file.                                        */
597
    /* -------------------------------------------------------------------- */
598
599
4
    if (poOpenInfo->eAccess == GA_ReadOnly)
600
4
        poDS->fpImage = VSIFOpenL(osTargetFile, "rb");
601
0
    else
602
0
        poDS->fpImage = VSIFOpenL(osTargetFile, "r+b");
603
604
4
    if (poDS->fpImage == nullptr)
605
0
    {
606
0
        CPLError(CE_Failure, CPLE_OpenFailed,
607
0
                 "Failed to open %s with write permission.\n%s",
608
0
                 osTargetFile.c_str(), VSIStrerror(errno));
609
0
        return nullptr;
610
0
    }
611
612
4
    poDS->eAccess = poOpenInfo->eAccess;
613
614
    /* -------------------------------------------------------------------- */
615
    /*      Compute the line offset.                                        */
616
    /* -------------------------------------------------------------------- */
617
4
    int nItemSize = GDALGetDataTypeSizeBytes(eDataType);
618
4
    int nLineOffset, nPixelOffset;
619
4
    vsi_l_offset nBandOffset;
620
621
4
    if (EQUAL(szLayout, "BIP"))
622
0
    {
623
0
        nPixelOffset = nItemSize * nBands;
624
0
        if (nPixelOffset > INT_MAX / nBands)
625
0
        {
626
0
            return nullptr;
627
0
        }
628
0
        nLineOffset = nPixelOffset * nCols;
629
0
        nBandOffset = nItemSize;
630
0
    }
631
4
    else if (EQUAL(szLayout, "BSQ"))
632
4
    {
633
4
        nPixelOffset = nItemSize;
634
4
        if (nPixelOffset > INT_MAX / nCols)
635
0
        {
636
0
            return nullptr;
637
0
        }
638
4
        nLineOffset = nPixelOffset * nCols;
639
4
        nBandOffset = static_cast<vsi_l_offset>(nLineOffset) * nRows;
640
4
    }
641
0
    else /* assume BIL */
642
0
    {
643
0
        nPixelOffset = nItemSize;
644
0
        if (nPixelOffset > INT_MAX / nBands ||
645
0
            nPixelOffset * nBands > INT_MAX / nCols)
646
0
        {
647
0
            return nullptr;
648
0
        }
649
0
        nLineOffset = nItemSize * nBands * nCols;
650
0
        nBandOffset = static_cast<vsi_l_offset>(nItemSize) * nCols;
651
0
    }
652
653
    /* -------------------------------------------------------------------- */
654
    /*      Create band information objects.                                */
655
    /* -------------------------------------------------------------------- */
656
8
    for (int i = 0; i < nBands; i++)
657
4
    {
658
4
        auto poBand = RawRasterBand::Create(
659
4
            poDS.get(), i + 1, poDS->fpImage, nSkipBytes + nBandOffset * i,
660
4
            nPixelOffset, nLineOffset, eDataType,
661
4
            chByteOrder == 'I' || chByteOrder == 'L'
662
4
                ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
663
4
                : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN,
664
4
            RawRasterBand::OwnFP::NO);
665
4
        if (!poBand)
666
0
            return nullptr;
667
668
4
        if (bNoDataSet)
669
4
            poBand->SetNoDataValue(dfNoData);
670
671
        // Set offset/scale values at the PAM level.
672
4
        poBand->SetOffset(CPLAtofM(poDS->GetKeyword("QUBE.CORE_BASE", "0.0")));
673
4
        poBand->SetScale(
674
4
            CPLAtofM(poDS->GetKeyword("QUBE.CORE_MULTIPLIER", "1.0")));
675
676
4
        poDS->SetBand(i + 1, std::move(poBand));
677
4
    }
678
679
    /* -------------------------------------------------------------------- */
680
    /*      Check for a .prj file. For isis2 I would like to keep this in   */
681
    /* -------------------------------------------------------------------- */
682
4
    const CPLString osPath = CPLGetPathSafe(poOpenInfo->pszFilename);
683
4
    const CPLString osName = CPLGetBasenameSafe(poOpenInfo->pszFilename);
684
4
    const std::string osPrjFile = CPLFormCIFilenameSafe(osPath, osName, "prj");
685
686
4
    VSILFILE *fp = VSIFOpenL(osPrjFile.c_str(), "r");
687
4
    if (fp != nullptr)
688
0
    {
689
0
        VSIFCloseL(fp);
690
691
0
        char **papszLines = CSLLoad(osPrjFile.c_str());
692
693
0
        poDS->m_oSRS.importFromESRI(papszLines);
694
695
0
        CSLDestroy(papszLines);
696
0
    }
697
698
4
    if (dfULXMap != 0.5 || dfULYMap != 0.5 || dfXDim != 1.0 || dfYDim != 1.0)
699
3
    {
700
3
        poDS->bGotTransform = TRUE;
701
3
        poDS->m_gt[0] = dfULXMap;
702
3
        poDS->m_gt[1] = dfXDim;
703
3
        poDS->m_gt[2] = 0.0;
704
3
        poDS->m_gt[3] = dfULYMap;
705
3
        poDS->m_gt[4] = 0.0;
706
3
        poDS->m_gt[5] = dfYDim;
707
3
    }
708
709
4
    if (!poDS->bGotTransform)
710
1
        poDS->bGotTransform = GDALReadWorldFile(poOpenInfo->pszFilename, "cbw",
711
1
                                                poDS->m_gt.data());
712
713
4
    if (!poDS->bGotTransform)
714
1
        poDS->bGotTransform = GDALReadWorldFile(poOpenInfo->pszFilename, "wld",
715
1
                                                poDS->m_gt.data());
716
717
    /* -------------------------------------------------------------------- */
718
    /*      Initialize any PAM information.                                 */
719
    /* -------------------------------------------------------------------- */
720
4
    poDS->SetDescription(poOpenInfo->pszFilename);
721
4
    poDS->TryLoadXML();
722
723
    /* -------------------------------------------------------------------- */
724
    /*      Check for overviews.                                            */
725
    /* -------------------------------------------------------------------- */
726
4
    poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
727
728
4
    return poDS.release();
729
4
}
730
731
/************************************************************************/
732
/*                             GetKeyword()                             */
733
/************************************************************************/
734
735
const char *ISIS2Dataset::GetKeyword(const char *pszPath,
736
                                     const char *pszDefault)
737
738
2.17k
{
739
2.17k
    return oKeywords.GetKeyword(pszPath, pszDefault);
740
2.17k
}
741
742
/************************************************************************/
743
/*                            GetKeywordSub()                           */
744
/************************************************************************/
745
746
const char *ISIS2Dataset::GetKeywordSub(const char *pszPath, int iSubscript,
747
                                        const char *pszDefault)
748
749
2.52k
{
750
2.52k
    const char *pszResult = oKeywords.GetKeyword(pszPath, nullptr);
751
752
2.52k
    if (pszResult == nullptr)
753
2.50k
        return pszDefault;
754
755
24
    if (pszResult[0] != '(')
756
0
        return pszDefault;
757
758
24
    char **papszTokens =
759
24
        CSLTokenizeString2(pszResult, "(,)", CSLT_HONOURSTRINGS);
760
761
24
    if (iSubscript <= CSLCount(papszTokens))
762
24
    {
763
24
        oTempResult = papszTokens[iSubscript - 1];
764
24
        CSLDestroy(papszTokens);
765
24
        return oTempResult.c_str();
766
24
    }
767
768
0
    CSLDestroy(papszTokens);
769
0
    return pszDefault;
770
24
}
771
772
/************************************************************************/
773
/*                            CleanString()                             */
774
/*                                                                      */
775
/* Removes single or double quotes, and converts spaces to underscores. */
776
/* The change is made in-place to CPLString.                            */
777
/************************************************************************/
778
779
void ISIS2Dataset::CleanString(CPLString &osInput)
780
781
36
{
782
36
    if ((osInput.size() < 2) ||
783
36
        ((osInput.at(0) != '"' || osInput.back() != '"') &&
784
35
         (osInput.at(0) != '\'' || osInput.back() != '\'')))
785
5
        return;
786
787
31
    char *pszWrk = CPLStrdup(osInput.c_str() + 1);
788
789
31
    pszWrk[strlen(pszWrk) - 1] = '\0';
790
791
4.25M
    for (int i = 0; pszWrk[i] != '\0'; i++)
792
4.25M
    {
793
4.25M
        if (pszWrk[i] == ' ')
794
1.12k
            pszWrk[i] = '_';
795
4.25M
    }
796
797
31
    osInput = pszWrk;
798
31
    CPLFree(pszWrk);
799
31
}
800
801
/************************************************************************/
802
/*                         GDALRegister_ISIS2()                         */
803
/************************************************************************/
804
805
void GDALRegister_ISIS2()
806
807
22
{
808
22
    if (GDALGetDriverByName(ISIS2_DRIVER_NAME) != nullptr)
809
0
        return;
810
811
22
    GDALDriver *poDriver = new GDALDriver();
812
22
    ISIS2DriverSetCommonMetadata(poDriver);
813
814
22
    poDriver->pfnOpen = ISIS2Dataset::Open;
815
816
22
    GetGDALDriverManager()->RegisterDriver(poDriver);
817
22
}