Coverage Report

Created: 2025-06-09 07:02

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