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