/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 | } |