Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  ACE2 Driver
       4             :  * Purpose:  Implementation of ACE2 elevation format read support.
       5             :  *           http://tethys.eaprs.cse.dmu.ac.uk/ACE2/shared/documentation
       6             :  * Author:   Even Rouault, <even dot rouault at spatialys.com>
       7             :  *
       8             :  ******************************************************************************
       9             :  * Copyright (c) 2011-2012, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "gdal_frmts.h"
      15             : #include "ogr_spatialref.h"
      16             : #include "rawdataset.h"
      17             : 
      18             : static const char *const apszCategorySource[] = {
      19             :     "Pure SRTM (above 60deg N pure GLOBE data, below 60S pure ACE [original] "
      20             :     "data)",
      21             :     "SRTM voids filled by interpolation and/or altimeter data",
      22             :     "SRTM data warped using the ERS-1 Geodetic Mission",
      23             :     "SRTM data warped using EnviSat & ERS-2 data",
      24             :     "Mean lake level data derived from Altimetry",
      25             :     "GLOBE/ACE data warped using combined altimetry (only above 60deg N)",
      26             :     "Pure altimetry data (derived from ERS-1 Geodetic Mission, ERS-2 and "
      27             :     "EnviSat data using Delaunay Triangulation",
      28             :     nullptr};
      29             : 
      30             : static const char *const apszCategoryQuality[] = {
      31             :     "Generic - use base datasets",
      32             :     "Accuracy of greater than +/- 16m",
      33             :     "Accuracy between +/- 16m - +/- 10m",
      34             :     "Accuracy between +/-10m - +/-5m",
      35             :     "Accuracy between +/-5m - +/-1m",
      36             :     "Accuracy between +/-1m",
      37             :     nullptr};
      38             : 
      39             : static const char *const apszCategoryConfidence[] = {
      40             :     "No confidence could be derived due to lack of data",
      41             :     "Heights generated by interpolation",
      42             :     "Low confidence",
      43             :     "Low confidence",
      44             :     "Low confidence",
      45             :     "Medium confidence",
      46             :     "Medium confidence",
      47             :     "Medium confidence",
      48             :     "Medium confidence",
      49             :     "Medium confidence",
      50             :     "Medium confidence",
      51             :     "Medium confidence",
      52             :     "Medium confidence",
      53             :     "High confidence",
      54             :     "High confidence",
      55             :     "High confidence",
      56             :     "High confidence",
      57             :     "Inland water confidence",
      58             :     "Inland water confidence",
      59             :     "Inland water confidence",
      60             :     "Inland water confidence",
      61             :     "Inland water confidence",
      62             :     nullptr};
      63             : 
      64             : /************************************************************************/
      65             : /* ==================================================================== */
      66             : /*                             ACE2Dataset                              */
      67             : /* ==================================================================== */
      68             : /************************************************************************/
      69             : 
      70             : class ACE2Dataset final : public GDALPamDataset
      71             : {
      72             :     friend class ACE2RasterBand;
      73             : 
      74             :     OGRSpatialReference m_oSRS{};
      75             :     double adfGeoTransform[6];
      76             : 
      77             :   public:
      78             :     ACE2Dataset();
      79             : 
      80           1 :     const OGRSpatialReference *GetSpatialRef() const override
      81             :     {
      82           1 :         return &m_oSRS;
      83             :     }
      84             : 
      85             :     CPLErr GetGeoTransform(double *) override;
      86             : 
      87             :     static GDALDataset *Open(GDALOpenInfo *);
      88             :     static int Identify(GDALOpenInfo *);
      89             : };
      90             : 
      91             : /************************************************************************/
      92             : /* ==================================================================== */
      93             : /*                            ACE2RasterBand                            */
      94             : /* ==================================================================== */
      95             : /************************************************************************/
      96             : 
      97             : class ACE2RasterBand final : public RawRasterBand
      98             : {
      99             :   public:
     100             :     ACE2RasterBand(VSILFILE *fpRaw, GDALDataType eDataType, int nXSize,
     101             :                    int nYSize);
     102             : 
     103             :     const char *GetUnitType() override;
     104             :     char **GetCategoryNames() override;
     105             : };
     106             : 
     107             : /************************************************************************/
     108             : /* ==================================================================== */
     109             : /*                             ACE2Dataset                              */
     110             : /* ==================================================================== */
     111             : /************************************************************************/
     112             : 
     113             : /************************************************************************/
     114             : /*                            ACE2Dataset()                             */
     115             : /************************************************************************/
     116             : 
     117           2 : ACE2Dataset::ACE2Dataset()
     118             : {
     119           2 :     m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
     120           2 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     121             : 
     122           2 :     adfGeoTransform[0] = 0.0;
     123           2 :     adfGeoTransform[1] = 1.0;
     124           2 :     adfGeoTransform[2] = 0.0;
     125           2 :     adfGeoTransform[3] = 0.0;
     126           2 :     adfGeoTransform[4] = 0.0;
     127           2 :     adfGeoTransform[5] = 1.0;
     128           2 : }
     129             : 
     130             : /************************************************************************/
     131             : /*                          GetGeoTransform()                           */
     132             : /************************************************************************/
     133             : 
     134           1 : CPLErr ACE2Dataset::GetGeoTransform(double *padfTransform)
     135             : 
     136             : {
     137           1 :     memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
     138           1 :     return CE_None;
     139             : }
     140             : 
     141             : /************************************************************************/
     142             : /*                          ACE2RasterBand()                            */
     143             : /************************************************************************/
     144             : 
     145           2 : ACE2RasterBand::ACE2RasterBand(VSILFILE *fpRawIn, GDALDataType eDataTypeIn,
     146           2 :                                int nXSize, int nYSize)
     147             :     : RawRasterBand(fpRawIn, 0, GDALGetDataTypeSizeBytes(eDataTypeIn),
     148           4 :                     nXSize * GDALGetDataTypeSizeBytes(eDataTypeIn), eDataTypeIn,
     149           2 :                     CPL_IS_LSB, nXSize, nYSize, RawRasterBand::OwnFP::YES)
     150             : {
     151           2 : }
     152             : 
     153             : /************************************************************************/
     154             : /*                             GetUnitType()                            */
     155             : /************************************************************************/
     156             : 
     157           0 : const char *ACE2RasterBand::GetUnitType()
     158             : {
     159           0 :     if (eDataType == GDT_Float32)
     160           0 :         return "m";
     161             : 
     162           0 :     return "";
     163             : }
     164             : 
     165             : /************************************************************************/
     166             : /*                         GetCategoryNames()                           */
     167             : /************************************************************************/
     168             : 
     169           0 : char **ACE2RasterBand::GetCategoryNames()
     170             : {
     171           0 :     if (eDataType != GDT_Int16)
     172           0 :         return nullptr;
     173             : 
     174           0 :     const char *pszName = poDS->GetDescription();
     175             : 
     176           0 :     if (strstr(pszName, "_SOURCE_"))
     177           0 :         return const_cast<char **>(apszCategorySource);
     178           0 :     if (strstr(pszName, "_QUALITY_"))
     179           0 :         return const_cast<char **>(apszCategoryQuality);
     180           0 :     if (strstr(pszName, "_CONF_"))
     181           0 :         return const_cast<char **>(apszCategoryConfidence);
     182             : 
     183           0 :     return nullptr;
     184             : }
     185             : 
     186             : /************************************************************************/
     187             : /*                             Identify()                               */
     188             : /************************************************************************/
     189             : 
     190       51977 : int ACE2Dataset::Identify(GDALOpenInfo *poOpenInfo)
     191             : 
     192             : {
     193      103950 :     if (!(poOpenInfo->IsExtensionEqualToCI("ACE2") ||
     194       51973 :           strstr(poOpenInfo->pszFilename, ".ACE2.gz") ||
     195       51974 :           strstr(poOpenInfo->pszFilename, ".ace2.gz")))
     196       51973 :         return FALSE;
     197             : 
     198           4 :     return TRUE;
     199             : }
     200             : 
     201             : /************************************************************************/
     202             : /*                                Open()                                */
     203             : /************************************************************************/
     204             : 
     205           2 : GDALDataset *ACE2Dataset::Open(GDALOpenInfo *poOpenInfo)
     206             : 
     207             : {
     208           2 :     if (!Identify(poOpenInfo))
     209           0 :         return nullptr;
     210             : 
     211           4 :     const std::string osBasename = CPLGetBasenameSafe(poOpenInfo->pszFilename);
     212           2 :     const char *pszBasename = osBasename.c_str();
     213             : 
     214           2 :     if (strlen(pszBasename) < 7)
     215           0 :         return nullptr;
     216             : 
     217             :     /* Determine southwest coordinates from filename */
     218             : 
     219             :     /* e.g. 30S120W_5M.ACE2 */
     220           2 :     char pszLatLonValueString[4] = {'\0'};
     221           2 :     memset(pszLatLonValueString, 0, 4);
     222             :     // cppcheck-suppress redundantCopy
     223           2 :     strncpy(pszLatLonValueString, &pszBasename[0], 2);
     224           2 :     int southWestLat = atoi(pszLatLonValueString);
     225           2 :     memset(pszLatLonValueString, 0, 4);
     226             :     // cppcheck-suppress redundantCopy
     227           2 :     strncpy(pszLatLonValueString, &pszBasename[3], 3);
     228           2 :     int southWestLon = atoi(pszLatLonValueString);
     229             : 
     230           2 :     if (pszBasename[2] == 'N' || pszBasename[2] == 'n')
     231             :         /*southWestLat = southWestLat*/;
     232           0 :     else if (pszBasename[2] == 'S' || pszBasename[2] == 's')
     233           0 :         southWestLat = southWestLat * -1;
     234             :     else
     235           0 :         return nullptr;
     236             : 
     237           2 :     if (pszBasename[6] == 'E' || pszBasename[6] == 'e')
     238             :         /*southWestLon = southWestLon*/;
     239           0 :     else if (pszBasename[6] == 'W' || pszBasename[6] == 'w')
     240           0 :         southWestLon = southWestLon * -1;
     241             :     else
     242           0 :         return nullptr;
     243             : 
     244           2 :     GDALDataType eDT = GDT_Unknown;
     245           2 :     if (strstr(pszBasename, "_CONF_") || strstr(pszBasename, "_QUALITY_") ||
     246           2 :         strstr(pszBasename, "_SOURCE_"))
     247           0 :         eDT = GDT_Int16;
     248             :     else
     249           2 :         eDT = GDT_Float32;
     250           2 :     int nWordSize = GDALGetDataTypeSize(eDT) / 8;
     251             : 
     252             :     VSIStatBufL sStat;
     253           2 :     if (strstr(pszBasename, "_5M"))
     254           2 :         sStat.st_size = 180 * 180 * nWordSize;
     255           0 :     else if (strstr(pszBasename, "_30S"))
     256           0 :         sStat.st_size = 1800 * 1800 * nWordSize;
     257           0 :     else if (strstr(pszBasename, "_9S"))
     258           0 :         sStat.st_size = 6000 * 6000 * nWordSize;
     259           0 :     else if (strstr(pszBasename, "_3S"))
     260           0 :         sStat.st_size = 18000 * 18000 * nWordSize;
     261             :     /* Check file size otherwise */
     262           0 :     else if (VSIStatL(poOpenInfo->pszFilename, &sStat) != 0)
     263             :     {
     264           0 :         return nullptr;
     265             :     }
     266             : 
     267           2 :     int nXSize = 0;
     268           2 :     int nYSize = 0;
     269             : 
     270           2 :     double dfPixelSize = 0;
     271           2 :     if (sStat.st_size == 180 * 180 * nWordSize)
     272             :     {
     273             :         /* 5 minute */
     274           2 :         nXSize = 180;
     275           2 :         nYSize = 180;
     276           2 :         dfPixelSize = 5.0 / 60;
     277             :     }
     278           0 :     else if (sStat.st_size == 1800 * 1800 * nWordSize)
     279             :     {
     280             :         /* 30 s */
     281           0 :         nXSize = 1800;
     282           0 :         nYSize = 1800;
     283           0 :         dfPixelSize = 30.0 / 3600;
     284             :     }
     285           0 :     else if (sStat.st_size == 6000 * 6000 * nWordSize)
     286             :     {
     287             :         /* 9 s */
     288           0 :         nXSize = 6000;
     289           0 :         nYSize = 6000;
     290           0 :         dfPixelSize = 9.0 / 3600;
     291             :     }
     292           0 :     else if (sStat.st_size == 18000 * 18000 * nWordSize)
     293             :     {
     294             :         /* 3 s */
     295           0 :         nXSize = 18000;
     296           0 :         nYSize = 18000;
     297           0 :         dfPixelSize = 3.0 / 3600;
     298             :     }
     299             :     else
     300           0 :         return nullptr;
     301             : 
     302             :     /* -------------------------------------------------------------------- */
     303             :     /*      Open file.                                                      */
     304             :     /* -------------------------------------------------------------------- */
     305             : 
     306           4 :     CPLString osFilename = poOpenInfo->pszFilename;
     307           2 :     if ((strstr(poOpenInfo->pszFilename, ".ACE2.gz") ||
     308           2 :          strstr(poOpenInfo->pszFilename, ".ace2.gz")) &&
     309           0 :         !STARTS_WITH(poOpenInfo->pszFilename, "/vsigzip/"))
     310           0 :         osFilename = "/vsigzip/" + osFilename;
     311             : 
     312           2 :     VSILFILE *fpImage = VSIFOpenL(osFilename, "rb+");
     313           2 :     if (fpImage == nullptr)
     314           0 :         return nullptr;
     315             : 
     316             :     /* -------------------------------------------------------------------- */
     317             :     /*      Create the dataset.                                             */
     318             :     /* -------------------------------------------------------------------- */
     319           4 :     auto poDS = std::make_unique<ACE2Dataset>();
     320             : 
     321           2 :     poDS->nRasterXSize = nXSize;
     322           2 :     poDS->nRasterYSize = nYSize;
     323             : 
     324           2 :     poDS->adfGeoTransform[0] = southWestLon;
     325           2 :     poDS->adfGeoTransform[1] = dfPixelSize;
     326           2 :     poDS->adfGeoTransform[2] = 0.0;
     327           2 :     poDS->adfGeoTransform[3] = southWestLat + nYSize * dfPixelSize;
     328           2 :     poDS->adfGeoTransform[4] = 0.0;
     329           2 :     poDS->adfGeoTransform[5] = -dfPixelSize;
     330             : 
     331             :     /* -------------------------------------------------------------------- */
     332             :     /*      Create band information objects                                 */
     333             :     /* -------------------------------------------------------------------- */
     334             :     auto poBand =
     335           4 :         std::make_unique<ACE2RasterBand>(fpImage, eDT, nXSize, nYSize);
     336           2 :     if (!poBand->IsValid())
     337           0 :         return nullptr;
     338           2 :     poDS->SetBand(1, std::move(poBand));
     339             : 
     340             :     /* -------------------------------------------------------------------- */
     341             :     /*      Initialize any PAM information.                                 */
     342             :     /* -------------------------------------------------------------------- */
     343           2 :     poDS->SetDescription(poOpenInfo->pszFilename);
     344           2 :     poDS->TryLoadXML();
     345             : 
     346             :     /* -------------------------------------------------------------------- */
     347             :     /*      Check for overviews.                                            */
     348             :     /* -------------------------------------------------------------------- */
     349           2 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     350             : 
     351           2 :     return poDS.release();
     352             : }
     353             : 
     354             : /************************************************************************/
     355             : /*                          GDALRegister_ACE2()                         */
     356             : /************************************************************************/
     357             : 
     358        1682 : void GDALRegister_ACE2()
     359             : 
     360             : {
     361        1682 :     if (GDALGetDriverByName("ACE2") != nullptr)
     362         301 :         return;
     363             : 
     364        1381 :     GDALDriver *poDriver = new GDALDriver();
     365             : 
     366        1381 :     poDriver->SetDescription("ACE2");
     367        1381 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     368        1381 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "ACE2");
     369        1381 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ace2.html");
     370        1381 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "ACE2");
     371        1381 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     372             : 
     373        1381 :     poDriver->pfnOpen = ACE2Dataset::Open;
     374        1381 :     poDriver->pfnIdentify = ACE2Dataset::Identify;
     375             : 
     376        1381 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     377             : }