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 >) 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 58862 : int NOAA_B_Dataset::IdentifyEx(GDALOpenInfo *poOpenInfo, bool &bBigEndianOut)
123 :
124 : {
125 58862 : if (poOpenInfo->nHeaderBytes < HEADER_SIZE)
126 55316 : return FALSE;
127 :
128 : #if !defined(FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION)
129 3546 : 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 58859 : int NOAA_B_Dataset::Identify(GDALOpenInfo *poOpenInfo)
175 :
176 : {
177 58859 : bool bBigEndian = false;
178 58859 : return IdentifyEx(poOpenInfo, bBigEndian);
179 : }
180 :
181 : /************************************************************************/
182 : /* GetGeoTransform() */
183 : /************************************************************************/
184 :
185 2 : CPLErr NOAA_B_Dataset::GetGeoTransform(GDALGeoTransform >) 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 : }
|