Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  NGSGEOID driver
       4             :  * Purpose:  GDALDataset driver for NGSGEOID dataset.
       5             :  * Author:   Even Rouault, <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2011, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "cpl_string.h"
      14             : #include "cpl_vsi_virtual.h"
      15             : #include "gdal_frmts.h"
      16             : #include "gdal_pam.h"
      17             : #include "ogr_srs_api.h"
      18             : 
      19             : #define HEADER_SIZE (4 * 8 + 3 * 4)
      20             : 
      21             : /************************************************************************/
      22             : /* ==================================================================== */
      23             : /*                            NGSGEOIDDataset                           */
      24             : /* ==================================================================== */
      25             : /************************************************************************/
      26             : 
      27             : class NGSGEOIDRasterBand;
      28             : 
      29             : class NGSGEOIDDataset final : public GDALPamDataset
      30             : {
      31             :     friend class NGSGEOIDRasterBand;
      32             : 
      33             :     VSILFILE *fp;
      34             :     double adfGeoTransform[6];
      35             :     int bIsLittleEndian;
      36             :     mutable OGRSpatialReference m_oSRS{};
      37             : 
      38             :     static int GetHeaderInfo(const GByte *pBuffer, double *padfGeoTransform,
      39             :                              int *pnRows, int *pnCols, int *pbIsLittleEndian);
      40             : 
      41             :   public:
      42             :     NGSGEOIDDataset();
      43             :     virtual ~NGSGEOIDDataset();
      44             : 
      45             :     virtual CPLErr GetGeoTransform(double *) override;
      46             :     const OGRSpatialReference *GetSpatialRef() const override;
      47             : 
      48             :     static GDALDataset *Open(GDALOpenInfo *);
      49             :     static int Identify(GDALOpenInfo *);
      50             : };
      51             : 
      52             : /************************************************************************/
      53             : /* ==================================================================== */
      54             : /*                          NGSGEOIDRasterBand                          */
      55             : /* ==================================================================== */
      56             : /************************************************************************/
      57             : 
      58             : class NGSGEOIDRasterBand final : public GDALPamRasterBand
      59             : {
      60             :     friend class NGSGEOIDDataset;
      61             : 
      62             :   public:
      63             :     explicit NGSGEOIDRasterBand(NGSGEOIDDataset *);
      64             : 
      65             :     virtual CPLErr IReadBlock(int, int, void *) override;
      66             : 
      67           0 :     virtual const char *GetUnitType() override
      68             :     {
      69           0 :         return "m";
      70             :     }
      71             : };
      72             : 
      73             : /************************************************************************/
      74             : /*                        NGSGEOIDRasterBand()                          */
      75             : /************************************************************************/
      76             : 
      77           4 : NGSGEOIDRasterBand::NGSGEOIDRasterBand(NGSGEOIDDataset *poDSIn)
      78             : 
      79             : {
      80           4 :     poDS = poDSIn;
      81           4 :     nBand = 1;
      82             : 
      83           4 :     eDataType = GDT_Float32;
      84             : 
      85           4 :     nBlockXSize = poDS->GetRasterXSize();
      86           4 :     nBlockYSize = 1;
      87           4 : }
      88             : 
      89             : /************************************************************************/
      90             : /*                             IReadBlock()                             */
      91             : /************************************************************************/
      92             : 
      93           2 : CPLErr NGSGEOIDRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
      94             :                                       void *pImage)
      95             : 
      96             : {
      97           2 :     NGSGEOIDDataset *poGDS = reinterpret_cast<NGSGEOIDDataset *>(poDS);
      98             : 
      99             :     /* First values in the file corresponds to the south-most line of the
     100             :      * imagery */
     101           2 :     VSIFSeekL(poGDS->fp,
     102           2 :               HEADER_SIZE +
     103           2 :                   static_cast<vsi_l_offset>(nRasterYSize - 1 - nBlockYOff) *
     104           2 :                       nRasterXSize * 4,
     105             :               SEEK_SET);
     106             : 
     107           2 :     if (static_cast<int>(VSIFReadL(pImage, 4, nRasterXSize, poGDS->fp)) !=
     108           2 :         nRasterXSize)
     109           0 :         return CE_Failure;
     110             : 
     111             : #ifdef CPL_MSB
     112             :     if (poGDS->bIsLittleEndian)
     113             :     {
     114             :         GDALSwapWords(pImage, 4, nRasterXSize, 4);
     115             :     }
     116             : #endif
     117             : 
     118             : #ifdef CPL_LSB
     119           2 :     if (!poGDS->bIsLittleEndian)
     120             :     {
     121           1 :         GDALSwapWords(pImage, 4, nRasterXSize, 4);
     122             :     }
     123             : #endif
     124             : 
     125           2 :     return CE_None;
     126             : }
     127             : 
     128             : /************************************************************************/
     129             : /*                          ~NGSGEOIDDataset()                          */
     130             : /************************************************************************/
     131             : 
     132           4 : NGSGEOIDDataset::NGSGEOIDDataset() : fp(nullptr), bIsLittleEndian(TRUE)
     133             : {
     134           4 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     135           4 :     adfGeoTransform[0] = 0;
     136           4 :     adfGeoTransform[1] = 1;
     137           4 :     adfGeoTransform[2] = 0;
     138           4 :     adfGeoTransform[3] = 0;
     139           4 :     adfGeoTransform[4] = 0;
     140           4 :     adfGeoTransform[5] = 1;
     141           4 : }
     142             : 
     143             : /************************************************************************/
     144             : /*                           ~NGSGEOIDDataset()                         */
     145             : /************************************************************************/
     146             : 
     147           8 : NGSGEOIDDataset::~NGSGEOIDDataset()
     148             : 
     149             : {
     150           4 :     FlushCache(true);
     151           4 :     if (fp)
     152           4 :         VSIFCloseL(fp);
     153           8 : }
     154             : 
     155             : /************************************************************************/
     156             : /*                            GetHeaderInfo()                           */
     157             : /************************************************************************/
     158             : 
     159        2861 : int NGSGEOIDDataset::GetHeaderInfo(const GByte *pBuffer,
     160             :                                    double *padfGeoTransform, int *pnRows,
     161             :                                    int *pnCols, int *pbIsLittleEndian)
     162             : {
     163             :     /* First check IKIND marker to determine if the file */
     164             :     /* is in little or big-endian order, and if it is a valid */
     165             :     /* NGSGEOID dataset */
     166             :     int nIKIND;
     167        2861 :     memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
     168        2861 :     CPL_LSBPTR32(&nIKIND);
     169        2861 :     if (nIKIND == 1)
     170             :     {
     171          27 :         *pbIsLittleEndian = TRUE;
     172             :     }
     173             :     else
     174             :     {
     175        2834 :         memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
     176        2834 :         CPL_MSBPTR32(&nIKIND);
     177        2834 :         if (nIKIND == 1)
     178             :         {
     179           7 :             *pbIsLittleEndian = FALSE;
     180             :         }
     181             :         else
     182             :         {
     183        2827 :             return FALSE;
     184             :         }
     185             :     }
     186             : 
     187             :     double dfSLAT;
     188          34 :     memcpy(&dfSLAT, pBuffer, 8);
     189          34 :     if (*pbIsLittleEndian)
     190             :     {
     191          27 :         CPL_LSBPTR64(&dfSLAT);
     192             :     }
     193             :     else
     194             :     {
     195           7 :         CPL_MSBPTR64(&dfSLAT);
     196             :     }
     197          34 :     pBuffer += 8;
     198             : 
     199             :     double dfWLON;
     200          34 :     memcpy(&dfWLON, pBuffer, 8);
     201          34 :     if (*pbIsLittleEndian)
     202             :     {
     203          27 :         CPL_LSBPTR64(&dfWLON);
     204             :     }
     205             :     else
     206             :     {
     207           7 :         CPL_MSBPTR64(&dfWLON);
     208             :     }
     209          34 :     pBuffer += 8;
     210             : 
     211             :     double dfDLAT;
     212          34 :     memcpy(&dfDLAT, pBuffer, 8);
     213          34 :     if (*pbIsLittleEndian)
     214             :     {
     215          27 :         CPL_LSBPTR64(&dfDLAT);
     216             :     }
     217             :     else
     218             :     {
     219           7 :         CPL_MSBPTR64(&dfDLAT);
     220             :     }
     221          34 :     pBuffer += 8;
     222             : 
     223             :     double dfDLON;
     224          34 :     memcpy(&dfDLON, pBuffer, 8);
     225          34 :     if (*pbIsLittleEndian)
     226             :     {
     227          27 :         CPL_LSBPTR64(&dfDLON);
     228             :     }
     229             :     else
     230             :     {
     231           7 :         CPL_MSBPTR64(&dfDLON);
     232             :     }
     233          34 :     pBuffer += 8;
     234             : 
     235             :     int nNLAT;
     236          34 :     memcpy(&nNLAT, pBuffer, 4);
     237          34 :     if (*pbIsLittleEndian)
     238             :     {
     239          27 :         CPL_LSBPTR32(&nNLAT);
     240             :     }
     241             :     else
     242             :     {
     243           7 :         CPL_MSBPTR32(&nNLAT);
     244             :     }
     245          34 :     pBuffer += 4;
     246             : 
     247             :     int nNLON;
     248          34 :     memcpy(&nNLON, pBuffer, 4);
     249          34 :     if (*pbIsLittleEndian)
     250             :     {
     251          27 :         CPL_LSBPTR32(&nNLON);
     252             :     }
     253             :     else
     254             :     {
     255           7 :         CPL_MSBPTR32(&nNLON);
     256             :     }
     257             :     /*pBuffer += 4;*/
     258             : 
     259             :     /*CPLDebug("NGSGEOID", "SLAT=%f, WLON=%f, DLAT=%f, DLON=%f, NLAT=%d,
     260             :        NLON=%d, IKIND=%d", dfSLAT, dfWLON, dfDLAT, dfDLON, nNLAT, nNLON,
     261             :        nIKIND);*/
     262             : 
     263          34 :     if (nNLAT <= 0 || nNLON <= 0 || dfDLAT <= 1e-15 || dfDLON <= 1e-15)
     264          22 :         return FALSE;
     265             : 
     266             :     /* Grids go over +180 in longitude */
     267             :     // Test written that way to be robust to NaN values
     268          12 :     if (!(dfSLAT >= -90.0 && dfSLAT + nNLAT * dfDLAT <= 90.0 &&
     269          12 :           dfWLON >= -180.0 && dfWLON + nNLON * dfDLON <= 360.0))
     270           0 :         return FALSE;
     271             : 
     272          12 :     padfGeoTransform[0] = dfWLON - dfDLON / 2;
     273          12 :     padfGeoTransform[1] = dfDLON;
     274          12 :     padfGeoTransform[2] = 0.0;
     275          12 :     padfGeoTransform[3] = dfSLAT + nNLAT * dfDLAT - dfDLAT / 2;
     276          12 :     padfGeoTransform[4] = 0.0;
     277          12 :     padfGeoTransform[5] = -dfDLAT;
     278             : 
     279          12 :     *pnRows = nNLAT;
     280          12 :     *pnCols = nNLON;
     281             : 
     282          12 :     return TRUE;
     283             : }
     284             : 
     285             : /************************************************************************/
     286             : /*                             Identify()                               */
     287             : /************************************************************************/
     288             : 
     289       51458 : int NGSGEOIDDataset::Identify(GDALOpenInfo *poOpenInfo)
     290             : {
     291       51458 :     if (poOpenInfo->nHeaderBytes < HEADER_SIZE)
     292       48601 :         return FALSE;
     293             : 
     294             :     double adfGeoTransform[6];
     295             :     int nRows, nCols;
     296             :     int bIsLittleEndian;
     297        2857 :     if (!GetHeaderInfo(poOpenInfo->pabyHeader, adfGeoTransform, &nRows, &nCols,
     298             :                        &bIsLittleEndian))
     299        2849 :         return FALSE;
     300             : 
     301           8 :     return TRUE;
     302             : }
     303             : 
     304             : /************************************************************************/
     305             : /*                                Open()                                */
     306             : /************************************************************************/
     307             : 
     308           4 : GDALDataset *NGSGEOIDDataset::Open(GDALOpenInfo *poOpenInfo)
     309             : 
     310             : {
     311           4 :     if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
     312           0 :         return nullptr;
     313             : 
     314           4 :     if (poOpenInfo->eAccess == GA_Update)
     315             :     {
     316           0 :         ReportUpdateNotSupportedByDriver("NGSGEOID");
     317           0 :         return nullptr;
     318             :     }
     319             : 
     320             :     /* -------------------------------------------------------------------- */
     321             :     /*      Create a corresponding GDALDataset.                             */
     322             :     /* -------------------------------------------------------------------- */
     323           4 :     NGSGEOIDDataset *poDS = new NGSGEOIDDataset();
     324           4 :     poDS->fp = poOpenInfo->fpL;
     325           4 :     poOpenInfo->fpL = nullptr;
     326             : 
     327           4 :     int nRows = 0, nCols = 0;
     328           4 :     GetHeaderInfo(poOpenInfo->pabyHeader, poDS->adfGeoTransform, &nRows, &nCols,
     329             :                   &poDS->bIsLittleEndian);
     330           4 :     poDS->nRasterXSize = nCols;
     331           4 :     poDS->nRasterYSize = nRows;
     332             : 
     333             :     /* -------------------------------------------------------------------- */
     334             :     /*      Create band information objects.                                */
     335             :     /* -------------------------------------------------------------------- */
     336           4 :     poDS->nBands = 1;
     337           4 :     poDS->SetBand(1, new NGSGEOIDRasterBand(poDS));
     338             : 
     339             :     /* -------------------------------------------------------------------- */
     340             :     /*      Initialize any PAM information.                                 */
     341             :     /* -------------------------------------------------------------------- */
     342           4 :     poDS->SetDescription(poOpenInfo->pszFilename);
     343           4 :     poDS->TryLoadXML();
     344             : 
     345             :     /* -------------------------------------------------------------------- */
     346             :     /*      Support overviews.                                              */
     347             :     /* -------------------------------------------------------------------- */
     348           4 :     poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
     349           4 :     return poDS;
     350             : }
     351             : 
     352             : /************************************************************************/
     353             : /*                          GetGeoTransform()                           */
     354             : /************************************************************************/
     355             : 
     356           2 : CPLErr NGSGEOIDDataset::GetGeoTransform(double *padfTransform)
     357             : 
     358             : {
     359           2 :     memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));
     360             : 
     361           2 :     return CE_None;
     362             : }
     363             : 
     364             : /************************************************************************/
     365             : /*                         GetSpatialRef()                              */
     366             : /************************************************************************/
     367             : 
     368           2 : const OGRSpatialReference *NGSGEOIDDataset::GetSpatialRef() const
     369             : {
     370           2 :     if (!m_oSRS.IsEmpty())
     371             :     {
     372           0 :         return &m_oSRS;
     373             :     }
     374             : 
     375             :     const CPLString osFilename =
     376           6 :         CPLString(CPLGetBasenameSafe(GetDescription())).tolower();
     377             : 
     378             :     // See https://www.ngs.noaa.gov/GEOID/GEOID12B/faq_2012B.shtml
     379             : 
     380             :     // GEOID2012 files ?
     381           2 :     if (STARTS_WITH(osFilename, "g2012") && osFilename.size() >= 7)
     382             :     {
     383           0 :         if (osFilename[6] == 'h' /* Hawai */ ||
     384           0 :             osFilename[6] == 's' /* Samoa */)
     385             :         {
     386             :             // NAD83 (PA11)
     387           0 :             m_oSRS.importFromEPSG(6322);
     388             :         }
     389           0 :         else if (osFilename[6] == 'g' /* Guam */)
     390             :         {
     391             :             // NAD83 (MA11)
     392           0 :             m_oSRS.importFromEPSG(6325);
     393             :         }
     394             :         else
     395             :         {
     396             :             // NAD83 (2011)
     397           0 :             m_oSRS.importFromEPSG(6318);
     398             :         }
     399             : 
     400           0 :         return &m_oSRS;
     401             :     }
     402             : 
     403             :     // USGG2012 files ? We should return IGS08, but there is only a
     404             :     // geocentric CRS in EPSG, so manually forge a geographic one from it
     405           2 :     if (STARTS_WITH(osFilename, "s2012"))
     406             :     {
     407           0 :         m_oSRS.importFromWkt(
     408             :             "GEOGCS[\"IGS08\",\n"
     409             :             "    DATUM[\"IGS08\",\n"
     410             :             "        SPHEROID[\"GRS 1980\",6378137,298.257222101,\n"
     411             :             "            AUTHORITY[\"EPSG\",\"7019\"]],\n"
     412             :             "        AUTHORITY[\"EPSG\",\"1141\"]],\n"
     413             :             "    PRIMEM[\"Greenwich\",0,\n"
     414             :             "        AUTHORITY[\"EPSG\",\"8901\"]],\n"
     415             :             "    UNIT[\"degree\",0.0174532925199433,\n"
     416             :             "        AUTHORITY[\"EPSG\",\"9122\"]]]");
     417           0 :         return &m_oSRS;
     418             :     }
     419             : 
     420           2 :     m_oSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
     421           2 :     return &m_oSRS;
     422             : }
     423             : 
     424             : /************************************************************************/
     425             : /*                       GDALRegister_NGSGEOID()                        */
     426             : /************************************************************************/
     427             : 
     428        1686 : void GDALRegister_NGSGEOID()
     429             : 
     430             : {
     431        1686 :     if (GDALGetDriverByName("NGSGEOID") != nullptr)
     432         302 :         return;
     433             : 
     434        1384 :     GDALDriver *poDriver = new GDALDriver();
     435             : 
     436        1384 :     poDriver->SetDescription("NGSGEOID");
     437        1384 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     438        1384 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NOAA NGS Geoid Height Grids");
     439        1384 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
     440        1384 :                               "drivers/raster/ngsgeoid.html");
     441        1384 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "bin");
     442             : 
     443        1384 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     444             : 
     445        1384 :     poDriver->pfnOpen = NGSGEOIDDataset::Open;
     446        1384 :     poDriver->pfnIdentify = NGSGEOIDDataset::Identify;
     447             : 
     448        1384 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     449             : }