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