LCOV - code coverage report
Current view: top level - frmts/raw - hkvdataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 277 382 72.5 %
Date: 2025-07-06 22:29:48 Functions: 14 19 73.7 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GView
       4             :  * Purpose:  Implementation of Atlantis HKV labelled blob support
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2000, Frank Warmerdam
       9             :  * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include <ctype.h>
      15             : 
      16             : #include "atlsci_spheroid.h"
      17             : #include "cpl_string.h"
      18             : #include "gdal_frmts.h"
      19             : #include "ogr_spatialref.h"
      20             : #include "rawdataset.h"
      21             : 
      22             : #include <cmath>
      23             : 
      24             : #include <algorithm>
      25             : 
      26             : /************************************************************************/
      27             : /* ==================================================================== */
      28             : /*                            HKVRasterBand                             */
      29             : /* ==================================================================== */
      30             : /************************************************************************/
      31             : 
      32             : class HKVDataset;
      33             : 
      34           2 : class HKVRasterBand final : public RawRasterBand
      35             : {
      36             :     friend class HKVDataset;
      37             : 
      38             :   public:
      39             :     HKVRasterBand(HKVDataset *poDS, int nBand, VSILFILE *fpRaw,
      40             :                   unsigned int nImgOffset, int nPixelOffset, int nLineOffset,
      41             :                   GDALDataType eDataType, int bNativeOrder);
      42             : 
      43             :     ~HKVRasterBand() override;
      44             : };
      45             : 
      46             : HKVRasterBand::~HKVRasterBand() = default;
      47             : 
      48             : /************************************************************************/
      49             : /*                      HKV Spheroids                                   */
      50             : /************************************************************************/
      51             : 
      52             : class HKVSpheroidList : public SpheroidList
      53             : {
      54             :   public:
      55             :     HKVSpheroidList();
      56             : 
      57           1 :     ~HKVSpheroidList()
      58           1 :     {
      59           1 :     }
      60             : };
      61             : 
      62           1 : HKVSpheroidList ::HKVSpheroidList()
      63             : {
      64           1 :     num_spheroids = 58;
      65           1 :     epsilonR = 0.1;
      66           1 :     epsilonI = 0.000001;
      67             : 
      68           1 :     spheroids[0].SetValuesByEqRadiusAndInvFlattening("airy-1830", 6377563.396,
      69             :                                                      299.3249646);
      70           1 :     spheroids[1].SetValuesByEqRadiusAndInvFlattening("modified-airy",
      71             :                                                      6377340.189, 299.3249646);
      72           1 :     spheroids[2].SetValuesByEqRadiusAndInvFlattening("australian-national",
      73             :                                                      6378160, 298.25);
      74           1 :     spheroids[3].SetValuesByEqRadiusAndInvFlattening("bessel-1841-namibia",
      75             :                                                      6377483.865, 299.1528128);
      76           1 :     spheroids[4].SetValuesByEqRadiusAndInvFlattening("bessel-1841", 6377397.155,
      77             :                                                      299.1528128);
      78           1 :     spheroids[5].SetValuesByEqRadiusAndInvFlattening("clarke-1858", 6378294.0,
      79             :                                                      294.297);
      80           1 :     spheroids[6].SetValuesByEqRadiusAndInvFlattening("clarke-1866", 6378206.4,
      81             :                                                      294.9786982);
      82           1 :     spheroids[7].SetValuesByEqRadiusAndInvFlattening("clarke-1880", 6378249.145,
      83             :                                                      293.465);
      84           1 :     spheroids[8].SetValuesByEqRadiusAndInvFlattening("everest-india-1830",
      85             :                                                      6377276.345, 300.8017);
      86           1 :     spheroids[9].SetValuesByEqRadiusAndInvFlattening("everest-sabah-sarawak",
      87             :                                                      6377298.556, 300.8017);
      88           1 :     spheroids[10].SetValuesByEqRadiusAndInvFlattening("everest-india-1956",
      89             :                                                       6377301.243, 300.8017);
      90           1 :     spheroids[11].SetValuesByEqRadiusAndInvFlattening("everest-malaysia-1969",
      91             :                                                       6377295.664, 300.8017);
      92           1 :     spheroids[12].SetValuesByEqRadiusAndInvFlattening("everest-malay-sing",
      93             :                                                       6377304.063, 300.8017);
      94           1 :     spheroids[13].SetValuesByEqRadiusAndInvFlattening("everest-pakistan",
      95             :                                                       6377309.613, 300.8017);
      96           1 :     spheroids[14].SetValuesByEqRadiusAndInvFlattening("modified-fisher-1960",
      97             :                                                       6378155, 298.3);
      98           1 :     spheroids[15].SetValuesByEqRadiusAndInvFlattening("helmert-1906", 6378200,
      99             :                                                       298.3);
     100           1 :     spheroids[16].SetValuesByEqRadiusAndInvFlattening("hough-1960", 6378270,
     101             :                                                       297);
     102           1 :     spheroids[17].SetValuesByEqRadiusAndInvFlattening("hughes", 6378273.0,
     103             :                                                       298.279);
     104           1 :     spheroids[18].SetValuesByEqRadiusAndInvFlattening("indonesian-1974",
     105             :                                                       6378160, 298.247);
     106           1 :     spheroids[19].SetValuesByEqRadiusAndInvFlattening("international-1924",
     107             :                                                       6378388, 297);
     108           1 :     spheroids[20].SetValuesByEqRadiusAndInvFlattening("iugc-67", 6378160.0,
     109             :                                                       298.254);
     110           1 :     spheroids[21].SetValuesByEqRadiusAndInvFlattening("iugc-75", 6378140.0,
     111             :                                                       298.25298);
     112           1 :     spheroids[22].SetValuesByEqRadiusAndInvFlattening("krassovsky-1940",
     113             :                                                       6378245, 298.3);
     114           1 :     spheroids[23].SetValuesByEqRadiusAndInvFlattening("kaula", 6378165.0,
     115             :                                                       292.308);
     116           1 :     spheroids[24].SetValuesByEqRadiusAndInvFlattening("grs-80", 6378137,
     117             :                                                       298.257222101);
     118           1 :     spheroids[25].SetValuesByEqRadiusAndInvFlattening("south-american-1969",
     119             :                                                       6378160, 298.25);
     120           1 :     spheroids[26].SetValuesByEqRadiusAndInvFlattening("wgs-72", 6378135,
     121             :                                                       298.26);
     122           1 :     spheroids[27].SetValuesByEqRadiusAndInvFlattening("wgs-84", 6378137,
     123             :                                                       298.257223563);
     124           1 :     spheroids[28].SetValuesByEqRadiusAndInvFlattening("ev-wgs-84", 6378137.0,
     125             :                                                       298.252841);
     126           1 :     spheroids[29].SetValuesByEqRadiusAndInvFlattening("ev-bessel", 6377397.0,
     127             :                                                       299.1976073);
     128             : 
     129           1 :     spheroids[30].SetValuesByEqRadiusAndInvFlattening("airy_1830", 6377563.396,
     130             :                                                       299.3249646);
     131           1 :     spheroids[31].SetValuesByEqRadiusAndInvFlattening("modified_airy",
     132             :                                                       6377340.189, 299.3249646);
     133           1 :     spheroids[32].SetValuesByEqRadiusAndInvFlattening("australian_national",
     134             :                                                       6378160, 298.25);
     135           1 :     spheroids[33].SetValuesByEqRadiusAndInvFlattening("bessel_1841_namibia",
     136             :                                                       6377483.865, 299.1528128);
     137           1 :     spheroids[34].SetValuesByEqRadiusAndInvFlattening("bessel_1841",
     138             :                                                       6377397.155, 299.1528128);
     139           1 :     spheroids[35].SetValuesByEqRadiusAndInvFlattening("clarke_1858", 6378294.0,
     140             :                                                       294.297);
     141           1 :     spheroids[36].SetValuesByEqRadiusAndInvFlattening("clarke_1866", 6378206.4,
     142             :                                                       294.9786982);
     143           1 :     spheroids[37].SetValuesByEqRadiusAndInvFlattening("clarke_1880",
     144             :                                                       6378249.145, 293.465);
     145           1 :     spheroids[38].SetValuesByEqRadiusAndInvFlattening("everest_india_1830",
     146             :                                                       6377276.345, 300.8017);
     147           1 :     spheroids[39].SetValuesByEqRadiusAndInvFlattening("everest_sabah_sarawak",
     148             :                                                       6377298.556, 300.8017);
     149           1 :     spheroids[40].SetValuesByEqRadiusAndInvFlattening("everest_india_1956",
     150             :                                                       6377301.243, 300.8017);
     151           1 :     spheroids[41].SetValuesByEqRadiusAndInvFlattening("everest_malaysia_1969",
     152             :                                                       6377295.664, 300.8017);
     153           1 :     spheroids[42].SetValuesByEqRadiusAndInvFlattening("everest_malay_sing",
     154             :                                                       6377304.063, 300.8017);
     155           1 :     spheroids[43].SetValuesByEqRadiusAndInvFlattening("everest_pakistan",
     156             :                                                       6377309.613, 300.8017);
     157           1 :     spheroids[44].SetValuesByEqRadiusAndInvFlattening("modified_fisher_1960",
     158             :                                                       6378155, 298.3);
     159           1 :     spheroids[45].SetValuesByEqRadiusAndInvFlattening("helmert_1906", 6378200,
     160             :                                                       298.3);
     161           1 :     spheroids[46].SetValuesByEqRadiusAndInvFlattening("hough_1960", 6378270,
     162             :                                                       297);
     163           1 :     spheroids[47].SetValuesByEqRadiusAndInvFlattening("indonesian_1974",
     164             :                                                       6378160, 298.247);
     165           1 :     spheroids[48].SetValuesByEqRadiusAndInvFlattening("international_1924",
     166             :                                                       6378388, 297);
     167           1 :     spheroids[49].SetValuesByEqRadiusAndInvFlattening("iugc_67", 6378160.0,
     168             :                                                       298.254);
     169           1 :     spheroids[50].SetValuesByEqRadiusAndInvFlattening("iugc_75", 6378140.0,
     170             :                                                       298.25298);
     171           1 :     spheroids[51].SetValuesByEqRadiusAndInvFlattening("krassovsky_1940",
     172             :                                                       6378245, 298.3);
     173           1 :     spheroids[52].SetValuesByEqRadiusAndInvFlattening("grs_80", 6378137,
     174             :                                                       298.257222101);
     175           1 :     spheroids[53].SetValuesByEqRadiusAndInvFlattening("south_american_1969",
     176             :                                                       6378160, 298.25);
     177           1 :     spheroids[54].SetValuesByEqRadiusAndInvFlattening("wgs_72", 6378135,
     178             :                                                       298.26);
     179           1 :     spheroids[55].SetValuesByEqRadiusAndInvFlattening("wgs_84", 6378137,
     180             :                                                       298.257223563);
     181           1 :     spheroids[56].SetValuesByEqRadiusAndInvFlattening("ev_wgs_84", 6378137.0,
     182             :                                                       298.252841);
     183           1 :     spheroids[57].SetValuesByEqRadiusAndInvFlattening("ev_bessel", 6377397.0,
     184             :                                                       299.1976073);
     185           1 : }
     186             : 
     187             : /************************************************************************/
     188             : /* ==================================================================== */
     189             : /*                              HKVDataset                              */
     190             : /* ==================================================================== */
     191             : /************************************************************************/
     192             : 
     193             : class HKVDataset final : public RawDataset
     194             : {
     195             :     friend class HKVRasterBand;
     196             : 
     197             :     char *pszPath;
     198             :     VSILFILE *fpBlob;
     199             : 
     200             :     int nGCPCount;
     201             :     GDAL_GCP *pasGCPList;
     202             : 
     203             :     void ProcessGeoref(const char *);
     204             :     void ProcessGeorefGCP(char **, const char *, double, double);
     205             : 
     206           1 :     void SetVersion(float version_number)
     207             :     {
     208             :         // Update stored info.
     209           1 :         MFF2version = version_number;
     210           1 :     }
     211             : 
     212             :     float MFF2version;
     213             : 
     214             :     GDALDataType eRasterType;
     215             : 
     216             :     void SetNoDataValue(double);
     217             : 
     218             :     OGRSpatialReference m_oSRS{};
     219             :     OGRSpatialReference m_oGCPSRS{};
     220             :     GDALGeoTransform m_gt{};
     221             : 
     222             :     char **papszAttrib;
     223             : 
     224             :     char **papszGeoref;
     225             : 
     226             :     // NOTE: The MFF2 format goes against GDAL's API in that nodata values are
     227             :     // set per-dataset rather than per-band.  To compromise, for writing out,
     228             :     // the dataset's nodata value will be set to the last value set on any of
     229             :     // the raster bands.
     230             : 
     231             :     bool bNoDataSet;
     232             :     double dfNoDataValue;
     233             : 
     234             :     CPL_DISALLOW_COPY_ASSIGN(HKVDataset)
     235             : 
     236             :     CPLErr Close() override;
     237             : 
     238             :   public:
     239             :     HKVDataset();
     240             :     ~HKVDataset() override;
     241             : 
     242           0 :     int GetGCPCount() override /* const */
     243             :     {
     244           0 :         return nGCPCount;
     245             :     }
     246             : 
     247           0 :     const OGRSpatialReference *GetGCPSpatialRef() const override
     248             :     {
     249           0 :         return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
     250             :     }
     251             : 
     252             :     const GDAL_GCP *GetGCPs() override;
     253             : 
     254           0 :     const OGRSpatialReference *GetSpatialRef() const override
     255             :     {
     256           0 :         return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
     257             :     }
     258             : 
     259             :     CPLErr GetGeoTransform(GDALGeoTransform &) const override;
     260             : 
     261             :     static GDALDataset *Open(GDALOpenInfo *);
     262             : };
     263             : 
     264             : /************************************************************************/
     265             : /* ==================================================================== */
     266             : /*                            HKVRasterBand                             */
     267             : /* ==================================================================== */
     268             : /************************************************************************/
     269             : 
     270             : /************************************************************************/
     271             : /*                           HKVRasterBand()                            */
     272             : /************************************************************************/
     273             : 
     274           1 : HKVRasterBand::HKVRasterBand(HKVDataset *poDSIn, int nBandIn, VSILFILE *fpRawIn,
     275             :                              unsigned int nImgOffsetIn, int nPixelOffsetIn,
     276             :                              int nLineOffsetIn, GDALDataType eDataTypeIn,
     277           1 :                              int bNativeOrderIn)
     278             :     : RawRasterBand(GDALDataset::FromHandle(poDSIn), nBandIn, fpRawIn,
     279             :                     nImgOffsetIn, nPixelOffsetIn, nLineOffsetIn, eDataTypeIn,
     280           1 :                     bNativeOrderIn, RawRasterBand::OwnFP::NO)
     281             : 
     282             : {
     283           1 :     poDS = poDSIn;
     284           1 :     nBand = nBandIn;
     285             : 
     286           1 :     nBlockXSize = poDS->GetRasterXSize();
     287           1 :     nBlockYSize = 1;
     288           1 : }
     289             : 
     290             : /************************************************************************/
     291             : /* ==================================================================== */
     292             : /*                              HKVDataset                              */
     293             : /* ==================================================================== */
     294             : /************************************************************************/
     295             : 
     296             : /************************************************************************/
     297             : /*                            HKVDataset()                             */
     298             : /************************************************************************/
     299             : 
     300           1 : HKVDataset::HKVDataset()
     301             :     : pszPath(nullptr), fpBlob(nullptr), nGCPCount(0), pasGCPList(nullptr),
     302             :       // Initialize datasets to new version; change if necessary.
     303             :       MFF2version(1.1f), eRasterType(GDT_Unknown), papszAttrib(nullptr),
     304           1 :       papszGeoref(nullptr), bNoDataSet(false), dfNoDataValue(0.0)
     305             : {
     306           1 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     307           1 :     m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     308           1 : }
     309             : 
     310             : /************************************************************************/
     311             : /*                            ~HKVDataset()                            */
     312             : /************************************************************************/
     313             : 
     314           2 : HKVDataset::~HKVDataset()
     315             : 
     316             : {
     317           1 :     HKVDataset::Close();
     318           2 : }
     319             : 
     320             : /************************************************************************/
     321             : /*                              Close()                                 */
     322             : /************************************************************************/
     323             : 
     324           2 : CPLErr HKVDataset::Close()
     325             : {
     326           2 :     CPLErr eErr = CE_None;
     327           2 :     if (nOpenFlags != OPEN_FLAGS_CLOSED)
     328             :     {
     329           1 :         if (HKVDataset::FlushCache(true) != CE_None)
     330           0 :             eErr = CE_Failure;
     331             : 
     332           1 :         if (fpBlob)
     333             :         {
     334           1 :             if (VSIFCloseL(fpBlob) != 0)
     335             :             {
     336           0 :                 eErr = CE_Failure;
     337           0 :                 CPLError(CE_Failure, CPLE_FileIO, "I/O error");
     338             :             }
     339             :         }
     340             : 
     341           1 :         if (nGCPCount > 0)
     342             :         {
     343           1 :             GDALDeinitGCPs(nGCPCount, pasGCPList);
     344           1 :             CPLFree(pasGCPList);
     345             :         }
     346             : 
     347           1 :         CPLFree(pszPath);
     348           1 :         CSLDestroy(papszGeoref);
     349           1 :         CSLDestroy(papszAttrib);
     350             : 
     351           1 :         if (GDALPamDataset::Close() != CE_None)
     352           0 :             eErr = CE_Failure;
     353             :     }
     354           2 :     return eErr;
     355             : }
     356             : 
     357             : /************************************************************************/
     358             : /*                          GetGeoTransform()                           */
     359             : /************************************************************************/
     360             : 
     361           0 : CPLErr HKVDataset::GetGeoTransform(GDALGeoTransform &gt) const
     362             : 
     363             : {
     364           0 :     gt = m_gt;
     365           0 :     return CE_None;
     366             : }
     367             : 
     368             : /************************************************************************/
     369             : /*                               GetGCP()                               */
     370             : /************************************************************************/
     371             : 
     372           0 : const GDAL_GCP *HKVDataset::GetGCPs()
     373             : 
     374             : {
     375           0 :     return pasGCPList;
     376             : }
     377             : 
     378             : /************************************************************************/
     379             : /*                          ProcessGeorefGCP()                          */
     380             : /************************************************************************/
     381             : 
     382           5 : void HKVDataset::ProcessGeorefGCP(char **papszGeorefIn, const char *pszBase,
     383             :                                   double dfRasterX, double dfRasterY)
     384             : 
     385             : {
     386             :     /* -------------------------------------------------------------------- */
     387             :     /*      Fetch the GCP from the string list.                             */
     388             :     /* -------------------------------------------------------------------- */
     389           5 :     char szFieldName[128] = {'\0'};
     390           5 :     snprintf(szFieldName, sizeof(szFieldName), "%s.latitude", pszBase);
     391           5 :     double dfLat = 0.0;
     392           5 :     if (CSLFetchNameValue(papszGeorefIn, szFieldName) == nullptr)
     393           0 :         return;
     394             :     else
     395           5 :         dfLat = CPLAtof(CSLFetchNameValue(papszGeorefIn, szFieldName));
     396             : 
     397           5 :     snprintf(szFieldName, sizeof(szFieldName), "%s.longitude", pszBase);
     398           5 :     double dfLong = 0.0;
     399           5 :     if (CSLFetchNameValue(papszGeorefIn, szFieldName) == nullptr)
     400           0 :         return;
     401             :     else
     402           5 :         dfLong = CPLAtof(CSLFetchNameValue(papszGeorefIn, szFieldName));
     403             : 
     404             :     /* -------------------------------------------------------------------- */
     405             :     /*      Add the gcp to the internal list.                               */
     406             :     /* -------------------------------------------------------------------- */
     407           5 :     GDALInitGCPs(1, pasGCPList + nGCPCount);
     408             : 
     409           5 :     CPLFree(pasGCPList[nGCPCount].pszId);
     410             : 
     411           5 :     pasGCPList[nGCPCount].pszId = CPLStrdup(pszBase);
     412             : 
     413           5 :     pasGCPList[nGCPCount].dfGCPX = dfLong;
     414           5 :     pasGCPList[nGCPCount].dfGCPY = dfLat;
     415           5 :     pasGCPList[nGCPCount].dfGCPZ = 0.0;
     416             : 
     417           5 :     pasGCPList[nGCPCount].dfGCPPixel = dfRasterX;
     418           5 :     pasGCPList[nGCPCount].dfGCPLine = dfRasterY;
     419             : 
     420           5 :     nGCPCount++;
     421             : }
     422             : 
     423             : /************************************************************************/
     424             : /*                           ProcessGeoref()                            */
     425             : /************************************************************************/
     426             : 
     427           1 : void HKVDataset::ProcessGeoref(const char *pszFilename)
     428             : 
     429             : {
     430             :     /* -------------------------------------------------------------------- */
     431             :     /*      Load the georef file, and boil white space away from around     */
     432             :     /*      the equal sign.                                                 */
     433             :     /* -------------------------------------------------------------------- */
     434           1 :     CSLDestroy(papszGeoref);
     435           1 :     papszGeoref = CSLLoad(pszFilename);
     436           1 :     if (papszGeoref == nullptr)
     437           0 :         return;
     438             : 
     439           1 :     HKVSpheroidList *hkvEllipsoids = new HKVSpheroidList;
     440             : 
     441          14 :     for (int i = 0; papszGeoref[i] != nullptr; i++)
     442             :     {
     443          13 :         int iDst = 0;
     444          13 :         char *pszLine = papszGeoref[i];
     445             : 
     446         433 :         for (int iSrc = 0; pszLine[iSrc] != '\0'; iSrc++)
     447             :         {
     448         420 :             if (pszLine[iSrc] != ' ')
     449             :             {
     450         420 :                 pszLine[iDst++] = pszLine[iSrc];
     451             :             }
     452             :         }
     453          13 :         pszLine[iDst] = '\0';
     454             :     }
     455             : 
     456             :     /* -------------------------------------------------------------------- */
     457             :     /*      Try to get GCPs, in lat/longs                     .             */
     458             :     /* -------------------------------------------------------------------- */
     459           1 :     nGCPCount = 0;
     460           1 :     pasGCPList = reinterpret_cast<GDAL_GCP *>(CPLCalloc(sizeof(GDAL_GCP), 5));
     461             : 
     462           1 :     if (MFF2version > 1.0)
     463             :     {
     464           1 :         ProcessGeorefGCP(papszGeoref, "top_left", 0, 0);
     465           1 :         ProcessGeorefGCP(papszGeoref, "top_right", GetRasterXSize(), 0);
     466           1 :         ProcessGeorefGCP(papszGeoref, "bottom_left", 0, GetRasterYSize());
     467           1 :         ProcessGeorefGCP(papszGeoref, "bottom_right", GetRasterXSize(),
     468           1 :                          GetRasterYSize());
     469           1 :         ProcessGeorefGCP(papszGeoref, "centre", GetRasterXSize() / 2.0,
     470           1 :                          GetRasterYSize() / 2.0);
     471             :     }
     472             :     else
     473             :     {
     474           0 :         ProcessGeorefGCP(papszGeoref, "top_left", 0.5, 0.5);
     475           0 :         ProcessGeorefGCP(papszGeoref, "top_right", GetRasterXSize() - 0.5, 0.5);
     476           0 :         ProcessGeorefGCP(papszGeoref, "bottom_left", 0.5,
     477           0 :                          GetRasterYSize() - 0.5);
     478           0 :         ProcessGeorefGCP(papszGeoref, "bottom_right", GetRasterXSize() - 0.5,
     479           0 :                          GetRasterYSize() - 0.5);
     480           0 :         ProcessGeorefGCP(papszGeoref, "centre", GetRasterXSize() / 2.0,
     481           0 :                          GetRasterYSize() / 2.0);
     482             :     }
     483             : 
     484           1 :     if (nGCPCount == 0)
     485             :     {
     486           0 :         CPLFree(pasGCPList);
     487           0 :         pasGCPList = nullptr;
     488             :     }
     489             : 
     490             :     /* -------------------------------------------------------------------- */
     491             :     /*      Do we have a recognised projection?                             */
     492             :     /* -------------------------------------------------------------------- */
     493           1 :     const char *pszProjName = CSLFetchNameValue(papszGeoref, "projection.name");
     494             :     const char *pszOriginLong =
     495           1 :         CSLFetchNameValue(papszGeoref, "projection.origin_longitude");
     496             :     const char *pszSpheroidName =
     497           1 :         CSLFetchNameValue(papszGeoref, "spheroid.name");
     498             : 
     499           2 :     if (pszSpheroidName != nullptr &&
     500           1 :         hkvEllipsoids->SpheroidInList(pszSpheroidName))
     501             :     {
     502             : #if 0
     503             :       // TODO(schwehr): Enable in trunk after 2.1 branch and fix.
     504             :       // Breaks tests on some platforms.
     505             :       CPLError( CE_Failure, CPLE_AppDefined,
     506             :                 "Unrecognized ellipsoid.  Not handled.  "
     507             :                 "Spheroid name not in spheroid list: '%s'",
     508             :                 pszSpheroidName );
     509             : #endif
     510             :         // Why were eq_radius and inv_flattening never used?
     511             :         // eq_radius = hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName);
     512             :         // inv_flattening =
     513             :         //     hkvEllipsoids->GetSpheroidInverseFlattening(pszSpheroidName);
     514             :     }
     515           0 :     else if (pszProjName != nullptr)
     516             :     {
     517           0 :         CPLError(CE_Warning, CPLE_AppDefined,
     518             :                  "Unrecognized ellipsoid.  Not handled.");
     519             :         // TODO(schwehr): This error is was never what was happening.
     520             :         // CPLError( CE_Warning, CPLE_AppDefined,
     521             :         //           "Unrecognized ellipsoid.  Using wgs-84 parameters.");
     522             :         // eq_radius=hkvEllipsoids->GetSpheroidEqRadius("wgs-84"); */
     523             :         // inv_flattening=hkvEllipsoids->GetSpheroidInverseFlattening("wgs-84");
     524             :     }
     525             : 
     526           1 :     if (pszProjName != nullptr && EQUAL(pszProjName, "utm") && nGCPCount == 5)
     527             :     {
     528             :         // int nZone = (int)((CPLAtof(pszOriginLong)+184.5) / 6.0);
     529           1 :         int nZone = 31;  // TODO(schwehr): Where does 31 come from?
     530             : 
     531           1 :         if (pszOriginLong == nullptr)
     532             :         {
     533             :             // If origin not specified, assume 0.0.
     534           0 :             CPLError(
     535             :                 CE_Warning, CPLE_AppDefined,
     536             :                 "No projection origin longitude specified.  Assuming 0.0.");
     537             :         }
     538             :         else
     539             :         {
     540           1 :             nZone = 31 + static_cast<int>(floor(CPLAtof(pszOriginLong) / 6.0));
     541             :         }
     542             : 
     543           2 :         OGRSpatialReference oUTM;
     544             : 
     545           1 :         if (pasGCPList[4].dfGCPY < 0)
     546           0 :             oUTM.SetUTM(nZone, 0);
     547             :         else
     548           1 :             oUTM.SetUTM(nZone, 1);
     549             : 
     550           2 :         OGRSpatialReference oLL;
     551           1 :         oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     552           1 :         if (pszOriginLong != nullptr)
     553             :         {
     554           1 :             oUTM.SetProjParm(SRS_PP_CENTRAL_MERIDIAN, CPLAtof(pszOriginLong));
     555           1 :             oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN, CPLAtof(pszOriginLong));
     556             :         }
     557             : 
     558           1 :         if ((pszSpheroidName == nullptr) ||
     559           1 :             (EQUAL(pszSpheroidName, "wgs-84")) ||
     560           1 :             (EQUAL(pszSpheroidName, "wgs_84")))
     561             :         {
     562           0 :             oUTM.SetWellKnownGeogCS("WGS84");
     563           0 :             oLL.SetWellKnownGeogCS("WGS84");
     564             :         }
     565             :         else
     566             :         {
     567           1 :             if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
     568             :             {
     569           1 :                 oUTM.SetGeogCS(
     570             :                     "unknown", "unknown", pszSpheroidName,
     571             :                     hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
     572             :                     hkvEllipsoids->GetSpheroidInverseFlattening(
     573             :                         pszSpheroidName));
     574           1 :                 oLL.SetGeogCS(
     575             :                     "unknown", "unknown", pszSpheroidName,
     576             :                     hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
     577             :                     hkvEllipsoids->GetSpheroidInverseFlattening(
     578             :                         pszSpheroidName));
     579             :             }
     580             :             else
     581             :             {
     582           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
     583             :                          "Unrecognized ellipsoid.  Using wgs-84 parameters.");
     584           0 :                 oUTM.SetWellKnownGeogCS("WGS84");
     585           0 :                 oLL.SetWellKnownGeogCS("WGS84");
     586             :             }
     587             :         }
     588             : 
     589             :         OGRCoordinateTransformation *poTransform =
     590           1 :             OGRCreateCoordinateTransformation(&oLL, &oUTM);
     591             : 
     592           1 :         bool bSuccess = true;
     593           1 :         if (poTransform == nullptr)
     594             :         {
     595           0 :             CPLErrorReset();
     596           0 :             bSuccess = false;
     597             :         }
     598             : 
     599           1 :         double dfUtmX[5] = {0.0};
     600           1 :         double dfUtmY[5] = {0.0};
     601             : 
     602           1 :         if (poTransform != nullptr)
     603             :         {
     604           6 :             for (int gcp_index = 0; gcp_index < 5; gcp_index++)
     605             :             {
     606           5 :                 dfUtmX[gcp_index] = pasGCPList[gcp_index].dfGCPX;
     607           5 :                 dfUtmY[gcp_index] = pasGCPList[gcp_index].dfGCPY;
     608             : 
     609           5 :                 if (bSuccess && !poTransform->Transform(1, &(dfUtmX[gcp_index]),
     610             :                                                         &(dfUtmY[gcp_index])))
     611           0 :                     bSuccess = false;
     612             :             }
     613             :         }
     614             : 
     615           1 :         if (bSuccess)
     616             :         {
     617             :             // Update GCPS to proper projection.
     618           6 :             for (int gcp_index = 0; gcp_index < 5; gcp_index++)
     619             :             {
     620           5 :                 pasGCPList[gcp_index].dfGCPX = dfUtmX[gcp_index];
     621           5 :                 pasGCPList[gcp_index].dfGCPY = dfUtmY[gcp_index];
     622             :             }
     623             : 
     624           1 :             m_oGCPSRS = oUTM;
     625             : 
     626           2 :             bool transform_ok = CPL_TO_BOOL(
     627           1 :                 GDALGCPsToGeoTransform(5, pasGCPList, m_gt.data(), 0));
     628             : 
     629           1 :             if (!transform_ok)
     630             :             {
     631             :                 // Transform may not be sufficient in all cases (slant range
     632             :                 // projection).
     633           0 :                 m_gt = GDALGeoTransform();
     634           0 :                 m_oGCPSRS.Clear();
     635             :             }
     636             :             else
     637             :             {
     638           1 :                 m_oSRS = std::move(oUTM);
     639             :             }
     640             :         }
     641             : 
     642           1 :         if (poTransform != nullptr)
     643           2 :             delete poTransform;
     644             :     }
     645           0 :     else if (pszProjName != nullptr && nGCPCount == 5)
     646             :     {
     647           0 :         OGRSpatialReference oLL;
     648           0 :         oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     649             : 
     650           0 :         if (pszOriginLong != nullptr)
     651             :         {
     652           0 :             oLL.SetProjParm(SRS_PP_LONGITUDE_OF_ORIGIN, CPLAtof(pszOriginLong));
     653             :         }
     654             : 
     655           0 :         if (pszSpheroidName == nullptr ||
     656           0 :             EQUAL(pszSpheroidName, "wgs-84") ||  // Dash.
     657           0 :             EQUAL(pszSpheroidName, "wgs_84"))    // Underscore.
     658             :         {
     659           0 :             oLL.SetWellKnownGeogCS("WGS84");
     660             :         }
     661             :         else
     662             :         {
     663           0 :             if (hkvEllipsoids->SpheroidInList(pszSpheroidName))
     664             :             {
     665           0 :                 oLL.SetGeogCS(
     666             :                     "", "", pszSpheroidName,
     667             :                     hkvEllipsoids->GetSpheroidEqRadius(pszSpheroidName),
     668             :                     hkvEllipsoids->GetSpheroidInverseFlattening(
     669             :                         pszSpheroidName));
     670             :             }
     671             :             else
     672             :             {
     673           0 :                 CPLError(CE_Warning, CPLE_AppDefined,
     674             :                          "Unrecognized ellipsoid.  "
     675             :                          "Using wgs-84 parameters.");
     676           0 :                 oLL.SetWellKnownGeogCS("WGS84");
     677             :             }
     678             :         }
     679             : 
     680             :         const bool transform_ok =
     681           0 :             CPL_TO_BOOL(GDALGCPsToGeoTransform(5, pasGCPList, m_gt.data(), 0));
     682             : 
     683           0 :         m_oSRS.Clear();
     684             : 
     685           0 :         if (!transform_ok)
     686             :         {
     687           0 :             m_gt = GDALGeoTransform();
     688             :         }
     689             :         else
     690             :         {
     691           0 :             m_oSRS = oLL;
     692             :         }
     693             : 
     694           0 :         m_oGCPSRS = std::move(oLL);
     695             :     }
     696             : 
     697           1 :     delete hkvEllipsoids;
     698             : }
     699             : 
     700             : /************************************************************************/
     701             : /*                                Open()                                */
     702             : /************************************************************************/
     703             : 
     704       33531 : GDALDataset *HKVDataset::Open(GDALOpenInfo *poOpenInfo)
     705             : 
     706             : {
     707             :     /* -------------------------------------------------------------------- */
     708             :     /*      We assume the dataset is passed as a directory.  Check for      */
     709             :     /*      an attrib and blob file as a minimum.                           */
     710             :     /* -------------------------------------------------------------------- */
     711       33531 :     if (!poOpenInfo->bIsDirectory)
     712       33178 :         return nullptr;
     713             : 
     714             :     std::string osFilename =
     715         706 :         CPLFormFilenameSafe(poOpenInfo->pszFilename, "image_data", nullptr);
     716             :     VSIStatBuf sStat;
     717         353 :     if (VSIStat(osFilename.c_str(), &sStat) != 0)
     718             :         osFilename =
     719         352 :             CPLFormFilenameSafe(poOpenInfo->pszFilename, "blob", nullptr);
     720         353 :     if (VSIStat(osFilename.c_str(), &sStat) != 0)
     721         352 :         return nullptr;
     722             : 
     723             :     osFilename =
     724           1 :         CPLFormFilenameSafe(poOpenInfo->pszFilename, "attrib", nullptr);
     725           1 :     if (VSIStat(osFilename.c_str(), &sStat) != 0)
     726           0 :         return nullptr;
     727             : 
     728             :     /* -------------------------------------------------------------------- */
     729             :     /*      Load the attrib file, and boil white space away from around     */
     730             :     /*      the equal sign.                                                 */
     731             :     /* -------------------------------------------------------------------- */
     732           1 :     char **papszAttrib = CSLLoad(osFilename.c_str());
     733           1 :     if (papszAttrib == nullptr)
     734           0 :         return nullptr;
     735             : 
     736          10 :     for (int i = 0; papszAttrib[i] != nullptr; i++)
     737             :     {
     738           9 :         int iDst = 0;
     739           9 :         char *pszLine = papszAttrib[i];
     740             : 
     741         252 :         for (int iSrc = 0; pszLine[iSrc] != '\0'; iSrc++)
     742             :         {
     743         243 :             if (pszLine[iSrc] != ' ')
     744             :             {
     745         211 :                 pszLine[iDst++] = pszLine[iSrc];
     746             :             }
     747             :         }
     748           9 :         pszLine[iDst] = '\0';
     749             :     }
     750             : 
     751             :     /* -------------------------------------------------------------------- */
     752             :     /*      Create a corresponding GDALDataset.                             */
     753             :     /* -------------------------------------------------------------------- */
     754           2 :     auto poDS = std::make_unique<HKVDataset>();
     755             : 
     756           1 :     poDS->pszPath = CPLStrdup(poOpenInfo->pszFilename);
     757           1 :     poDS->papszAttrib = papszAttrib;
     758             : 
     759           1 :     poDS->eAccess = poOpenInfo->eAccess;
     760             : 
     761             :     /* -------------------------------------------------------------------- */
     762             :     /*      Set some dataset wide information.                              */
     763             :     /* -------------------------------------------------------------------- */
     764           1 :     bool bNative = false;
     765           1 :     bool bComplex = false;
     766           1 :     int nRawBands = 0;
     767             : 
     768           2 :     if (CSLFetchNameValue(papszAttrib, "extent.cols") == nullptr ||
     769           1 :         CSLFetchNameValue(papszAttrib, "extent.rows") == nullptr)
     770             :     {
     771           0 :         return nullptr;
     772             :     }
     773             : 
     774           1 :     poDS->nRasterXSize = atoi(CSLFetchNameValue(papszAttrib, "extent.cols"));
     775           1 :     poDS->nRasterYSize = atoi(CSLFetchNameValue(papszAttrib, "extent.rows"));
     776             : 
     777           1 :     if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
     778             :     {
     779           0 :         return nullptr;
     780             :     }
     781             : 
     782           1 :     const char *pszValue = CSLFetchNameValue(papszAttrib, "pixel.order");
     783           1 :     if (pszValue == nullptr)
     784           0 :         bNative = true;
     785             :     else
     786             :     {
     787             : #ifdef CPL_MSB
     788             :         bNative = strstr(pszValue, "*msbf") != NULL;
     789             : #else
     790           1 :         bNative = strstr(pszValue, "*lsbf") != nullptr;
     791             : #endif
     792             :     }
     793             : 
     794           1 :     bool bNoDataSet = false;
     795           1 :     double dfNoDataValue = 0.0;
     796           1 :     pszValue = CSLFetchNameValue(papszAttrib, "pixel.no_data");
     797           1 :     if (pszValue != nullptr)
     798             :     {
     799           0 :         bNoDataSet = true;
     800           0 :         dfNoDataValue = CPLAtof(pszValue);
     801             :     }
     802             : 
     803           1 :     pszValue = CSLFetchNameValue(papszAttrib, "channel.enumeration");
     804           1 :     if (pszValue != nullptr)
     805           1 :         nRawBands = atoi(pszValue);
     806             :     else
     807           0 :         nRawBands = 1;
     808             : 
     809           1 :     if (!GDALCheckBandCount(nRawBands, TRUE))
     810             :     {
     811           0 :         return nullptr;
     812             :     }
     813             : 
     814           1 :     pszValue = CSLFetchNameValue(papszAttrib, "pixel.field");
     815           1 :     if (pszValue != nullptr && strstr(pszValue, "*complex") != nullptr)
     816           0 :         bComplex = true;
     817             :     else
     818           1 :         bComplex = false;
     819             : 
     820             :     /* Get the version number, if present (if not, assume old version. */
     821             :     /* Versions differ in their interpretation of corner coordinates.  */
     822             : 
     823           1 :     if (CSLFetchNameValue(papszAttrib, "version") != nullptr)
     824           1 :         poDS->SetVersion(static_cast<float>(
     825           1 :             CPLAtof(CSLFetchNameValue(papszAttrib, "version"))));
     826             :     else
     827           0 :         poDS->SetVersion(1.0);
     828             : 
     829             :     /* -------------------------------------------------------------------- */
     830             :     /*      Figure out the datatype                                         */
     831             :     /* -------------------------------------------------------------------- */
     832           1 :     const char *pszEncoding = CSLFetchNameValue(papszAttrib, "pixel.encoding");
     833           1 :     if (pszEncoding == nullptr)
     834           0 :         pszEncoding = "{ *unsigned }";
     835             : 
     836           1 :     int nSize = 1;
     837           1 :     if (CSLFetchNameValue(papszAttrib, "pixel.size") != nullptr)
     838           1 :         nSize = atoi(CSLFetchNameValue(papszAttrib, "pixel.size")) / 8;
     839             : #if 0
     840             :     int nPseudoBands;
     841             :     if( bComplex )
     842             :         nPseudoBands = 2;
     843             :     else
     844             :         nPseudoBands = 1;
     845             : #endif
     846             : 
     847             :     GDALDataType eType;
     848           1 :     if (nSize == 1)
     849           1 :         eType = GDT_Byte;
     850           0 :     else if (nSize == 2 && strstr(pszEncoding, "*unsigned") != nullptr)
     851           0 :         eType = GDT_UInt16;
     852           0 :     else if (nSize == 4 && bComplex)
     853           0 :         eType = GDT_CInt16;
     854           0 :     else if (nSize == 2)
     855           0 :         eType = GDT_Int16;
     856           0 :     else if (nSize == 4 && strstr(pszEncoding, "*unsigned") != nullptr)
     857           0 :         eType = GDT_UInt32;
     858           0 :     else if (nSize == 8 && strstr(pszEncoding, "*two") != nullptr && bComplex)
     859           0 :         eType = GDT_CInt32;
     860           0 :     else if (nSize == 4 && strstr(pszEncoding, "*two") != nullptr)
     861           0 :         eType = GDT_Int32;
     862           0 :     else if (nSize == 8 && bComplex)
     863           0 :         eType = GDT_CFloat32;
     864           0 :     else if (nSize == 4)
     865           0 :         eType = GDT_Float32;
     866           0 :     else if (nSize == 16 && bComplex)
     867           0 :         eType = GDT_CFloat64;
     868           0 :     else if (nSize == 8)
     869           0 :         eType = GDT_Float64;
     870             :     else
     871             :     {
     872           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     873             :                  "Unsupported pixel data type in %s.\n"
     874             :                  "pixel.size=%d pixel.encoding=%s",
     875           0 :                  poDS->pszPath, nSize, pszEncoding);
     876           0 :         return nullptr;
     877             :     }
     878             : 
     879             :     /* -------------------------------------------------------------------- */
     880             :     /*      Open the blob file.                                             */
     881             :     /* -------------------------------------------------------------------- */
     882           1 :     osFilename = CPLFormFilenameSafe(poDS->pszPath, "image_data", nullptr);
     883           1 :     if (VSIStat(osFilename.c_str(), &sStat) != 0)
     884           0 :         osFilename = CPLFormFilenameSafe(poDS->pszPath, "blob", nullptr);
     885           1 :     if (poOpenInfo->eAccess == GA_ReadOnly)
     886             :     {
     887           1 :         poDS->fpBlob = VSIFOpenL(osFilename.c_str(), "rb");
     888           1 :         if (poDS->fpBlob == nullptr)
     889             :         {
     890           0 :             CPLError(CE_Failure, CPLE_OpenFailed,
     891             :                      "Unable to open file %s for read access.",
     892             :                      osFilename.c_str());
     893           0 :             return nullptr;
     894             :         }
     895             :     }
     896             :     else
     897             :     {
     898           0 :         poDS->fpBlob = VSIFOpenL(osFilename.c_str(), "rb+");
     899           0 :         if (poDS->fpBlob == nullptr)
     900             :         {
     901           0 :             CPLError(CE_Failure, CPLE_OpenFailed,
     902             :                      "Unable to open file %s for update access.",
     903             :                      osFilename.c_str());
     904           0 :             return nullptr;
     905             :         }
     906             :     }
     907             : 
     908             :     /* -------------------------------------------------------------------- */
     909             :     /*      Build the overview filename, as blob file = "_ovr".             */
     910             :     /* -------------------------------------------------------------------- */
     911           2 :     std::string osOvrFilename(osFilename);
     912           1 :     osOvrFilename += "_ovr";
     913             : 
     914             :     /* -------------------------------------------------------------------- */
     915             :     /*      Define the bands.                                               */
     916             :     /* -------------------------------------------------------------------- */
     917           1 :     const int nPixelOffset = nRawBands * nSize;
     918           1 :     const int nLineOffset = nPixelOffset * poDS->GetRasterXSize();
     919           1 :     int nOffset = 0;
     920             : 
     921           2 :     for (int iRawBand = 0; iRawBand < nRawBands; iRawBand++)
     922             :     {
     923             :         auto poBand = std::make_unique<HKVRasterBand>(
     924           1 :             poDS.get(), poDS->GetRasterCount() + 1, poDS->fpBlob, nOffset,
     925           1 :             nPixelOffset, nLineOffset, eType, bNative);
     926           1 :         if (!poBand->IsValid())
     927           0 :             return nullptr;
     928             : 
     929           1 :         if (bNoDataSet)
     930           0 :             poBand->SetNoDataValue(dfNoDataValue);
     931           1 :         poDS->SetBand(poDS->GetRasterCount() + 1, std::move(poBand));
     932           1 :         nOffset += GDALGetDataTypeSizeBytes(eType);
     933             :     }
     934             : 
     935           1 :     poDS->eRasterType = eType;
     936             : 
     937             :     /* -------------------------------------------------------------------- */
     938             :     /*      Process the georef file if there is one.                        */
     939             :     /* -------------------------------------------------------------------- */
     940           1 :     osFilename = CPLFormFilenameSafe(poDS->pszPath, "georef", nullptr);
     941           1 :     if (VSIStat(osFilename.c_str(), &sStat) == 0)
     942           1 :         poDS->ProcessGeoref(osFilename.c_str());
     943             : 
     944             :     /* -------------------------------------------------------------------- */
     945             :     /*      Initialize any PAM information.                                 */
     946             :     /* -------------------------------------------------------------------- */
     947           1 :     poDS->SetDescription(osOvrFilename.c_str());
     948           1 :     poDS->TryLoadXML();
     949             : 
     950             :     /* -------------------------------------------------------------------- */
     951             :     /*      Handle overviews.                                               */
     952             :     /* -------------------------------------------------------------------- */
     953           1 :     poDS->oOvManager.Initialize(poDS.get(), osOvrFilename.c_str(), nullptr,
     954             :                                 TRUE);
     955             : 
     956           1 :     return poDS.release();
     957             : }
     958             : 
     959             : /************************************************************************/
     960             : /*                         GDALRegister_HKV()                           */
     961             : /************************************************************************/
     962             : 
     963        1928 : void GDALRegister_HKV()
     964             : 
     965             : {
     966        1928 :     if (GDALGetDriverByName("MFF2") != nullptr)
     967         282 :         return;
     968             : 
     969        1646 :     GDALDriver *poDriver = new GDALDriver();
     970             : 
     971        1646 :     poDriver->SetDescription("MFF2");
     972        1646 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     973        1646 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Vexcel MFF2 (HKV) Raster");
     974        1646 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/mff2.html");
     975             : 
     976        1646 :     poDriver->pfnOpen = HKVDataset::Open;
     977             : 
     978        1646 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     979             : }

Generated by: LCOV version 1.14