Line data    Source code

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