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