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