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 2791 : 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 2791 : memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
168 2791 : CPL_LSBPTR32(&nIKIND);
169 2791 : if (nIKIND == 1)
170 : {
171 27 : *pbIsLittleEndian = TRUE;
172 : }
173 : else
174 : {
175 2764 : memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
176 2764 : CPL_MSBPTR32(&nIKIND);
177 2764 : if (nIKIND == 1)
178 : {
179 7 : *pbIsLittleEndian = FALSE;
180 : }
181 : else
182 : {
183 2757 : 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 50643 : int NGSGEOIDDataset::Identify(GDALOpenInfo *poOpenInfo)
290 : {
291 50643 : if (poOpenInfo->nHeaderBytes < HEADER_SIZE)
292 47853 : return FALSE;
293 :
294 : double adfGeoTransform[6];
295 : int nRows, nCols;
296 : int bIsLittleEndian;
297 2790 : if (!GetHeaderInfo(poOpenInfo->pabyHeader, adfGeoTransform, &nRows, &nCols,
298 : &bIsLittleEndian))
299 2779 : 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 : CPLError(
317 : CE_Failure, CPLE_NotSupported,
318 : "The NGSGEOID driver does not support update access to existing"
319 : " datasets.\n");
320 0 : return nullptr;
321 : }
322 :
323 : /* -------------------------------------------------------------------- */
324 : /* Create a corresponding GDALDataset. */
325 : /* -------------------------------------------------------------------- */
326 4 : NGSGEOIDDataset *poDS = new NGSGEOIDDataset();
327 4 : poDS->fp = poOpenInfo->fpL;
328 4 : poOpenInfo->fpL = nullptr;
329 :
330 4 : int nRows = 0, nCols = 0;
331 4 : GetHeaderInfo(poOpenInfo->pabyHeader, poDS->adfGeoTransform, &nRows, &nCols,
332 : &poDS->bIsLittleEndian);
333 4 : poDS->nRasterXSize = nCols;
334 4 : poDS->nRasterYSize = nRows;
335 :
336 : /* -------------------------------------------------------------------- */
337 : /* Create band information objects. */
338 : /* -------------------------------------------------------------------- */
339 4 : poDS->nBands = 1;
340 4 : poDS->SetBand(1, new NGSGEOIDRasterBand(poDS));
341 :
342 : /* -------------------------------------------------------------------- */
343 : /* Initialize any PAM information. */
344 : /* -------------------------------------------------------------------- */
345 4 : poDS->SetDescription(poOpenInfo->pszFilename);
346 4 : poDS->TryLoadXML();
347 :
348 : /* -------------------------------------------------------------------- */
349 : /* Support overviews. */
350 : /* -------------------------------------------------------------------- */
351 4 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
352 4 : return poDS;
353 : }
354 :
355 : /************************************************************************/
356 : /* GetGeoTransform() */
357 : /************************************************************************/
358 :
359 2 : CPLErr NGSGEOIDDataset::GetGeoTransform(double *padfTransform)
360 :
361 : {
362 2 : memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));
363 :
364 2 : return CE_None;
365 : }
366 :
367 : /************************************************************************/
368 : /* GetSpatialRef() */
369 : /************************************************************************/
370 :
371 2 : const OGRSpatialReference *NGSGEOIDDataset::GetSpatialRef() const
372 : {
373 2 : if (!m_oSRS.IsEmpty())
374 : {
375 0 : return &m_oSRS;
376 : }
377 :
378 : const CPLString osFilename =
379 4 : CPLString(CPLGetBasename(GetDescription())).tolower();
380 :
381 : // See https://www.ngs.noaa.gov/GEOID/GEOID12B/faq_2012B.shtml
382 :
383 : // GEOID2012 files ?
384 2 : if (STARTS_WITH(osFilename, "g2012") && osFilename.size() >= 7)
385 : {
386 0 : if (osFilename[6] == 'h' /* Hawai */ ||
387 0 : osFilename[6] == 's' /* Samoa */)
388 : {
389 : // NAD83 (PA11)
390 0 : m_oSRS.importFromEPSG(6322);
391 : }
392 0 : else if (osFilename[6] == 'g' /* Guam */)
393 : {
394 : // NAD83 (MA11)
395 0 : m_oSRS.importFromEPSG(6325);
396 : }
397 : else
398 : {
399 : // NAD83 (2011)
400 0 : m_oSRS.importFromEPSG(6318);
401 : }
402 :
403 0 : return &m_oSRS;
404 : }
405 :
406 : // USGG2012 files ? We should return IGS08, but there is only a
407 : // geocentric CRS in EPSG, so manually forge a geographic one from it
408 2 : if (STARTS_WITH(osFilename, "s2012"))
409 : {
410 0 : m_oSRS.importFromWkt(
411 : "GEOGCS[\"IGS08\",\n"
412 : " DATUM[\"IGS08\",\n"
413 : " SPHEROID[\"GRS 1980\",6378137,298.257222101,\n"
414 : " AUTHORITY[\"EPSG\",\"7019\"]],\n"
415 : " AUTHORITY[\"EPSG\",\"1141\"]],\n"
416 : " PRIMEM[\"Greenwich\",0,\n"
417 : " AUTHORITY[\"EPSG\",\"8901\"]],\n"
418 : " UNIT[\"degree\",0.0174532925199433,\n"
419 : " AUTHORITY[\"EPSG\",\"9122\"]]]");
420 0 : return &m_oSRS;
421 : }
422 :
423 2 : m_oSRS.importFromWkt(SRS_WKT_WGS84_LAT_LONG);
424 2 : return &m_oSRS;
425 : }
426 :
427 : /************************************************************************/
428 : /* GDALRegister_NGSGEOID() */
429 : /************************************************************************/
430 :
431 1595 : void GDALRegister_NGSGEOID()
432 :
433 : {
434 1595 : if (GDALGetDriverByName("NGSGEOID") != nullptr)
435 302 : return;
436 :
437 1293 : GDALDriver *poDriver = new GDALDriver();
438 :
439 1293 : poDriver->SetDescription("NGSGEOID");
440 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
441 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NOAA NGS Geoid Height Grids");
442 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
443 1293 : "drivers/raster/ngsgeoid.html");
444 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "bin");
445 :
446 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
447 :
448 1293 : poDriver->pfnOpen = NGSGEOIDDataset::Open;
449 1293 : poDriver->pfnIdentify = NGSGEOIDDataset::Identify;
450 :
451 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
452 : }
|
