1 : /******************************************************************************
2 : * $Id: bagdataset.cpp 25774 2013-03-20 20:19:13Z rouault $
3 : *
4 : * Project: Hierarchical Data Format Release 5 (HDF5)
5 : * Purpose: Read BAG datasets.
6 : * Author: Frank Warmerdam <warmerdam@pobox.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2009, Frank Warmerdam <warmerdam@pobox.com>
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 "gh5_convenience.h"
31 :
32 : #include "gdal_pam.h"
33 : #include "gdal_priv.h"
34 : #include "ogr_spatialref.h"
35 : #include "cpl_string.h"
36 :
37 : CPL_CVSID("$Id: bagdataset.cpp 25774 2013-03-20 20:19:13Z rouault $");
38 :
39 : CPL_C_START
40 : void GDALRegister_BAG(void);
41 : CPL_C_END
42 :
43 : OGRErr OGR_SRS_ImportFromISO19115( OGRSpatialReference *poThis,
44 : const char *pszISOXML );
45 :
46 : /************************************************************************/
47 : /* ==================================================================== */
48 : /* BAGDataset */
49 : /* ==================================================================== */
50 : /************************************************************************/
51 : class BAGDataset : public GDALPamDataset
52 : {
53 :
54 : friend class BAGRasterBand;
55 :
56 : hid_t hHDF5;
57 :
58 : char *pszProjection;
59 : double adfGeoTransform[6];
60 :
61 : void LoadMetadata();
62 :
63 : char *pszXMLMetadata;
64 : char *apszMDList[2];
65 :
66 : public:
67 : BAGDataset();
68 : ~BAGDataset();
69 :
70 : virtual CPLErr GetGeoTransform( double * );
71 : virtual const char *GetProjectionRef(void);
72 : virtual char **GetMetadata( const char * pszDomain = "" );
73 :
74 : static GDALDataset *Open( GDALOpenInfo * );
75 : static int Identify( GDALOpenInfo * );
76 :
77 : OGRErr ParseWKTFromXML( const char *pszISOXML );
78 : };
79 :
80 : /************************************************************************/
81 : /* ==================================================================== */
82 : /* BAGRasterBand */
83 : /* ==================================================================== */
84 : /************************************************************************/
85 : class BAGRasterBand : public GDALPamRasterBand
86 : {
87 : friend class BAGDataset;
88 :
89 : hid_t hDatasetID;
90 : hid_t native;
91 : hid_t dataspace;
92 :
93 : bool bMinMaxSet;
94 : double dfMinimum;
95 : double dfMaximum;
96 :
97 : public:
98 :
99 : BAGRasterBand( BAGDataset *, int );
100 : ~BAGRasterBand();
101 :
102 : bool Initialize( hid_t hDataset, const char *pszName );
103 :
104 : virtual CPLErr IReadBlock( int, int, void * );
105 : virtual double GetNoDataValue( int * );
106 :
107 : virtual double GetMinimum( int *pbSuccess = NULL );
108 : virtual double GetMaximum(int *pbSuccess = NULL );
109 : };
110 :
111 : /************************************************************************/
112 : /* BAGRasterBand() */
113 : /************************************************************************/
114 3 : BAGRasterBand::BAGRasterBand( BAGDataset *poDS, int nBand )
115 :
116 : {
117 3 : this->poDS = poDS;
118 3 : this->nBand = nBand;
119 :
120 3 : hDatasetID = -1;
121 3 : dataspace = -1;
122 3 : native = -1;
123 3 : bMinMaxSet = false;
124 3 : }
125 :
126 : /************************************************************************/
127 : /* ~BAGRasterBand() */
128 : /************************************************************************/
129 :
130 3 : BAGRasterBand::~BAGRasterBand()
131 : {
132 3 : if( dataspace > 0 )
133 3 : H5Sclose(dataspace);
134 :
135 3 : if( native > 0 )
136 3 : H5Tclose( native );
137 :
138 3 : if( hDatasetID > 0 )
139 3 : H5Dclose( hDatasetID );
140 3 : }
141 :
142 : /************************************************************************/
143 : /* Initialize() */
144 : /************************************************************************/
145 :
146 3 : bool BAGRasterBand::Initialize( hid_t hDatasetID, const char *pszName )
147 :
148 : {
149 3 : SetDescription( pszName );
150 :
151 3 : this->hDatasetID = hDatasetID;
152 :
153 3 : hid_t datatype = H5Dget_type( hDatasetID );
154 3 : dataspace = H5Dget_space( hDatasetID );
155 3 : hid_t n_dims = H5Sget_simple_extent_ndims( dataspace );
156 3 : native = H5Tget_native_type( datatype, H5T_DIR_ASCEND );
157 : hsize_t dims[3], maxdims[3];
158 :
159 3 : eDataType = GH5_GetDataType( native );
160 :
161 3 : if( n_dims == 2 )
162 : {
163 3 : H5Sget_simple_extent_dims( dataspace, dims, maxdims );
164 :
165 3 : nRasterXSize = (int) dims[1];
166 3 : nRasterYSize = (int) dims[0];
167 : }
168 : else
169 : {
170 : CPLError( CE_Failure, CPLE_AppDefined,
171 0 : "Dataset not of rank 2." );
172 0 : return false;
173 : }
174 :
175 3 : nBlockXSize = nRasterXSize;
176 3 : nBlockYSize = 1;
177 :
178 : /* -------------------------------------------------------------------- */
179 : /* Check for chunksize, and use it as blocksize for optimized */
180 : /* reading. */
181 : /* -------------------------------------------------------------------- */
182 3 : hid_t listid = H5Dget_create_plist( hDatasetID );
183 3 : if (listid>0)
184 : {
185 3 : if(H5Pget_layout(listid) == H5D_CHUNKED)
186 : {
187 : hsize_t panChunkDims[3];
188 0 : int nDimSize = H5Pget_chunk(listid, 3, panChunkDims);
189 0 : nBlockXSize = (int) panChunkDims[nDimSize-1];
190 0 : nBlockYSize = (int) panChunkDims[nDimSize-2];
191 : }
192 :
193 3 : int nfilters = H5Pget_nfilters( listid );
194 :
195 : H5Z_filter_t filter;
196 : char name[120];
197 3 : size_t cd_nelmts = 20;
198 : unsigned int cd_values[20];
199 : unsigned int flags;
200 3 : for (int i = 0; i < nfilters; i++)
201 : {
202 0 : filter = H5Pget_filter(listid, i, &flags, (size_t *)&cd_nelmts, cd_values, 120, name);
203 0 : if (filter == H5Z_FILTER_DEFLATE)
204 0 : poDS->SetMetadataItem( "COMPRESSION", "DEFLATE", "IMAGE_STRUCTURE" );
205 0 : else if (filter == H5Z_FILTER_NBIT)
206 0 : poDS->SetMetadataItem( "COMPRESSION", "NBIT", "IMAGE_STRUCTURE" );
207 0 : else if (filter == H5Z_FILTER_SCALEOFFSET)
208 0 : poDS->SetMetadataItem( "COMPRESSION", "SCALEOFFSET", "IMAGE_STRUCTURE" );
209 0 : else if (filter == H5Z_FILTER_SZIP)
210 0 : poDS->SetMetadataItem( "COMPRESSION", "SZIP", "IMAGE_STRUCTURE" );
211 : }
212 :
213 3 : H5Pclose(listid);
214 : }
215 :
216 : /* -------------------------------------------------------------------- */
217 : /* Load min/max information. */
218 : /* -------------------------------------------------------------------- */
219 3 : if( EQUAL(pszName,"elevation")
220 : && GH5_FetchAttribute( hDatasetID, "Maximum Elevation Value",
221 : dfMaximum )
222 : && GH5_FetchAttribute( hDatasetID, "Minimum Elevation Value",
223 : dfMinimum ) )
224 1 : bMinMaxSet = true;
225 2 : else if( EQUAL(pszName,"uncertainty")
226 : && GH5_FetchAttribute( hDatasetID, "Maximum Uncertainty Value",
227 : dfMaximum )
228 : && GH5_FetchAttribute( hDatasetID, "Minimum Uncertainty Value",
229 : dfMinimum ) )
230 1 : bMinMaxSet = true;
231 1 : else if( EQUAL(pszName,"nominal_elevation")
232 : && GH5_FetchAttribute( hDatasetID, "max_value",
233 : dfMaximum )
234 : && GH5_FetchAttribute( hDatasetID, "min_value",
235 : dfMinimum ) )
236 1 : bMinMaxSet = true;
237 :
238 3 : return true;
239 : }
240 :
241 : /************************************************************************/
242 : /* GetMinimum() */
243 : /************************************************************************/
244 :
245 1 : double BAGRasterBand::GetMinimum( int * pbSuccess )
246 :
247 : {
248 1 : if( bMinMaxSet )
249 : {
250 1 : if( pbSuccess )
251 1 : *pbSuccess = TRUE;
252 1 : return dfMinimum;
253 : }
254 : else
255 0 : return GDALRasterBand::GetMinimum( pbSuccess );
256 : }
257 :
258 : /************************************************************************/
259 : /* GetMaximum() */
260 : /************************************************************************/
261 :
262 1 : double BAGRasterBand::GetMaximum( int * pbSuccess )
263 :
264 : {
265 1 : if( bMinMaxSet )
266 : {
267 1 : if( pbSuccess )
268 1 : *pbSuccess = TRUE;
269 1 : return dfMaximum;
270 : }
271 : else
272 0 : return GDALRasterBand::GetMaximum( pbSuccess );
273 : }
274 :
275 : /************************************************************************/
276 : /* GetNoDataValue() */
277 : /************************************************************************/
278 3 : double BAGRasterBand::GetNoDataValue( int * pbSuccess )
279 :
280 : {
281 3 : if( pbSuccess )
282 3 : *pbSuccess = TRUE;
283 :
284 3 : if( EQUAL(GetDescription(),"elevation") )
285 1 : return 1000000.0;
286 2 : else if( EQUAL(GetDescription(),"uncertainty") )
287 1 : return 0.0;
288 1 : else if( EQUAL(GetDescription(),"nominal_elevation") )
289 1 : return 1000000.0;
290 : else
291 0 : return GDALPamRasterBand::GetNoDataValue( pbSuccess );
292 : }
293 :
294 : /************************************************************************/
295 : /* IReadBlock() */
296 : /************************************************************************/
297 30 : CPLErr BAGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
298 : void * pImage )
299 : {
300 : herr_t status;
301 : hsize_t count[3];
302 : H5OFFSET_TYPE offset[3];
303 : int nSizeOfData;
304 : hid_t memspace;
305 : hsize_t col_dims[3];
306 : hsize_t rank;
307 :
308 30 : rank=2;
309 :
310 30 : offset[0] = MAX(0,nRasterYSize - (nBlockYOff+1)*nBlockYSize);
311 30 : offset[1] = nBlockXOff*nBlockXSize;
312 30 : count[0] = nBlockYSize;
313 30 : count[1] = nBlockXSize;
314 :
315 30 : nSizeOfData = H5Tget_size( native );
316 30 : memset( pImage,0,nBlockXSize*nBlockYSize*nSizeOfData );
317 :
318 : /* blocksize may not be a multiple of imagesize */
319 30 : count[0] = MIN( size_t(nBlockYSize), GetYSize() - offset[0]);
320 30 : count[1] = MIN( size_t(nBlockXSize), GetXSize() - offset[1]);
321 :
322 30 : if( nRasterYSize - (nBlockYOff+1)*nBlockYSize < 0 )
323 : {
324 0 : count[0] += (nRasterYSize - (nBlockYOff+1)*nBlockYSize);
325 : }
326 :
327 : /* -------------------------------------------------------------------- */
328 : /* Select block from file space */
329 : /* -------------------------------------------------------------------- */
330 : status = H5Sselect_hyperslab( dataspace,
331 : H5S_SELECT_SET,
332 : offset, NULL,
333 30 : count, NULL );
334 :
335 : /* -------------------------------------------------------------------- */
336 : /* Create memory space to receive the data */
337 : /* -------------------------------------------------------------------- */
338 30 : col_dims[0]=nBlockYSize;
339 30 : col_dims[1]=nBlockXSize;
340 30 : memspace = H5Screate_simple( (int) rank, col_dims, NULL );
341 30 : H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
342 : status = H5Sselect_hyperslab(memspace,
343 : H5S_SELECT_SET,
344 : mem_offset, NULL,
345 30 : count, NULL);
346 :
347 : status = H5Dread ( hDatasetID,
348 : native,
349 : memspace,
350 : dataspace,
351 : H5P_DEFAULT,
352 30 : pImage );
353 :
354 30 : H5Sclose( memspace );
355 :
356 : /* -------------------------------------------------------------------- */
357 : /* Y flip the data. */
358 : /* -------------------------------------------------------------------- */
359 30 : int nLinesToFlip = count[0];
360 30 : int nLineSize = nSizeOfData * nBlockXSize;
361 30 : GByte *pabyTemp = (GByte *) CPLMalloc(nLineSize);
362 :
363 30 : for( int iY = 0; iY < nLinesToFlip/2; iY++ )
364 : {
365 : memcpy( pabyTemp,
366 : ((GByte *)pImage) + iY * nLineSize,
367 0 : nLineSize );
368 : memcpy( ((GByte *)pImage) + iY * nLineSize,
369 : ((GByte *)pImage) + (nLinesToFlip-iY-1) * nLineSize,
370 0 : nLineSize );
371 : memcpy( ((GByte *)pImage) + (nLinesToFlip-iY-1) * nLineSize,
372 : pabyTemp,
373 0 : nLineSize );
374 : }
375 :
376 30 : CPLFree( pabyTemp );
377 :
378 : /* -------------------------------------------------------------------- */
379 : /* Return success or failure. */
380 : /* -------------------------------------------------------------------- */
381 30 : if( status < 0 )
382 : {
383 : CPLError( CE_Failure, CPLE_AppDefined,
384 0 : "H5Dread() failed for block." );
385 0 : return CE_Failure;
386 : }
387 : else
388 30 : return CE_None;
389 : }
390 :
391 : /************************************************************************/
392 : /* ==================================================================== */
393 : /* BAGDataset */
394 : /* ==================================================================== */
395 : /************************************************************************/
396 :
397 : /************************************************************************/
398 : /* BAGDataset() */
399 : /************************************************************************/
400 :
401 1 : BAGDataset::BAGDataset()
402 : {
403 1 : hHDF5 = -1;
404 1 : pszXMLMetadata = NULL;
405 1 : pszProjection = NULL;
406 :
407 1 : adfGeoTransform[0] = 0.0;
408 1 : adfGeoTransform[1] = 1.0;
409 1 : adfGeoTransform[2] = 0.0;
410 1 : adfGeoTransform[3] = 0.0;
411 1 : adfGeoTransform[4] = 0.0;
412 1 : adfGeoTransform[5] = 1.0;
413 1 : }
414 :
415 : /************************************************************************/
416 : /* ~BAGDataset() */
417 : /************************************************************************/
418 1 : BAGDataset::~BAGDataset( )
419 : {
420 1 : FlushCache();
421 :
422 1 : if( hHDF5 >= 0 )
423 1 : H5Fclose( hHDF5 );
424 :
425 1 : CPLFree( pszXMLMetadata );
426 1 : CPLFree( pszProjection );
427 1 : }
428 :
429 : /************************************************************************/
430 : /* Identify() */
431 : /************************************************************************/
432 :
433 11867 : int BAGDataset::Identify( GDALOpenInfo * poOpenInfo )
434 :
435 : {
436 : /* -------------------------------------------------------------------- */
437 : /* Is it an HDF5 file? */
438 : /* -------------------------------------------------------------------- */
439 : static const char achSignature[] = "\211HDF\r\n\032\n";
440 :
441 11867 : if( poOpenInfo->pabyHeader == NULL
442 : || memcmp(poOpenInfo->pabyHeader,achSignature,8) != 0 )
443 11862 : return FALSE;
444 :
445 : /* -------------------------------------------------------------------- */
446 : /* Does it have the extension .bag? */
447 : /* -------------------------------------------------------------------- */
448 5 : if( !EQUAL(CPLGetExtension(poOpenInfo->pszFilename),"bag") )
449 4 : return FALSE;
450 :
451 1 : return TRUE;
452 : }
453 :
454 : /************************************************************************/
455 : /* Open() */
456 : /************************************************************************/
457 :
458 1738 : GDALDataset *BAGDataset::Open( GDALOpenInfo * poOpenInfo )
459 :
460 : {
461 : /* -------------------------------------------------------------------- */
462 : /* Confirm that this appears to be a BAG file. */
463 : /* -------------------------------------------------------------------- */
464 1738 : if( !Identify( poOpenInfo ) )
465 1737 : return NULL;
466 :
467 : /* -------------------------------------------------------------------- */
468 : /* Confirm the requested access is supported. */
469 : /* -------------------------------------------------------------------- */
470 1 : if( poOpenInfo->eAccess == GA_Update )
471 : {
472 : CPLError( CE_Failure, CPLE_NotSupported,
473 0 : "The BAG driver does not support update access." );
474 0 : return NULL;
475 : }
476 :
477 : /* -------------------------------------------------------------------- */
478 : /* Open the file as an HDF5 file. */
479 : /* -------------------------------------------------------------------- */
480 : hid_t hHDF5 = H5Fopen( poOpenInfo->pszFilename,
481 1 : H5F_ACC_RDONLY, H5P_DEFAULT );
482 :
483 1 : if( hHDF5 < 0 )
484 0 : return NULL;
485 :
486 : /* -------------------------------------------------------------------- */
487 : /* Confirm it is a BAG dataset by checking for the */
488 : /* BAG_Root/Bag Version attribute. */
489 : /* -------------------------------------------------------------------- */
490 1 : hid_t hBagRoot = H5Gopen( hHDF5, "/BAG_root" );
491 1 : hid_t hVersion = -1;
492 :
493 1 : if( hBagRoot >= 0 )
494 1 : hVersion = H5Aopen_name( hBagRoot, "Bag Version" );
495 :
496 1 : if( hVersion < 0 )
497 : {
498 0 : H5Fclose( hHDF5 );
499 0 : return NULL;
500 : }
501 1 : H5Aclose( hVersion );
502 :
503 : /* -------------------------------------------------------------------- */
504 : /* Create a corresponding dataset. */
505 : /* -------------------------------------------------------------------- */
506 1 : BAGDataset *poDS = new BAGDataset();
507 :
508 1 : poDS->hHDF5 = hHDF5;
509 :
510 : /* -------------------------------------------------------------------- */
511 : /* Extract version as metadata. */
512 : /* -------------------------------------------------------------------- */
513 1 : CPLString osVersion;
514 :
515 2 : if( GH5_FetchAttribute( hBagRoot, "Bag Version", osVersion ) )
516 1 : poDS->SetMetadataItem( "BagVersion", osVersion );
517 :
518 1 : H5Gclose( hBagRoot );
519 :
520 : /* -------------------------------------------------------------------- */
521 : /* Fetch the elevation dataset and attach as a band. */
522 : /* -------------------------------------------------------------------- */
523 1 : int nNextBand = 1;
524 1 : hid_t hElevation = H5Dopen( hHDF5, "/BAG_root/elevation" );
525 1 : if( hElevation < 0 )
526 : {
527 0 : delete poDS;
528 0 : return NULL;
529 : }
530 :
531 1 : BAGRasterBand *poElevBand = new BAGRasterBand( poDS, nNextBand );
532 :
533 2 : if( !poElevBand->Initialize( hElevation, "elevation" ) )
534 : {
535 0 : delete poElevBand;
536 0 : delete poDS;
537 0 : return NULL;
538 : }
539 :
540 1 : poDS->nRasterXSize = poElevBand->nRasterXSize;
541 1 : poDS->nRasterYSize = poElevBand->nRasterYSize;
542 :
543 1 : poDS->SetBand( nNextBand++, poElevBand );
544 :
545 : /* -------------------------------------------------------------------- */
546 : /* Try to do the same for the uncertainty band. */
547 : /* -------------------------------------------------------------------- */
548 1 : hid_t hUncertainty = H5Dopen( hHDF5, "/BAG_root/uncertainty" );
549 1 : BAGRasterBand *poUBand = new BAGRasterBand( poDS, nNextBand );
550 :
551 2 : if( hUncertainty >= 0 && poUBand->Initialize( hUncertainty, "uncertainty") )
552 : {
553 1 : poDS->SetBand( nNextBand++, poUBand );
554 : }
555 : else
556 0 : delete poUBand;
557 :
558 : /* -------------------------------------------------------------------- */
559 : /* Try to do the same for the uncertainty band. */
560 : /* -------------------------------------------------------------------- */
561 1 : hid_t hNominal = -1;
562 :
563 1 : H5E_BEGIN_TRY {
564 1 : hNominal = H5Dopen( hHDF5, "/BAG_root/nominal_elevation" );
565 1 : } H5E_END_TRY;
566 :
567 1 : BAGRasterBand *poNBand = new BAGRasterBand( poDS, nNextBand );
568 2 : if( hNominal >= 0 && poNBand->Initialize( hNominal,
569 : "nominal_elevation" ) )
570 : {
571 1 : poDS->SetBand( nNextBand++, poNBand );
572 : }
573 : else
574 0 : delete poNBand;
575 :
576 : /* -------------------------------------------------------------------- */
577 : /* Load the XML metadata. */
578 : /* -------------------------------------------------------------------- */
579 1 : poDS->LoadMetadata();
580 :
581 : /* -------------------------------------------------------------------- */
582 : /* Setup/check for pam .aux.xml. */
583 : /* -------------------------------------------------------------------- */
584 1 : poDS->SetDescription( poOpenInfo->pszFilename );
585 1 : poDS->TryLoadXML();
586 :
587 : /* -------------------------------------------------------------------- */
588 : /* Setup overviews. */
589 : /* -------------------------------------------------------------------- */
590 1 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
591 :
592 1 : return( poDS );
593 : }
594 :
595 : /************************************************************************/
596 : /* LoadMetadata() */
597 : /************************************************************************/
598 :
599 1 : void BAGDataset::LoadMetadata()
600 :
601 : {
602 : /* -------------------------------------------------------------------- */
603 : /* Load the metadata from the file. */
604 : /* -------------------------------------------------------------------- */
605 1 : hid_t hMDDS = H5Dopen( hHDF5, "/BAG_root/metadata" );
606 1 : hid_t datatype = H5Dget_type( hMDDS );
607 1 : hid_t dataspace = H5Dget_space( hMDDS );
608 1 : hid_t native = H5Tget_native_type( datatype, H5T_DIR_ASCEND );
609 : hsize_t dims[3], maxdims[3];
610 :
611 1 : H5Sget_simple_extent_dims( dataspace, dims, maxdims );
612 :
613 1 : pszXMLMetadata = (char *) CPLCalloc((int) (dims[0]+1),1);
614 :
615 1 : H5Dread( hMDDS, native, H5S_ALL, dataspace, H5P_DEFAULT, pszXMLMetadata );
616 :
617 1 : H5Tclose( native );
618 1 : H5Sclose( dataspace );
619 1 : H5Tclose( datatype );
620 1 : H5Dclose( hMDDS );
621 :
622 1 : if( strlen(pszXMLMetadata) == 0 )
623 0 : return;
624 :
625 : /* -------------------------------------------------------------------- */
626 : /* Try to get the geotransform. */
627 : /* -------------------------------------------------------------------- */
628 1 : CPLXMLNode *psRoot = CPLParseXMLString( pszXMLMetadata );
629 :
630 1 : if( psRoot == NULL )
631 0 : return;
632 :
633 1 : CPLStripXMLNamespace( psRoot, NULL, TRUE );
634 :
635 1 : CPLXMLNode *psGeo = CPLSearchXMLNode( psRoot, "=MD_Georectified" );
636 :
637 1 : if( psGeo != NULL )
638 : {
639 : char **papszCornerTokens =
640 : CSLTokenizeStringComplex(
641 : CPLGetXMLValue( psGeo, "cornerPoints.Point.coordinates", "" ),
642 1 : " ,", FALSE, FALSE );
643 :
644 1 : if( CSLCount(papszCornerTokens ) == 4 )
645 : {
646 1 : double dfLLX = atof( papszCornerTokens[0] );
647 1 : double dfLLY = atof( papszCornerTokens[1] );
648 1 : double dfURX = atof( papszCornerTokens[2] );
649 1 : double dfURY = atof( papszCornerTokens[3] );
650 :
651 1 : adfGeoTransform[0] = dfLLX;
652 1 : adfGeoTransform[1] = (dfURX - dfLLX) / (GetRasterXSize()-1);
653 1 : adfGeoTransform[3] = dfURY;
654 1 : adfGeoTransform[5] = (dfLLY - dfURY) / (GetRasterYSize()-1);
655 :
656 1 : adfGeoTransform[0] -= adfGeoTransform[1] * 0.5;
657 1 : adfGeoTransform[3] -= adfGeoTransform[5] * 0.5;
658 : }
659 1 : CSLDestroy( papszCornerTokens );
660 : }
661 :
662 : /* -------------------------------------------------------------------- */
663 : /* Try to get the coordinate system. */
664 : /* -------------------------------------------------------------------- */
665 1 : OGRSpatialReference oSRS;
666 :
667 1 : if( OGR_SRS_ImportFromISO19115( &oSRS, pszXMLMetadata )
668 : == OGRERR_NONE )
669 : {
670 0 : oSRS.exportToWkt( &pszProjection );
671 : }
672 : else
673 : {
674 1 : ParseWKTFromXML( pszXMLMetadata );
675 : }
676 :
677 : /* -------------------------------------------------------------------- */
678 : /* Fetch acquisition date. */
679 : /* -------------------------------------------------------------------- */
680 1 : CPLXMLNode *psDateTime = CPLSearchXMLNode( psRoot, "=dateTime" );
681 1 : if( psDateTime != NULL )
682 : {
683 1 : const char *pszDateTimeValue = CPLGetXMLValue( psDateTime, NULL, "" );
684 1 : if( pszDateTimeValue )
685 1 : SetMetadataItem( "BAG_DATETIME", pszDateTimeValue );
686 : }
687 :
688 1 : CPLDestroyXMLNode( psRoot );
689 : }
690 :
691 : /************************************************************************/
692 : /* ParseWKTFromXML() */
693 : /************************************************************************/
694 1 : OGRErr BAGDataset::ParseWKTFromXML( const char *pszISOXML )
695 : {
696 1 : OGRSpatialReference oSRS;
697 1 : CPLXMLNode *psRoot = CPLParseXMLString( pszISOXML );
698 1 : OGRErr eOGRErr = OGRERR_FAILURE;
699 :
700 1 : if( psRoot == NULL )
701 0 : return eOGRErr;
702 :
703 1 : CPLStripXMLNamespace( psRoot, NULL, TRUE );
704 :
705 1 : CPLXMLNode *psRSI = CPLSearchXMLNode( psRoot, "=referenceSystemInfo" );
706 1 : if( psRSI == NULL )
707 : {
708 : CPLError( CE_Failure, CPLE_AppDefined,
709 0 : "Unable to find <referenceSystemInfo> in metadata." );
710 0 : CPLDestroyXMLNode( psRoot );
711 0 : return eOGRErr;
712 : }
713 :
714 1 : oSRS.Clear();
715 :
716 : const char *pszSRCodeString =
717 1 : CPLGetXMLValue( psRSI, "MD_ReferenceSystem.referenceSystemIdentifier.RS_Identifier.code.CharacterString", NULL );
718 1 : if( pszSRCodeString == NULL )
719 : {
720 : CPLDebug("BAG",
721 1 : "Unable to find /MI_Metadata/referenceSystemInfo[1]/MD_ReferenceSystem[1]/referenceSystemIdentifier[1]/RS_Identifier[1]/code[1]/CharacterString[1] in metadata." );
722 1 : CPLDestroyXMLNode( psRoot );
723 1 : return eOGRErr;
724 : }
725 :
726 : const char *pszSRCodeSpace =
727 0 : CPLGetXMLValue( psRSI, "MD_ReferenceSystem.referenceSystemIdentifier.RS_Identifier.codeSpace.CharacterString", "" );
728 0 : if( !EQUAL( pszSRCodeSpace, "WKT" ) )
729 : {
730 : CPLError( CE_Failure, CPLE_AppDefined,
731 0 : "Spatial reference string is not in WKT." );
732 0 : CPLDestroyXMLNode( psRoot );
733 0 : return eOGRErr;
734 : }
735 :
736 0 : char* pszWKT = const_cast< char* >( pszSRCodeString );
737 0 : if( oSRS.importFromWkt( &pszWKT ) != OGRERR_NONE )
738 : {
739 : CPLError( CE_Failure, CPLE_AppDefined,
740 0 : "Failed parsing WKT string \"%s\".", pszSRCodeString );
741 0 : CPLDestroyXMLNode( psRoot );
742 0 : return eOGRErr;
743 : }
744 :
745 0 : oSRS.exportToWkt( &pszProjection );
746 0 : eOGRErr = OGRERR_NONE;
747 :
748 0 : psRSI = CPLSearchXMLNode( psRSI->psNext, "=referenceSystemInfo" );
749 0 : if( psRSI == NULL )
750 : {
751 : CPLError( CE_Failure, CPLE_AppDefined,
752 0 : "Unable to find second instance of <referenceSystemInfo> in metadata." );
753 0 : CPLDestroyXMLNode( psRoot );
754 0 : return eOGRErr;
755 : }
756 :
757 : pszSRCodeString =
758 0 : CPLGetXMLValue( psRSI, "MD_ReferenceSystem.referenceSystemIdentifier.RS_Identifier.code.CharacterString", NULL );
759 0 : if( pszSRCodeString == NULL )
760 : {
761 : CPLError( CE_Failure, CPLE_AppDefined,
762 0 : "Unable to find /MI_Metadata/referenceSystemInfo[2]/MD_ReferenceSystem[1]/referenceSystemIdentifier[1]/RS_Identifier[1]/code[1]/CharacterString[1] in metadata." );
763 0 : CPLDestroyXMLNode( psRoot );
764 0 : return eOGRErr;
765 : }
766 :
767 : pszSRCodeSpace =
768 0 : CPLGetXMLValue( psRSI, "MD_ReferenceSystem.referenceSystemIdentifier.RS_Identifier.codeSpace.CharacterString", "" );
769 0 : if( !EQUAL( pszSRCodeSpace, "WKT" ) )
770 : {
771 : CPLError( CE_Failure, CPLE_AppDefined,
772 0 : "Spatial reference string is not in WKT." );
773 0 : CPLDestroyXMLNode( psRoot );
774 0 : return eOGRErr;
775 : }
776 :
777 0 : if( EQUALN(pszSRCodeString, "VERTCS", 6 ) )
778 : {
779 0 : CPLString oString( pszProjection );
780 0 : oString += ",";
781 0 : oString += pszSRCodeString;
782 0 : if ( pszProjection )
783 0 : CPLFree( pszProjection );
784 0 : pszProjection = CPLStrdup( oString );
785 : }
786 :
787 0 : CPLDestroyXMLNode( psRoot );
788 :
789 0 : return eOGRErr;
790 : }
791 :
792 : /************************************************************************/
793 : /* GetGeoTransform() */
794 : /************************************************************************/
795 :
796 0 : CPLErr BAGDataset::GetGeoTransform( double *padfGeoTransform )
797 :
798 : {
799 0 : if( adfGeoTransform[0] != 0.0 || adfGeoTransform[3] != 0.0 )
800 : {
801 0 : memcpy( padfGeoTransform, adfGeoTransform, sizeof(double)*6 );
802 0 : return CE_None;
803 : }
804 : else
805 0 : return GDALPamDataset::GetGeoTransform( padfGeoTransform );
806 : }
807 :
808 : /************************************************************************/
809 : /* GetProjectionRef() */
810 : /************************************************************************/
811 :
812 0 : const char *BAGDataset::GetProjectionRef()
813 :
814 : {
815 0 : if( pszProjection )
816 0 : return pszProjection;
817 : else
818 0 : return GDALPamDataset::GetProjectionRef();
819 : }
820 :
821 : /************************************************************************/
822 : /* GetMetadata() */
823 : /************************************************************************/
824 :
825 1 : char **BAGDataset::GetMetadata( const char *pszDomain )
826 :
827 : {
828 1 : if( pszDomain != NULL && EQUAL(pszDomain,"xml:BAG") )
829 : {
830 1 : apszMDList[0] = pszXMLMetadata;
831 1 : apszMDList[1] = NULL;
832 :
833 1 : return apszMDList;
834 : }
835 : else
836 0 : return GDALPamDataset::GetMetadata( pszDomain );
837 : }
838 :
839 : /************************************************************************/
840 : /* GDALRegister_BAG() */
841 : /************************************************************************/
842 610 : void GDALRegister_BAG( )
843 :
844 : {
845 : GDALDriver *poDriver;
846 :
847 610 : if (! GDAL_CHECK_VERSION("BAG"))
848 0 : return;
849 :
850 610 : if( GDALGetDriverByName( "BAG" ) == NULL )
851 : {
852 588 : poDriver = new GDALDriver();
853 :
854 588 : poDriver->SetDescription( "BAG" );
855 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
856 588 : "Bathymetry Attributed Grid" );
857 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
858 588 : "frmt_bag.html" );
859 588 : poDriver->pfnOpen = BAGDataset::Open;
860 588 : poDriver->pfnIdentify = BAGDataset::Identify;
861 :
862 588 : GetGDALDriverManager( )->RegisterDriver( poDriver );
863 : }
864 : }
|
