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 "gdal_driver.h"
18 : #include "gdal_drivermanager.h"
19 : #include "gdal_openinfo.h"
20 : #include "gdal_cpp_functions.h"
21 : #include "ogr_srs_api.h"
22 :
23 : #define HEADER_SIZE (4 * 8 + 3 * 4)
24 :
25 : /************************************************************************/
26 : /* ==================================================================== */
27 : /* NGSGEOIDDataset */
28 : /* ==================================================================== */
29 : /************************************************************************/
30 :
31 : class NGSGEOIDRasterBand;
32 :
33 : class NGSGEOIDDataset final : public GDALPamDataset
34 : {
35 : friend class NGSGEOIDRasterBand;
36 :
37 : VSILFILE *fp{};
38 : GDALGeoTransform m_gt{};
39 : int bIsLittleEndian{};
40 : mutable OGRSpatialReference m_oSRS{};
41 :
42 : static int GetHeaderInfo(const GByte *pBuffer, GDALGeoTransform >,
43 : int *pnRows, int *pnCols, int *pbIsLittleEndian);
44 :
45 : CPL_DISALLOW_COPY_ASSIGN(NGSGEOIDDataset)
46 :
47 : public:
48 : NGSGEOIDDataset();
49 : ~NGSGEOIDDataset() override;
50 :
51 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
52 : const OGRSpatialReference *GetSpatialRef() const override;
53 :
54 : static GDALDataset *Open(GDALOpenInfo *);
55 : static int Identify(GDALOpenInfo *);
56 : };
57 :
58 : /************************************************************************/
59 : /* ==================================================================== */
60 : /* NGSGEOIDRasterBand */
61 : /* ==================================================================== */
62 : /************************************************************************/
63 :
64 : class NGSGEOIDRasterBand final : public GDALPamRasterBand
65 : {
66 : friend class NGSGEOIDDataset;
67 :
68 : public:
69 : explicit NGSGEOIDRasterBand(NGSGEOIDDataset *);
70 :
71 : CPLErr IReadBlock(int, int, void *) override;
72 :
73 0 : const char *GetUnitType() override
74 : {
75 0 : return "m";
76 : }
77 : };
78 :
79 : /************************************************************************/
80 : /* NGSGEOIDRasterBand() */
81 : /************************************************************************/
82 :
83 4 : NGSGEOIDRasterBand::NGSGEOIDRasterBand(NGSGEOIDDataset *poDSIn)
84 :
85 : {
86 4 : poDS = poDSIn;
87 4 : nBand = 1;
88 :
89 4 : eDataType = GDT_Float32;
90 :
91 4 : nBlockXSize = poDS->GetRasterXSize();
92 4 : nBlockYSize = 1;
93 4 : }
94 :
95 : /************************************************************************/
96 : /* IReadBlock() */
97 : /************************************************************************/
98 :
99 2 : CPLErr NGSGEOIDRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
100 : void *pImage)
101 :
102 : {
103 2 : NGSGEOIDDataset *poGDS = cpl::down_cast<NGSGEOIDDataset *>(poDS);
104 :
105 : /* First values in the file corresponds to the south-most line of the
106 : * imagery */
107 2 : VSIFSeekL(poGDS->fp,
108 2 : HEADER_SIZE +
109 2 : static_cast<vsi_l_offset>(nRasterYSize - 1 - nBlockYOff) *
110 2 : nRasterXSize * 4,
111 : SEEK_SET);
112 :
113 2 : if (static_cast<int>(VSIFReadL(pImage, 4, nRasterXSize, poGDS->fp)) !=
114 2 : nRasterXSize)
115 0 : return CE_Failure;
116 :
117 : #ifdef CPL_MSB
118 : if (poGDS->bIsLittleEndian)
119 : {
120 : GDALSwapWords(pImage, 4, nRasterXSize, 4);
121 : }
122 : #endif
123 :
124 : #ifdef CPL_LSB
125 2 : if (!poGDS->bIsLittleEndian)
126 : {
127 1 : GDALSwapWords(pImage, 4, nRasterXSize, 4);
128 : }
129 : #endif
130 :
131 2 : return CE_None;
132 : }
133 :
134 : /************************************************************************/
135 : /* ~NGSGEOIDDataset() */
136 : /************************************************************************/
137 :
138 4 : NGSGEOIDDataset::NGSGEOIDDataset() : fp(nullptr), bIsLittleEndian(TRUE)
139 : {
140 4 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
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 3100 : int NGSGEOIDDataset::GetHeaderInfo(const GByte *pBuffer, GDALGeoTransform >,
160 : int *pnRows, int *pnCols,
161 : 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 3100 : memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
168 3100 : CPL_LSBPTR32(&nIKIND);
169 3100 : if (nIKIND == 1)
170 : {
171 27 : *pbIsLittleEndian = TRUE;
172 : }
173 : else
174 : {
175 3073 : memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
176 3073 : CPL_MSBPTR32(&nIKIND);
177 3073 : if (nIKIND == 1)
178 : {
179 7 : *pbIsLittleEndian = FALSE;
180 : }
181 : else
182 : {
183 3066 : 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 : gt[0] = dfWLON - dfDLON / 2;
273 12 : gt[1] = dfDLON;
274 12 : gt[2] = 0.0;
275 12 : gt[3] = dfSLAT + nNLAT * dfDLAT - dfDLAT / 2;
276 12 : gt[4] = 0.0;
277 12 : gt[5] = -dfDLAT;
278 :
279 12 : *pnRows = nNLAT;
280 12 : *pnCols = nNLON;
281 :
282 12 : return TRUE;
283 : }
284 :
285 : /************************************************************************/
286 : /* Identify() */
287 : /************************************************************************/
288 :
289 58290 : int NGSGEOIDDataset::Identify(GDALOpenInfo *poOpenInfo)
290 : {
291 58290 : if (poOpenInfo->nHeaderBytes < HEADER_SIZE)
292 55196 : return FALSE;
293 :
294 3094 : GDALGeoTransform gt;
295 : int nRows, nCols;
296 : int bIsLittleEndian;
297 3094 : if (!GetHeaderInfo(poOpenInfo->pabyHeader, gt, &nRows, &nCols,
298 : &bIsLittleEndian))
299 3088 : 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->m_gt, &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(GDALGeoTransform >) const
357 :
358 : {
359 2 : gt = m_gt;
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 2038 : void GDALRegister_NGSGEOID()
429 :
430 : {
431 2038 : if (GDALGetDriverByName("NGSGEOID") != nullptr)
432 283 : return;
433 :
434 1755 : GDALDriver *poDriver = new GDALDriver();
435 :
436 1755 : poDriver->SetDescription("NGSGEOID");
437 1755 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
438 1755 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NOAA NGS Geoid Height Grids");
439 1755 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
440 1755 : "drivers/raster/ngsgeoid.html");
441 1755 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "bin");
442 :
443 1755 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
444 :
445 1755 : poDriver->pfnOpen = NGSGEOIDDataset::Open;
446 1755 : poDriver->pfnIdentify = NGSGEOIDDataset::Identify;
447 :
448 1755 : GetGDALDriverManager()->RegisterDriver(poDriver);
449 : }
|
