Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: NGSGEOID driver
4 : * Purpose: GDALDataset driver for NGSGEOID dataset.
5 : * Author: Even Rouault, <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2011, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_string.h"
14 : #include "cpl_vsi_virtual.h"
15 : #include "gdal_frmts.h"
16 : #include "gdal_pam.h"
17 : #include "ogr_srs_api.h"
18 :
19 : #define HEADER_SIZE (4 * 8 + 3 * 4)
20 :
21 : /************************************************************************/
22 : /* ==================================================================== */
23 : /* NGSGEOIDDataset */
24 : /* ==================================================================== */
25 : /************************************************************************/
26 :
27 : class NGSGEOIDRasterBand;
28 :
29 : class NGSGEOIDDataset final : public GDALPamDataset
30 : {
31 : friend class NGSGEOIDRasterBand;
32 :
33 : VSILFILE *fp;
34 : double adfGeoTransform[6];
35 : int bIsLittleEndian;
36 : mutable OGRSpatialReference m_oSRS{};
37 :
38 : static int GetHeaderInfo(const GByte *pBuffer, double *padfGeoTransform,
39 : int *pnRows, int *pnCols, int *pbIsLittleEndian);
40 :
41 : public:
42 : NGSGEOIDDataset();
43 : virtual ~NGSGEOIDDataset();
44 :
45 : virtual CPLErr GetGeoTransform(double *) override;
46 : const OGRSpatialReference *GetSpatialRef() const override;
47 :
48 : static GDALDataset *Open(GDALOpenInfo *);
49 : static int Identify(GDALOpenInfo *);
50 : };
51 :
52 : /************************************************************************/
53 : /* ==================================================================== */
54 : /* NGSGEOIDRasterBand */
55 : /* ==================================================================== */
56 : /************************************************************************/
57 :
58 : class NGSGEOIDRasterBand final : public GDALPamRasterBand
59 : {
60 : friend class NGSGEOIDDataset;
61 :
62 : public:
63 : explicit NGSGEOIDRasterBand(NGSGEOIDDataset *);
64 :
65 : virtual CPLErr IReadBlock(int, int, void *) override;
66 :
67 0 : virtual const char *GetUnitType() override
68 : {
69 0 : return "m";
70 : }
71 : };
72 :
73 : /************************************************************************/
74 : /* NGSGEOIDRasterBand() */
75 : /************************************************************************/
76 :
77 4 : NGSGEOIDRasterBand::NGSGEOIDRasterBand(NGSGEOIDDataset *poDSIn)
78 :
79 : {
80 4 : poDS = poDSIn;
81 4 : nBand = 1;
82 :
83 4 : eDataType = GDT_Float32;
84 :
85 4 : nBlockXSize = poDS->GetRasterXSize();
86 4 : nBlockYSize = 1;
87 4 : }
88 :
89 : /************************************************************************/
90 : /* IReadBlock() */
91 : /************************************************************************/
92 :
93 2 : CPLErr NGSGEOIDRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
94 : void *pImage)
95 :
96 : {
97 2 : NGSGEOIDDataset *poGDS = reinterpret_cast<NGSGEOIDDataset *>(poDS);
98 :
99 : /* First values in the file corresponds to the south-most line of the
100 : * imagery */
101 2 : VSIFSeekL(poGDS->fp,
102 2 : HEADER_SIZE +
103 2 : static_cast<vsi_l_offset>(nRasterYSize - 1 - nBlockYOff) *
104 2 : nRasterXSize * 4,
105 : SEEK_SET);
106 :
107 2 : if (static_cast<int>(VSIFReadL(pImage, 4, nRasterXSize, poGDS->fp)) !=
108 2 : nRasterXSize)
109 0 : return CE_Failure;
110 :
111 : #ifdef CPL_MSB
112 : if (poGDS->bIsLittleEndian)
113 : {
114 : GDALSwapWords(pImage, 4, nRasterXSize, 4);
115 : }
116 : #endif
117 :
118 : #ifdef CPL_LSB
119 2 : if (!poGDS->bIsLittleEndian)
120 : {
121 1 : GDALSwapWords(pImage, 4, nRasterXSize, 4);
122 : }
123 : #endif
124 :
125 2 : return CE_None;
126 : }
127 :
128 : /************************************************************************/
129 : /* ~NGSGEOIDDataset() */
130 : /************************************************************************/
131 :
132 4 : NGSGEOIDDataset::NGSGEOIDDataset() : fp(nullptr), bIsLittleEndian(TRUE)
133 : {
134 4 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
135 4 : adfGeoTransform[0] = 0;
136 4 : adfGeoTransform[1] = 1;
137 4 : adfGeoTransform[2] = 0;
138 4 : adfGeoTransform[3] = 0;
139 4 : adfGeoTransform[4] = 0;
140 4 : adfGeoTransform[5] = 1;
141 4 : }
142 :
143 : /************************************************************************/
144 : /* ~NGSGEOIDDataset() */
145 : /************************************************************************/
146 :
147 8 : NGSGEOIDDataset::~NGSGEOIDDataset()
148 :
149 : {
150 4 : FlushCache(true);
151 4 : if (fp)
152 4 : VSIFCloseL(fp);
153 8 : }
154 :
155 : /************************************************************************/
156 : /* GetHeaderInfo() */
157 : /************************************************************************/
158 :
159 2861 : int NGSGEOIDDataset::GetHeaderInfo(const GByte *pBuffer,
160 : double *padfGeoTransform, int *pnRows,
161 : int *pnCols, int *pbIsLittleEndian)
162 : {
163 : /* First check IKIND marker to determine if the file */
164 : /* is in little or big-endian order, and if it is a valid */
165 : /* NGSGEOID dataset */
166 : int nIKIND;
167 2861 : memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
168 2861 : CPL_LSBPTR32(&nIKIND);
169 2861 : if (nIKIND == 1)
170 : {
171 27 : *pbIsLittleEndian = TRUE;
172 : }
173 : else
174 : {
175 2834 : memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
176 2834 : CPL_MSBPTR32(&nIKIND);
177 2834 : if (nIKIND == 1)
178 : {
179 7 : *pbIsLittleEndian = FALSE;
180 : }
181 : else
182 : {
183 2827 : return FALSE;
184 : }
185 : }
186 :
187 : double dfSLAT;
188 34 : memcpy(&dfSLAT, pBuffer, 8);
189 34 : if (*pbIsLittleEndian)
190 : {
191 27 : CPL_LSBPTR64(&dfSLAT);
192 : }
193 : else
194 : {
195 7 : CPL_MSBPTR64(&dfSLAT);
196 : }
197 34 : pBuffer += 8;
198 :
199 : double dfWLON;
200 34 : memcpy(&dfWLON, pBuffer, 8);
201 34 : if (*pbIsLittleEndian)
202 : {
203 27 : CPL_LSBPTR64(&dfWLON);
204 : }
205 : else
206 : {
207 7 : CPL_MSBPTR64(&dfWLON);
208 : }
209 34 : pBuffer += 8;
210 :
211 : double dfDLAT;
212 34 : memcpy(&dfDLAT, pBuffer, 8);
213 34 : if (*pbIsLittleEndian)
214 : {
215 27 : CPL_LSBPTR64(&dfDLAT);
216 : }
217 : else
218 : {
219 7 : CPL_MSBPTR64(&dfDLAT);
220 : }
221 34 : pBuffer += 8;
222 :
223 : double dfDLON;
224 34 : memcpy(&dfDLON, pBuffer, 8);
225 34 : if (*pbIsLittleEndian)
226 : {
227 27 : CPL_LSBPTR64(&dfDLON);
228 : }
229 : else
230 : {
231 7 : CPL_MSBPTR64(&dfDLON);
232 : }
233 34 : pBuffer += 8;
234 :
235 : int nNLAT;
236 34 : memcpy(&nNLAT, pBuffer, 4);
237 34 : if (*pbIsLittleEndian)
238 : {
239 27 : CPL_LSBPTR32(&nNLAT);
240 : }
241 : else
242 : {
243 7 : CPL_MSBPTR32(&nNLAT);
244 : }
245 34 : pBuffer += 4;
246 :
247 : int nNLON;
248 34 : memcpy(&nNLON, pBuffer, 4);
249 34 : if (*pbIsLittleEndian)
250 : {
251 27 : CPL_LSBPTR32(&nNLON);
252 : }
253 : else
254 : {
255 7 : CPL_MSBPTR32(&nNLON);
256 : }
257 : /*pBuffer += 4;*/
258 :
259 : /*CPLDebug("NGSGEOID", "SLAT=%f, WLON=%f, DLAT=%f, DLON=%f, NLAT=%d,
260 : NLON=%d, IKIND=%d", dfSLAT, dfWLON, dfDLAT, dfDLON, nNLAT, nNLON,
261 : nIKIND);*/
262 :
263 34 : if (nNLAT <= 0 || nNLON <= 0 || dfDLAT <= 1e-15 || dfDLON <= 1e-15)
264 22 : return FALSE;
265 :
266 : /* Grids go over +180 in longitude */
267 : // Test written that way to be robust to NaN values
268 12 : if (!(dfSLAT >= -90.0 && dfSLAT + nNLAT * dfDLAT <= 90.0 &&
269 12 : dfWLON >= -180.0 && dfWLON + nNLON * dfDLON <= 360.0))
270 0 : return FALSE;
271 :
272 12 : padfGeoTransform[0] = dfWLON - dfDLON / 2;
273 12 : padfGeoTransform[1] = dfDLON;
274 12 : padfGeoTransform[2] = 0.0;
275 12 : padfGeoTransform[3] = dfSLAT + nNLAT * dfDLAT - dfDLAT / 2;
276 12 : padfGeoTransform[4] = 0.0;
277 12 : padfGeoTransform[5] = -dfDLAT;
278 :
279 12 : *pnRows = nNLAT;
280 12 : *pnCols = nNLON;
281 :
282 12 : return TRUE;
283 : }
284 :
285 : /************************************************************************/
286 : /* Identify() */
287 : /************************************************************************/
288 :
289 51458 : int NGSGEOIDDataset::Identify(GDALOpenInfo *poOpenInfo)
290 : {
291 51458 : if (poOpenInfo->nHeaderBytes < HEADER_SIZE)
292 48601 : return FALSE;
293 :
294 : double adfGeoTransform[6];
295 : int nRows, nCols;
296 : int bIsLittleEndian;
297 2857 : if (!GetHeaderInfo(poOpenInfo->pabyHeader, adfGeoTransform, &nRows, &nCols,
298 : &bIsLittleEndian))
299 2849 : return FALSE;
300 :
301 8 : return TRUE;
302 : }
303 :
304 : /************************************************************************/
305 : /* Open() */
306 : /************************************************************************/
307 :
308 4 : GDALDataset *NGSGEOIDDataset::Open(GDALOpenInfo *poOpenInfo)
309 :
310 : {
311 4 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
312 0 : return nullptr;
313 :
314 4 : if (poOpenInfo->eAccess == GA_Update)
315 : {
316 0 : ReportUpdateNotSupportedByDriver("NGSGEOID");
317 0 : return nullptr;
318 : }
319 :
320 : /* -------------------------------------------------------------------- */
321 : /* Create a corresponding GDALDataset. */
322 : /* -------------------------------------------------------------------- */
323 4 : NGSGEOIDDataset *poDS = new NGSGEOIDDataset();
324 4 : poDS->fp = poOpenInfo->fpL;
325 4 : poOpenInfo->fpL = nullptr;
326 :
327 4 : int nRows = 0, nCols = 0;
328 4 : GetHeaderInfo(poOpenInfo->pabyHeader, poDS->adfGeoTransform, &nRows, &nCols,
329 : &poDS->bIsLittleEndian);
330 4 : poDS->nRasterXSize = nCols;
331 4 : poDS->nRasterYSize = nRows;
332 :
333 : /* -------------------------------------------------------------------- */
334 : /* Create band information objects. */
335 : /* -------------------------------------------------------------------- */
336 4 : poDS->nBands = 1;
337 4 : poDS->SetBand(1, new NGSGEOIDRasterBand(poDS));
338 :
339 : /* -------------------------------------------------------------------- */
340 : /* Initialize any PAM information. */
341 : /* -------------------------------------------------------------------- */
342 4 : poDS->SetDescription(poOpenInfo->pszFilename);
343 4 : poDS->TryLoadXML();
344 :
345 : /* -------------------------------------------------------------------- */
346 : /* Support overviews. */
347 : /* -------------------------------------------------------------------- */
348 4 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
349 4 : return poDS;
350 : }
351 :
352 : /************************************************************************/
353 : /* GetGeoTransform() */
354 : /************************************************************************/
355 :
356 2 : CPLErr NGSGEOIDDataset::GetGeoTransform(double *padfTransform)
357 :
358 : {
359 2 : memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));
360 :
361 2 : return CE_None;
362 : }
363 :
364 : /************************************************************************/
365 : /* GetSpatialRef() */
366 : /************************************************************************/
367 :
368 2 : const OGRSpatialReference *NGSGEOIDDataset::GetSpatialRef() const
369 : {
370 2 : if (!m_oSRS.IsEmpty())
371 : {
372 0 : return &m_oSRS;
373 : }
374 :
375 : const CPLString osFilename =
376 6 : CPLString(CPLGetBasenameSafe(GetDescription())).tolower();
377 :
378 : // See https://www.ngs.noaa.gov/GEOID/GEOID12B/faq_2012B.shtml
379 :
380 : // GEOID2012 files ?
381 2 : if (STARTS_WITH(osFilename, "g2012") && osFilename.size() >= 7)
382 : {
383 0 : if (osFilename[6] == 'h' /* Hawai */ ||
384 0 : osFilename[6] == 's' /* Samoa */)
385 : {
386 : // NAD83 (PA11)
387 0 : m_oSRS.importFromEPSG(6322);
388 : }
389 0 : else if (osFilename[6] == 'g' /* Guam */)
390 : {
391 : // NAD83 (MA11)
392 0 : m_oSRS.importFromEPSG(6325);
393 : }
394 : else
395 : {
396 : // NAD83 (2011)
397 0 : m_oSRS.importFromEPSG(6318);
398 : }
399 :
400 0 : return &m_oSRS;
401 : }
402 :
403 : // USGG2012 files ? We should return IGS08, but there is only a
404 : // geocentric CRS in EPSG, so manually forge a geographic one from it
405 2 : if (STARTS_WITH(osFilename, "s2012"))
406 : {
407 0 : m_oSRS.importFromWkt(
408 : "GEOGCS[\"IGS08\",\n"
409 : " DATUM[\"IGS08\",\n"
410 : " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n"
411 : " AUTHORITY[\"EPSG\",\"7019\"]],\n"
412 : " AUTHORITY[\"EPSG\",\"1141\"]],\n"
413 : " PRIMEM[\"Greenwich\",0,\n"
414 : " AUTHORITY[\"EPSG\",\"8901\"]],\n"
415 : " UNIT[\"degree\",0.0174532925199433,\n"
416 : " AUTHORITY[\"EPSG\",\"9122\"]]]");
417 0 : return &m_oSRS;
418 : }
419 :
420 2 : m_oSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
421 2 : return &m_oSRS;
422 : }
423 :
424 : /************************************************************************/
425 : /* GDALRegister_NGSGEOID() */
426 : /************************************************************************/
427 :
428 1686 : void GDALRegister_NGSGEOID()
429 :
430 : {
431 1686 : if (GDALGetDriverByName("NGSGEOID") != nullptr)
432 302 : return;
433 :
434 1384 : GDALDriver *poDriver = new GDALDriver();
435 :
436 1384 : poDriver->SetDescription("NGSGEOID");
437 1384 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
438 1384 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NOAA NGS Geoid Height Grids");
439 1384 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
440 1384 : "drivers/raster/ngsgeoid.html");
441 1384 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "bin");
442 :
443 1384 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
444 :
445 1384 : poDriver->pfnOpen = NGSGEOIDDataset::Open;
446 1384 : poDriver->pfnIdentify = NGSGEOIDDataset::Identify;
447 :
448 1384 : GetGDALDriverManager()->RegisterDriver(poDriver);
449 : }
|