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