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       51414 : int ACE2Dataset::Identify(GDALOpenInfo *poOpenInfo)
     191             : 
     192             : {
     193      102823 :     if (!(EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "ACE2") ||
     194       51409 :           strstr(poOpenInfo->pszFilename, ".ACE2.gz") ||
     195       51409 :           strstr(poOpenInfo->pszFilename, ".ace2.gz")))
     196       51407 :         return FALSE;
     197             : 
     198           6 :     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           2 :     const char *pszBasename = CPLGetBasename(poOpenInfo->pszFilename);
     212             : 
     213           2 :     if (strlen(pszBasename) < 7)
     214           0 :         return nullptr;
     215             : 
     216             :     /* Determine southwest coordinates from filename */
     217             : 
     218             :     /* e.g. 30S120W_5M.ACE2 */
     219           2 :     char pszLatLonValueString[4] = {'\0'};
     220           2 :     memset(pszLatLonValueString, 0, 4);
     221             :     // cppcheck-suppress redundantCopy
     222           2 :     strncpy(pszLatLonValueString, &pszBasename[0], 2);
     223           2 :     int southWestLat = atoi(pszLatLonValueString);
     224           2 :     memset(pszLatLonValueString, 0, 4);
     225             :     // cppcheck-suppress redundantCopy
     226           2 :     strncpy(pszLatLonValueString, &pszBasename[3], 3);
     227           2 :     int southWestLon = atoi(pszLatLonValueString);
     228             : 
     229           2 :     if (pszBasename[2] == 'N' || pszBasename[2] == 'n')
     230             :         /*southWestLat = southWestLat*/;
     231           0 :     else if (pszBasename[2] == 'S' || pszBasename[2] == 's')
     232           0 :         southWestLat = southWestLat * -1;
     233             :     else
     234           0 :         return nullptr;
     235             : 
     236           2 :     if (pszBasename[6] == 'E' || pszBasename[6] == 'e')
     237             :         /*southWestLon = southWestLon*/;
     238           0 :     else if (pszBasename[6] == 'W' || pszBasename[6] == 'w')
     239           0 :         southWestLon = southWestLon * -1;
     240             :     else
     241           0 :         return nullptr;
     242             : 
     243           2 :     GDALDataType eDT = GDT_Unknown;
     244           2 :     if (strstr(pszBasename, "_CONF_") || strstr(pszBasename, "_QUALITY_") ||
     245           2 :         strstr(pszBasename, "_SOURCE_"))
     246           0 :         eDT = GDT_Int16;
     247             :     else
     248           2 :         eDT = GDT_Float32;
     249           2 :     int nWordSize = GDALGetDataTypeSize(eDT) / 8;
     250             : 
     251             :     VSIStatBufL sStat;
     252           2 :     if (strstr(pszBasename, "_5M"))
     253           2 :         sStat.st_size = 180 * 180 * nWordSize;
     254           0 :     else if (strstr(pszBasename, "_30S"))
     255           0 :         sStat.st_size = 1800 * 1800 * nWordSize;
     256           0 :     else if (strstr(pszBasename, "_9S"))
     257           0 :         sStat.st_size = 6000 * 6000 * nWordSize;
     258           0 :     else if (strstr(pszBasename, "_3S"))
     259           0 :         sStat.st_size = 18000 * 18000 * nWordSize;
     260             :     /* Check file size otherwise */
     261           0 :     else if (VSIStatL(poOpenInfo->pszFilename, &sStat) != 0)
     262             :     {
     263           0 :         return nullptr;
     264             :     }
     265             : 
     266           2 :     int nXSize = 0;
     267           2 :     int nYSize = 0;
     268             : 
     269           2 :     double dfPixelSize = 0;
     270           2 :     if (sStat.st_size == 180 * 180 * nWordSize)
     271             :     {
     272             :         /* 5 minute */
     273           2 :         nXSize = 180;
     274           2 :         nYSize = 180;
     275           2 :         dfPixelSize = 5.0 / 60;
     276             :     }
     277           0 :     else if (sStat.st_size == 1800 * 1800 * nWordSize)
     278             :     {
     279             :         /* 30 s */
     280           0 :         nXSize = 1800;
     281           0 :         nYSize = 1800;
     282           0 :         dfPixelSize = 30.0 / 3600;
     283             :     }
     284           0 :     else if (sStat.st_size == 6000 * 6000 * nWordSize)
     285             :     {
     286             :         /* 9 s */
     287           0 :         nXSize = 6000;
     288           0 :         nYSize = 6000;
     289           0 :         dfPixelSize = 9.0 / 3600;
     290             :     }
     291           0 :     else if (sStat.st_size == 18000 * 18000 * nWordSize)
     292             :     {
     293             :         /* 3 s */
     294           0 :         nXSize = 18000;
     295           0 :         nYSize = 18000;
     296           0 :         dfPixelSize = 3.0 / 3600;
     297             :     }
     298             :     else
     299           0 :         return nullptr;
     300             : 
     301             :     /* -------------------------------------------------------------------- */
     302             :     /*      Open file.                                                      */
     303             :     /* -------------------------------------------------------------------- */
     304             : 
     305           4 :     CPLString osFilename = poOpenInfo->pszFilename;
     306           2 :     if ((strstr(poOpenInfo->pszFilename, ".ACE2.gz") ||
     307           2 :          strstr(poOpenInfo->pszFilename, ".ace2.gz")) &&
     308           0 :         !STARTS_WITH(poOpenInfo->pszFilename, "/vsigzip/"))
     309           0 :         osFilename = "/vsigzip/" + osFilename;
     310             : 
     311           2 :     VSILFILE *fpImage = VSIFOpenL(osFilename, "rb+");
     312           2 :     if (fpImage == nullptr)
     313           0 :         return nullptr;
     314             : 
     315             :     /* -------------------------------------------------------------------- */
     316             :     /*      Create the dataset.                                             */
     317             :     /* -------------------------------------------------------------------- */
     318           4 :     auto poDS = std::make_unique<ACE2Dataset>();
     319             : 
     320           2 :     poDS->nRasterXSize = nXSize;
     321           2 :     poDS->nRasterYSize = nYSize;
     322             : 
     323           2 :     poDS->adfGeoTransform[0] = southWestLon;
     324           2 :     poDS->adfGeoTransform[1] = dfPixelSize;
     325           2 :     poDS->adfGeoTransform[2] = 0.0;
     326           2 :     poDS->adfGeoTransform[3] = southWestLat + nYSize * dfPixelSize;
     327           2 :     poDS->adfGeoTransform[4] = 0.0;
     328           2 :     poDS->adfGeoTransform[5] = -dfPixelSize;
     329             : 
     330             :     /* -------------------------------------------------------------------- */
     331             :     /*      Create band information objects                                 */
     332             :     /* -------------------------------------------------------------------- */
     333             :     auto poBand =
     334           4 :         std::make_unique<ACE2RasterBand>(fpImage, eDT, nXSize, nYSize);
     335           2 :     if (!poBand->IsValid())
     336           0 :         return nullptr;
     337           2 :     poDS->SetBand(1, std::move(poBand));
     338             : 
     339             :     /* -------------------------------------------------------------------- */
     340             :     /*      Initialize any PAM information.                                 */
     341             :     /* -------------------------------------------------------------------- */
     342           2 :     poDS->SetDescription(poOpenInfo->pszFilename);
     343           2 :     poDS->TryLoadXML();
     344             : 
     345             :     /* -------------------------------------------------------------------- */
     346             :     /*      Check for overviews.                                            */
     347             :     /* -------------------------------------------------------------------- */
     348           2 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     349             : 
     350           2 :     return poDS.release();
     351             : }
     352             : 
     353             : /************************************************************************/
     354             : /*                          GDALRegister_ACE2()                         */
     355             : /************************************************************************/
     356             : 
     357        1595 : void GDALRegister_ACE2()
     358             : 
     359             : {
     360        1595 :     if (GDALGetDriverByName("ACE2") != nullptr)
     361         302 :         return;
     362             : 
     363        1293 :     GDALDriver *poDriver = new GDALDriver();
     364             : 
     365        1293 :     poDriver->SetDescription("ACE2");
     366        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     367        1293 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "ACE2");
     368        1293 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ace2.html");
     369        1293 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "ACE2");
     370        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     371             : 
     372        1293 :     poDriver->pfnOpen = ACE2Dataset::Open;
     373        1293 :     poDriver->pfnIdentify = ACE2Dataset::Identify;
     374             : 
     375        1293 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     376             : }