1 : /******************************************************************************
2 : * $Id: ngsgeoiddataset.cpp 23334 2011-11-05 22:40:46Z rouault $
3 : *
4 : * Project: NGSGEOID driver
5 : * Purpose: GDALDataset driver for NGSGEOID dataset.
6 : * Author: Even Rouault, <even dot rouault at mines dash paris dot org>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2011, Even Rouault
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "cpl_vsi_virtual.h"
31 : #include "cpl_string.h"
32 : #include "gdal_pam.h"
33 : #include "ogr_srs_api.h"
34 :
35 : CPL_CVSID("$Id: ngsgeoiddataset.cpp 23334 2011-11-05 22:40:46Z rouault $");
36 :
37 : CPL_C_START
38 : void GDALRegister_NGSGEOID(void);
39 : CPL_C_END
40 :
41 : #define HEADER_SIZE (4 * 8 + 3 * 4)
42 :
43 : /************************************************************************/
44 : /* ==================================================================== */
45 : /* NGSGEOIDDataset */
46 : /* ==================================================================== */
47 : /************************************************************************/
48 :
49 : class NGSGEOIDRasterBand;
50 :
51 : class NGSGEOIDDataset : public GDALPamDataset
52 : {
53 : friend class NGSGEOIDRasterBand;
54 :
55 : VSILFILE *fp;
56 : double adfGeoTransform[6];
57 : int bIsLittleEndian;
58 :
59 : static int GetHeaderInfo( const GByte* pBuffer,
60 : double* padfGeoTransform,
61 : int* pnRows,
62 : int* pnCols,
63 : int* pbIsLittleEndian );
64 :
65 : public:
66 : NGSGEOIDDataset();
67 : virtual ~NGSGEOIDDataset();
68 :
69 : virtual CPLErr GetGeoTransform( double * );
70 : virtual const char* GetProjectionRef();
71 :
72 : static GDALDataset *Open( GDALOpenInfo * );
73 : static int Identify( GDALOpenInfo * );
74 : };
75 :
76 : /************************************************************************/
77 : /* ==================================================================== */
78 : /* NGSGEOIDRasterBand */
79 : /* ==================================================================== */
80 : /************************************************************************/
81 :
82 : class NGSGEOIDRasterBand : public GDALPamRasterBand
83 2 : {
84 : friend class NGSGEOIDDataset;
85 :
86 : public:
87 :
88 : NGSGEOIDRasterBand( NGSGEOIDDataset * );
89 :
90 : virtual CPLErr IReadBlock( int, int, void * );
91 :
92 0 : virtual const char* GetUnitType() { return "m"; }
93 : };
94 :
95 :
96 : /************************************************************************/
97 : /* NGSGEOIDRasterBand() */
98 : /************************************************************************/
99 :
100 2 : NGSGEOIDRasterBand::NGSGEOIDRasterBand( NGSGEOIDDataset *poDS )
101 :
102 : {
103 2 : this->poDS = poDS;
104 2 : this->nBand = 1;
105 :
106 2 : eDataType = GDT_Float32;
107 :
108 2 : nBlockXSize = poDS->GetRasterXSize();
109 2 : nBlockYSize = 1;
110 2 : }
111 :
112 : /************************************************************************/
113 : /* IReadBlock() */
114 : /************************************************************************/
115 :
116 2 : CPLErr NGSGEOIDRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
117 : void * pImage )
118 :
119 : {
120 2 : NGSGEOIDDataset *poGDS = (NGSGEOIDDataset *) poDS;
121 :
122 : /* First values in the file corresponds to the south-most line of the imagery */
123 : VSIFSeekL(poGDS->fp,
124 : HEADER_SIZE + (nRasterYSize - 1 - nBlockYOff) * nRasterXSize * 4,
125 2 : SEEK_SET);
126 :
127 2 : if ((int)VSIFReadL(pImage, 4, nRasterXSize, poGDS->fp) != nRasterXSize)
128 0 : return CE_Failure;
129 :
130 2 : if (poGDS->bIsLittleEndian)
131 : {
132 : #ifdef CPL_MSB
133 : GDALSwapWords( pImage, 4, nRasterXSize, 4 );
134 : #endif
135 : }
136 : else
137 : {
138 : #ifdef CPL_LSB
139 1 : GDALSwapWords( pImage, 4, nRasterXSize, 4 );
140 : #endif
141 : }
142 :
143 2 : return CE_None;
144 : }
145 :
146 : /************************************************************************/
147 : /* ~NGSGEOIDDataset() */
148 : /************************************************************************/
149 :
150 2 : NGSGEOIDDataset::NGSGEOIDDataset()
151 : {
152 2 : fp = NULL;
153 2 : adfGeoTransform[0] = 0;
154 2 : adfGeoTransform[1] = 1;
155 2 : adfGeoTransform[2] = 0;
156 2 : adfGeoTransform[3] = 0;
157 2 : adfGeoTransform[4] = 0;
158 2 : adfGeoTransform[5] = 1;
159 2 : bIsLittleEndian = TRUE;
160 2 : }
161 :
162 : /************************************************************************/
163 : /* ~NGSGEOIDDataset() */
164 : /************************************************************************/
165 :
166 2 : NGSGEOIDDataset::~NGSGEOIDDataset()
167 :
168 : {
169 2 : FlushCache();
170 2 : if (fp)
171 2 : VSIFCloseL(fp);
172 2 : }
173 :
174 : /************************************************************************/
175 : /* GetHeaderInfo() */
176 : /************************************************************************/
177 :
178 173 : int NGSGEOIDDataset::GetHeaderInfo( const GByte* pBuffer,
179 : double* padfGeoTransform,
180 : int* pnRows,
181 : int* pnCols,
182 : int* pbIsLittleEndian )
183 : {
184 : double dfSLAT;
185 : double dfWLON;
186 : double dfDLAT;
187 : double dfDLON;
188 : int nNLAT;
189 : int nNLON;
190 : int nIKIND;
191 :
192 : /* First check IKIND marker to determine if the file */
193 : /* is in little or big-endian order, and if it is a valid */
194 : /* NGSGEOID dataset */
195 173 : memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
196 : CPL_LSBPTR32(&nIKIND);
197 173 : if (nIKIND == 1)
198 : {
199 2 : *pbIsLittleEndian = TRUE;
200 : }
201 : else
202 : {
203 171 : memcpy(&nIKIND, pBuffer + HEADER_SIZE - 4, 4);
204 171 : CPL_MSBPTR32(&nIKIND);
205 171 : if (nIKIND == 1)
206 : {
207 2 : *pbIsLittleEndian = FALSE;
208 : }
209 : else
210 : {
211 169 : return FALSE;
212 : }
213 : }
214 :
215 4 : memcpy(&dfSLAT, pBuffer, 8);
216 4 : if (*pbIsLittleEndian)
217 : {
218 : CPL_LSBPTR64(&dfSLAT);
219 : }
220 : else
221 : {
222 2 : CPL_MSBPTR64(&dfSLAT);
223 : }
224 4 : pBuffer += 8;
225 4 : memcpy(&dfWLON, pBuffer, 8);
226 4 : if (*pbIsLittleEndian)
227 : {
228 : CPL_LSBPTR64(&dfWLON);
229 : }
230 : else
231 : {
232 2 : CPL_MSBPTR64(&dfWLON);
233 : }
234 4 : pBuffer += 8;
235 4 : memcpy(&dfDLAT, pBuffer, 8);
236 4 : if (*pbIsLittleEndian)
237 : {
238 : CPL_LSBPTR64(&dfDLAT);
239 : }
240 : else
241 : {
242 2 : CPL_MSBPTR64(&dfDLAT);
243 : }
244 4 : pBuffer += 8;
245 4 : memcpy(&dfDLON, pBuffer, 8);
246 4 : if (*pbIsLittleEndian)
247 : {
248 : CPL_LSBPTR64(&dfDLON);
249 : }
250 : else
251 : {
252 2 : CPL_MSBPTR64(&dfDLON);
253 : }
254 4 : pBuffer += 8;
255 4 : memcpy(&nNLAT, pBuffer, 4);
256 4 : if (*pbIsLittleEndian)
257 : {
258 : CPL_LSBPTR32(&nNLAT);
259 : }
260 : else
261 : {
262 2 : CPL_MSBPTR32(&nNLAT);
263 : }
264 4 : pBuffer += 4;
265 4 : memcpy(&nNLON, pBuffer, 4);
266 4 : if (*pbIsLittleEndian)
267 : {
268 : CPL_LSBPTR32(&nNLON);
269 : }
270 : else
271 : {
272 2 : CPL_MSBPTR32(&nNLON);
273 : }
274 4 : pBuffer += 4;
275 :
276 : /*CPLDebug("NGSGEOID", "SLAT=%f, WLON=%f, DLAT=%f, DLON=%f, NLAT=%d, NLON=%d, IKIND=%d",
277 : dfSLAT, dfWLON, dfDLAT, dfDLON, nNLAT, nNLON, nIKIND);*/
278 :
279 4 : if (nNLAT <= 0 || nNLON <= 0 || dfDLAT <= 0.0 || dfDLON <= 0.0)
280 0 : return FALSE;
281 :
282 : /* Grids go over +180 in longitude */
283 4 : if (dfSLAT < -90.0 || dfSLAT + nNLAT * dfDLAT > 90.0 ||
284 : dfWLON < -180.0 || dfWLON + nNLON * dfDLON > 360.0)
285 0 : return FALSE;
286 :
287 4 : padfGeoTransform[0] = dfWLON - dfDLON / 2;
288 4 : padfGeoTransform[1] = dfDLON;
289 4 : padfGeoTransform[2] = 0.0;
290 4 : padfGeoTransform[3] = dfSLAT + nNLAT * dfDLAT - dfDLAT / 2;
291 4 : padfGeoTransform[4] = 0.0;
292 4 : padfGeoTransform[5] = -dfDLAT;
293 :
294 4 : *pnRows = nNLAT;
295 4 : *pnCols = nNLON;
296 :
297 4 : return TRUE;
298 : }
299 :
300 : /************************************************************************/
301 : /* Identify() */
302 : /************************************************************************/
303 :
304 11384 : int NGSGEOIDDataset::Identify( GDALOpenInfo * poOpenInfo )
305 : {
306 11384 : if (poOpenInfo->nHeaderBytes < HEADER_SIZE)
307 11213 : return FALSE;
308 :
309 : double adfGeoTransform[6];
310 : int nRows, nCols;
311 : int bIsLittleEndian;
312 171 : if ( !GetHeaderInfo( poOpenInfo->pabyHeader,
313 : adfGeoTransform,
314 : &nRows, &nCols, &bIsLittleEndian ) )
315 169 : return FALSE;
316 :
317 2 : return TRUE;
318 : }
319 :
320 :
321 : /************************************************************************/
322 : /* Open() */
323 : /************************************************************************/
324 :
325 1320 : GDALDataset *NGSGEOIDDataset::Open( GDALOpenInfo * poOpenInfo )
326 :
327 : {
328 1320 : if (!Identify(poOpenInfo))
329 1318 : return NULL;
330 :
331 2 : if (poOpenInfo->eAccess == GA_Update)
332 : {
333 : CPLError( CE_Failure, CPLE_NotSupported,
334 : "The NGSGEOID driver does not support update access to existing"
335 0 : " datasets.\n" );
336 0 : return NULL;
337 : }
338 :
339 2 : VSILFILE* fp = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
340 2 : if (fp == NULL)
341 0 : return NULL;
342 :
343 : /* -------------------------------------------------------------------- */
344 : /* Create a corresponding GDALDataset. */
345 : /* -------------------------------------------------------------------- */
346 : NGSGEOIDDataset *poDS;
347 :
348 2 : poDS = new NGSGEOIDDataset();
349 2 : poDS->fp = fp;
350 :
351 : int nRows, nCols;
352 : GetHeaderInfo( poOpenInfo->pabyHeader,
353 : poDS->adfGeoTransform,
354 : &nRows,
355 : &nCols,
356 2 : &poDS->bIsLittleEndian );
357 2 : poDS->nRasterXSize = nCols;
358 2 : poDS->nRasterYSize = nRows;
359 :
360 : /* -------------------------------------------------------------------- */
361 : /* Create band information objects. */
362 : /* -------------------------------------------------------------------- */
363 2 : poDS->nBands = 1;
364 4 : poDS->SetBand( 1, new NGSGEOIDRasterBand( poDS ) );
365 :
366 : /* -------------------------------------------------------------------- */
367 : /* Initialize any PAM information. */
368 : /* -------------------------------------------------------------------- */
369 2 : poDS->SetDescription( poOpenInfo->pszFilename );
370 2 : poDS->TryLoadXML();
371 :
372 : /* -------------------------------------------------------------------- */
373 : /* Support overviews. */
374 : /* -------------------------------------------------------------------- */
375 2 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
376 2 : return( poDS );
377 : }
378 :
379 : /************************************************************************/
380 : /* GetGeoTransform() */
381 : /************************************************************************/
382 :
383 2 : CPLErr NGSGEOIDDataset::GetGeoTransform( double * padfTransform )
384 :
385 : {
386 2 : memcpy(padfTransform, adfGeoTransform, 6 * sizeof(double));
387 :
388 2 : return( CE_None );
389 : }
390 :
391 : /************************************************************************/
392 : /* GetProjectionRef() */
393 : /************************************************************************/
394 :
395 2 : const char* NGSGEOIDDataset::GetProjectionRef()
396 : {
397 2 : return SRS_WKT_WGS84;
398 : }
399 :
400 : /************************************************************************/
401 : /* GDALRegister_NGSGEOID() */
402 : /************************************************************************/
403 :
404 582 : void GDALRegister_NGSGEOID()
405 :
406 : {
407 : GDALDriver *poDriver;
408 :
409 582 : if( GDALGetDriverByName( "NGSGEOID" ) == NULL )
410 : {
411 561 : poDriver = new GDALDriver();
412 :
413 561 : poDriver->SetDescription( "NGSGEOID" );
414 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
415 561 : "NOAA NGS Geoid Height Grids" );
416 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
417 561 : "frmt_ngsgeoid.html" );
418 561 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "bin" );
419 :
420 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
421 :
422 561 : poDriver->pfnOpen = NGSGEOIDDataset::Open;
423 561 : poDriver->pfnIdentify = NGSGEOIDDataset::Identify;
424 :
425 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
426 : }
427 582 : }
428 :
|
