Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  Implementation of NOAA .b format used for GEOCON / NADCON5 grids
       5             :  * Author:   Even Rouault <even dot rouault at spatialys.com>
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2022, Even Rouault <even dot rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "cpl_conv.h"
      14             : #include "cpl_string.h"
      15             : #include "gdal_frmts.h"
      16             : #include "rawdataset.h"
      17             : #include "ogr_srs_api.h"
      18             : 
      19             : #include <limits>
      20             : 
      21             : // Specification of the format is at "paragraph 10.2 ".b" grids (GEOCON and
      22             : // NADCON 5.0)" of "NOAA Technical Report NOS NGS 63" at
      23             : // https://geodesy.noaa.gov/library/pdfs/NOAA_TR_NOS_NGS_0063.pdf
      24             : 
      25             : constexpr int HEADER_SIZE = 52;
      26             : constexpr int FORTRAN_HEADER_SIZE = 4;
      27             : constexpr int FORTRAN_TRAILER_SIZE = 4;
      28             : 
      29             : /************************************************************************/
      30             : /* ==================================================================== */
      31             : /*                          NOAA_B_Dataset                              */
      32             : /* ==================================================================== */
      33             : /************************************************************************/
      34             : 
      35             : class NOAA_B_Dataset final : public RawDataset
      36             : {
      37             :     OGRSpatialReference m_oSRS{};
      38             :     double m_adfGeoTransform[6];
      39             : 
      40             :     CPL_DISALLOW_COPY_ASSIGN(NOAA_B_Dataset)
      41             : 
      42             :     static int IdentifyEx(GDALOpenInfo *poOpenInfo, bool &bBigEndianOut);
      43             : 
      44           4 :     CPLErr Close() override
      45             :     {
      46           4 :         return GDALPamDataset::Close();
      47             :     }
      48             : 
      49             :   public:
      50           4 :     NOAA_B_Dataset()
      51           4 :     {
      52           4 :         m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
      53             : 
      54           4 :         m_adfGeoTransform[0] = 0.0;
      55           4 :         m_adfGeoTransform[1] = 1.0;
      56           4 :         m_adfGeoTransform[2] = 0.0;
      57           4 :         m_adfGeoTransform[3] = 0.0;
      58           4 :         m_adfGeoTransform[4] = 0.0;
      59           4 :         m_adfGeoTransform[5] = 1.0;
      60           4 :     }
      61             : 
      62             :     CPLErr GetGeoTransform(double *padfTransform) override;
      63             : 
      64           0 :     const OGRSpatialReference *GetSpatialRef() const override
      65             :     {
      66           0 :         return &m_oSRS;
      67             :     }
      68             : 
      69             :     static GDALDataset *Open(GDALOpenInfo *);
      70             :     static int Identify(GDALOpenInfo *);
      71             : };
      72             : 
      73             : /************************************************************************/
      74             : /* ==================================================================== */
      75             : /*                          NOAA_B_Dataset                              */
      76             : /* ==================================================================== */
      77             : /************************************************************************/
      78             : 
      79             : /************************************************************************/
      80             : /*                           GetHeaderValues()                          */
      81             : /************************************************************************/
      82             : 
      83          16 : static void GetHeaderValues(const GDALOpenInfo *poOpenInfo, double &dfSWLat,
      84             :                             double &dfSWLon, double &dfDeltaLat,
      85             :                             double &dfDeltaLon, int32_t &nRows, int32_t &nCols,
      86             :                             int32_t &iKind, bool bBigEndian)
      87             : {
      88          64 :     const auto ReadFloat64 = [bBigEndian](const GByte *&ptr)
      89             :     {
      90             :         double v;
      91          64 :         memcpy(&v, ptr, sizeof(v));
      92          64 :         ptr += sizeof(v);
      93          64 :         if (bBigEndian)
      94          40 :             CPL_MSBPTR64(&v);
      95             :         else
      96          24 :             CPL_LSBPTR64(&v);
      97          64 :         return v;
      98          16 :     };
      99             : 
     100          48 :     const auto ReadInt32 = [bBigEndian](const GByte *&ptr)
     101             :     {
     102             :         int32_t v;
     103          48 :         memcpy(&v, ptr, sizeof(v));
     104          48 :         ptr += sizeof(v);
     105          48 :         if (bBigEndian)
     106          30 :             CPL_MSBPTR32(&v);
     107             :         else
     108          18 :             CPL_LSBPTR32(&v);
     109          48 :         return v;
     110          16 :     };
     111             : 
     112          16 :     const GByte *ptr = poOpenInfo->pabyHeader + FORTRAN_HEADER_SIZE;
     113             : 
     114          16 :     dfSWLat = ReadFloat64(ptr);
     115          16 :     dfSWLon = ReadFloat64(ptr);
     116          16 :     dfDeltaLat = ReadFloat64(ptr);
     117          16 :     dfDeltaLon = ReadFloat64(ptr);
     118             : 
     119          16 :     nRows = ReadInt32(ptr);
     120          16 :     nCols = ReadInt32(ptr);
     121          16 :     iKind = ReadInt32(ptr);
     122          16 : }
     123             : 
     124             : /************************************************************************/
     125             : /*                              Identify()                              */
     126             : /************************************************************************/
     127             : 
     128       51266 : int NOAA_B_Dataset::IdentifyEx(GDALOpenInfo *poOpenInfo, bool &bBigEndianOut)
     129             : 
     130             : {
     131       51266 :     if (poOpenInfo->nHeaderBytes < HEADER_SIZE)
     132       48020 :         return FALSE;
     133             : 
     134             : #if !defined(FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION)
     135        3246 :     if (!EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "b"))
     136        3239 :         return FALSE;
     137             : #endif
     138             : 
     139             :     // Sanity checks on header
     140             :     double dfSWLat;
     141             :     double dfSWLon;
     142             :     double dfDeltaLat;
     143             :     double dfDeltaLon;
     144             :     int32_t nRows;
     145             :     int32_t nCols;
     146             :     int32_t iKind;
     147             : 
     148             :     // Fun... nadcon5 files are encoded in big-endian, but vertcon3 files...
     149             :     // in little-endian. We could probably figure that out directly from the
     150             :     // 4 bytes which are 0x00 0x00 0x00 0x2C for nadcon5, and the reverse for
     151             :     // vertcon3, but the semantics of those 4 bytes is undocumented.
     152             :     // So try both possibilities and rely on sanity checks.
     153          12 :     for (int i = 0; i < 2; ++i)
     154             :     {
     155          12 :         const bool bBigEndian = i == 0 ? true : false;
     156          12 :         GetHeaderValues(poOpenInfo, dfSWLat, dfSWLon, dfDeltaLat, dfDeltaLon,
     157             :                         nRows, nCols, iKind, bBigEndian);
     158          12 :         if (!(fabs(dfSWLat) <= 90))
     159           0 :             continue;
     160          12 :         if (!(fabs(dfSWLon) <=
     161             :               360))  // NADCON5 grids typically have SWLon > 180
     162           0 :             continue;
     163          12 :         if (!(dfDeltaLat > 0 && dfDeltaLat <= 1))
     164           0 :             continue;
     165          12 :         if (!(dfDeltaLon > 0 && dfDeltaLon <= 1))
     166           0 :             continue;
     167          12 :         if (!(nRows > 0 && dfSWLat + (nRows - 1) * dfDeltaLat <= 90))
     168           0 :             continue;
     169          12 :         if (!(nCols > 0 && (nCols - 1) * dfDeltaLon <= 360))
     170           0 :             continue;
     171          12 :         if (!(iKind >= -1 && iKind <= 2))
     172           4 :             continue;
     173             : 
     174           8 :         bBigEndianOut = bBigEndian;
     175           8 :         return TRUE;
     176             :     }
     177           0 :     return FALSE;
     178             : }
     179             : 
     180       51263 : int NOAA_B_Dataset::Identify(GDALOpenInfo *poOpenInfo)
     181             : 
     182             : {
     183       51263 :     bool bBigEndian = false;
     184      102521 :     return IdentifyEx(poOpenInfo, bBigEndian);
     185             : }
     186             : 
     187             : /************************************************************************/
     188             : /*                          GetGeoTransform()                           */
     189             : /************************************************************************/
     190             : 
     191           2 : CPLErr NOAA_B_Dataset::GetGeoTransform(double *padfTransform)
     192             : 
     193             : {
     194           2 :     memcpy(padfTransform, m_adfGeoTransform, sizeof(double) * 6);
     195           2 :     return CE_None;
     196             : }
     197             : 
     198             : /************************************************************************/
     199             : /*                                Open()                                */
     200             : /************************************************************************/
     201             : 
     202           4 : GDALDataset *NOAA_B_Dataset::Open(GDALOpenInfo *poOpenInfo)
     203             : 
     204             : {
     205           4 :     bool bBigEndian = false;
     206           8 :     if (!IdentifyEx(poOpenInfo, bBigEndian) || poOpenInfo->fpL == nullptr ||
     207           4 :         poOpenInfo->eAccess == GA_Update)
     208             :     {
     209           0 :         return nullptr;
     210             :     }
     211             : 
     212             :     /* -------------------------------------------------------------------- */
     213             :     /*      Read the header.                                                */
     214             :     /* -------------------------------------------------------------------- */
     215             :     double dfSWLat;
     216             :     double dfSWLon;
     217             :     double dfDeltaLat;
     218             :     double dfDeltaLon;
     219             :     int32_t nRows;
     220             :     int32_t nCols;
     221             :     int32_t iKind;
     222           4 :     GetHeaderValues(poOpenInfo, dfSWLat, dfSWLon, dfDeltaLat, dfDeltaLon, nRows,
     223             :                     nCols, iKind, bBigEndian);
     224             : 
     225           4 :     if (iKind == -1)
     226             :     {
     227           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     228             :                  "KIND = -1 in NOAA .b dataset not supported");
     229           0 :         return nullptr;
     230             :     }
     231             : 
     232           4 :     const GDALDataType eDT =
     233             :         // iKind == -1 ? GDT_Int16 :
     234           8 :         iKind == 0   ? GDT_Int32
     235           4 :         : iKind == 1 ? GDT_Float32
     236             :                      : GDT_Int16;
     237           4 :     const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
     238           8 :     if (!GDALCheckDatasetDimensions(nCols, nRows) ||
     239           8 :         (nDTSize > 0 && static_cast<vsi_l_offset>(nCols) * nRows >
     240           4 :                             std::numeric_limits<vsi_l_offset>::max() / nDTSize))
     241             :     {
     242           0 :         return nullptr;
     243             :     }
     244           8 :     if (nDTSize > 0 && nCols > (std::numeric_limits<int>::max() -
     245           4 :                                 FORTRAN_HEADER_SIZE - FORTRAN_TRAILER_SIZE) /
     246             :                                    nDTSize)
     247             :     {
     248           0 :         return nullptr;
     249             :     }
     250           4 :     const int nLineSize =
     251           4 :         FORTRAN_HEADER_SIZE + nCols * nDTSize + FORTRAN_TRAILER_SIZE;
     252             : 
     253             :     /* -------------------------------------------------------------------- */
     254             :     /*      Create a corresponding GDALDataset.                             */
     255             :     /* -------------------------------------------------------------------- */
     256           8 :     auto poDS = std::make_unique<NOAA_B_Dataset>();
     257             : 
     258           4 :     poDS->nRasterXSize = nCols;
     259           4 :     poDS->nRasterYSize = nRows;
     260             : 
     261             :     // Adjust longitude > 180 to [-180, 180] range
     262           4 :     if (dfSWLon > 180)
     263           0 :         dfSWLon -= 360;
     264             : 
     265             :     // Convert from south-west center-of-pixel convention to
     266             :     // north-east pixel-corner convention
     267           4 :     poDS->m_adfGeoTransform[0] = dfSWLon - dfDeltaLon / 2;
     268           4 :     poDS->m_adfGeoTransform[1] = dfDeltaLon;
     269           4 :     poDS->m_adfGeoTransform[2] = 0.0;
     270           8 :     poDS->m_adfGeoTransform[3] =
     271           4 :         dfSWLat + (nRows - 1) * dfDeltaLat + dfDeltaLat / 2;
     272           4 :     poDS->m_adfGeoTransform[4] = 0.0;
     273           4 :     poDS->m_adfGeoTransform[5] = -dfDeltaLat;
     274             : 
     275             :     /* -------------------------------------------------------------------- */
     276             :     /*      Create band information object.                                 */
     277             :     /* -------------------------------------------------------------------- */
     278             : 
     279             :     // Borrow file handle
     280           4 :     VSILFILE *fpImage = poOpenInfo->fpL;
     281           4 :     poOpenInfo->fpL = nullptr;
     282             : 
     283             :     // Records are presented from the southern-most to the northern-most
     284             :     auto poBand = RawRasterBand::Create(
     285           4 :         poDS.get(), 1, fpImage,
     286             :         // skip to beginning of northern-most line
     287             :         HEADER_SIZE +
     288           8 :             static_cast<vsi_l_offset>(poDS->nRasterYSize - 1) * nLineSize +
     289             :             FORTRAN_HEADER_SIZE,
     290             :         nDTSize, -nLineSize, eDT,
     291             :         bBigEndian ? RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN
     292             :                    : RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
     293           8 :         RawRasterBand::OwnFP::YES);
     294           4 :     if (!poBand)
     295           0 :         return nullptr;
     296           4 :     poDS->SetBand(1, std::move(poBand));
     297             : 
     298             :     /* -------------------------------------------------------------------- */
     299             :     /*      Guess CRS from filename.                                        */
     300             :     /* -------------------------------------------------------------------- */
     301           8 :     const std::string osFilename(CPLGetFilename(poOpenInfo->pszFilename));
     302             : 
     303             :     static const struct
     304             :     {
     305             :         const char *pszPrefix;
     306             :         int nEPSGCode;
     307             :     }
     308             :     // Cf https://geodesy.noaa.gov/pub/nadcon5/20160901release/Builds/
     309             :     asFilenameToCRS[] = {
     310             :         {"nadcon5.nad27.", 4267},       // NAD27
     311             :         {"nadcon5.pr40.", 4139},        // Puerto Rico (1940)
     312             :         {"nadcon5.ohd.", 4135},         // Old Hawaian
     313             :         {"nadcon5.sl1952.", 4136},      // Saint Lawrence Island (1952)
     314             :         {"nadcon5.sp1952.", 4137},      // Saint Paul Island (1952)
     315             :         {"nadcon5.sg1952.", 4138},      // Saint George Island (1952)
     316             :         {"nadcon5.as62.", 4169},        // American Samoa 1962
     317             :         {"nadcon5.gu63.", 4675},        // Guam 1963
     318             :         {"nadcon5.nad83_1986.", 4269},  // NAD83
     319             :         {"nadcon5.nad83_harn.", 4152},  // NAD83(HARN)
     320             :         {"nadcon5.nad83_1992.",
     321             :          4152},  // NAD83(1992) for Alaska is NAD83(HARN) in EPSG
     322             :         {"nadcon5.nad83_1993.",
     323             :          4152},  // NAD83(1993) for American Samoa, PRVI, Guam and Hawaii is
     324             :                  // NAD83(HARN) in EPSG
     325             :         {"nadcon5.nad83_1997.", 8545},  // NAD83(HARN Corrected)
     326             :         {"nadcon5.nad83_fbn.", 8860},   // NAD83(FBN)
     327             :         {"nadcon5.nad83_2002.",
     328             :          8860},  // NAD83(2002) for Alaska, PRVI and Guam is NAD83(FBN) in EPSG
     329             :         {"nadcon5.nad83_2007.", 4759},  // NAD83(NSRS2007)
     330             :     };
     331             : 
     332          68 :     for (const auto &sPair : asFilenameToCRS)
     333             :     {
     334          64 :         if (STARTS_WITH_CI(osFilename.c_str(), sPair.pszPrefix))
     335             :         {
     336           0 :             poDS->m_oSRS.importFromEPSG(sPair.nEPSGCode);
     337           0 :             break;
     338             :         }
     339             :     }
     340           4 :     if (poDS->m_oSRS.IsEmpty())
     341             :     {
     342           4 :         poDS->m_oSRS.importFromWkt(
     343             :             "GEOGCRS[\"Unspecified geographic CRS\",DATUM[\"Unspecified datum "
     344             :             "based on GRS80 ellipsoid\",ELLIPSOID[\"GRS "
     345             :             "1980\",6378137,298.257222101]],CS[ellipsoidal,2],AXIS[\"geodetic "
     346             :             "latitude (Lat)\",north,ANGLEUNIT[\"degree\",0.0174532925199433]], "
     347             :             "       AXIS[\"geodetic longitude "
     348             :             "(Lon)\",east,ORDER[2],ANGLEUNIT[\"degree\",0.0174532925199433]]]");
     349             :     }
     350             : 
     351             :     /* -------------------------------------------------------------------- */
     352             :     /*      Initialize any PAM information.                                 */
     353             :     /* -------------------------------------------------------------------- */
     354           4 :     poDS->SetDescription(poOpenInfo->pszFilename);
     355           4 :     poDS->TryLoadXML();
     356             : 
     357             :     /* -------------------------------------------------------------------- */
     358             :     /*      Check for overviews.                                            */
     359             :     /* -------------------------------------------------------------------- */
     360           4 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     361             : 
     362           4 :     return poDS.release();
     363             : }
     364             : 
     365             : /************************************************************************/
     366             : /*                        GDALRegister_NOAA_B()                         */
     367             : /************************************************************************/
     368             : 
     369        1595 : void GDALRegister_NOAA_B()
     370             : {
     371        1595 :     if (GDALGetDriverByName("NOAA_B") != nullptr)
     372         302 :         return;
     373             : 
     374        1293 :     GDALDriver *poDriver = new GDALDriver();
     375             : 
     376        1293 :     poDriver->SetDescription("NOAA_B");
     377        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     378        1293 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
     379        1293 :                               "NOAA GEOCON/NADCON5 .b format");
     380        1293 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "b");
     381        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     382        1293 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/noaa_b.html");
     383             : 
     384        1293 :     poDriver->pfnIdentify = NOAA_B_Dataset::Identify;
     385        1293 :     poDriver->pfnOpen = NOAA_B_Dataset::Open;
     386             : 
     387        1293 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     388             : }