1 : /******************************************************************************
2 : * $Id: bagdataset.cpp 24756 2012-08-11 15:21:02Z 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 24756 2012-08-11 15:21:02Z 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 3 : H5Pclose(listid);
193 : }
194 :
195 : /* -------------------------------------------------------------------- */
196 : /* Load min/max information. */
197 : /* -------------------------------------------------------------------- */
198 3 : if( EQUAL(pszName,"elevation")
199 : && GH5_FetchAttribute( hDatasetID, "Maximum Elevation Value",
200 : dfMaximum )
201 : && GH5_FetchAttribute( hDatasetID, "Minimum Elevation Value",
202 : dfMinimum ) )
203 1 : bMinMaxSet = true;
204 2 : else if( EQUAL(pszName,"uncertainty")
205 : && GH5_FetchAttribute( hDatasetID, "Maximum Uncertainty Value",
206 : dfMaximum )
207 : && GH5_FetchAttribute( hDatasetID, "Minimum Uncertainty Value",
208 : dfMinimum ) )
209 1 : bMinMaxSet = true;
210 1 : else if( EQUAL(pszName,"nominal_elevation")
211 : && GH5_FetchAttribute( hDatasetID, "max_value",
212 : dfMaximum )
213 : && GH5_FetchAttribute( hDatasetID, "min_value",
214 : dfMinimum ) )
215 1 : bMinMaxSet = true;
216 :
217 3 : return true;
218 : }
219 :
220 : /************************************************************************/
221 : /* GetMinimum() */
222 : /************************************************************************/
223 :
224 1 : double BAGRasterBand::GetMinimum( int * pbSuccess )
225 :
226 : {
227 1 : if( bMinMaxSet )
228 : {
229 1 : if( pbSuccess )
230 1 : *pbSuccess = TRUE;
231 1 : return dfMinimum;
232 : }
233 : else
234 0 : return GDALRasterBand::GetMinimum( pbSuccess );
235 : }
236 :
237 : /************************************************************************/
238 : /* GetMaximum() */
239 : /************************************************************************/
240 :
241 1 : double BAGRasterBand::GetMaximum( int * pbSuccess )
242 :
243 : {
244 1 : if( bMinMaxSet )
245 : {
246 1 : if( pbSuccess )
247 1 : *pbSuccess = TRUE;
248 1 : return dfMaximum;
249 : }
250 : else
251 0 : return GDALRasterBand::GetMaximum( pbSuccess );
252 : }
253 :
254 : /************************************************************************/
255 : /* GetNoDataValue() */
256 : /************************************************************************/
257 3 : double BAGRasterBand::GetNoDataValue( int * pbSuccess )
258 :
259 : {
260 3 : if( pbSuccess )
261 3 : *pbSuccess = TRUE;
262 :
263 3 : if( EQUAL(GetDescription(),"elevation") )
264 1 : return 1000000.0;
265 2 : else if( EQUAL(GetDescription(),"uncertainty") )
266 1 : return 0.0;
267 1 : else if( EQUAL(GetDescription(),"nominal_elevation") )
268 1 : return 1000000.0;
269 : else
270 0 : return GDALPamRasterBand::GetNoDataValue( pbSuccess );
271 : }
272 :
273 : /************************************************************************/
274 : /* IReadBlock() */
275 : /************************************************************************/
276 30 : CPLErr BAGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
277 : void * pImage )
278 : {
279 : herr_t status;
280 : hsize_t count[3];
281 : H5OFFSET_TYPE offset[3];
282 : int nSizeOfData;
283 : hid_t memspace;
284 : hsize_t col_dims[3];
285 : hsize_t rank;
286 :
287 30 : rank=2;
288 :
289 30 : offset[0] = MAX(0,nRasterYSize - (nBlockYOff+1)*nBlockYSize);
290 30 : offset[1] = nBlockXOff*nBlockXSize;
291 30 : count[0] = nBlockYSize;
292 30 : count[1] = nBlockXSize;
293 :
294 30 : nSizeOfData = H5Tget_size( native );
295 30 : memset( pImage,0,nBlockXSize*nBlockYSize*nSizeOfData );
296 :
297 : /* blocksize may not be a multiple of imagesize */
298 30 : count[0] = MIN( size_t(nBlockYSize), GetYSize() - offset[0]);
299 30 : count[1] = MIN( size_t(nBlockXSize), GetXSize() - offset[1]);
300 :
301 30 : if( nRasterYSize - (nBlockYOff+1)*nBlockYSize < 0 )
302 : {
303 0 : count[0] += (nRasterYSize - (nBlockYOff+1)*nBlockYSize);
304 : }
305 :
306 : /* -------------------------------------------------------------------- */
307 : /* Select block from file space */
308 : /* -------------------------------------------------------------------- */
309 : status = H5Sselect_hyperslab( dataspace,
310 : H5S_SELECT_SET,
311 : offset, NULL,
312 30 : count, NULL );
313 :
314 : /* -------------------------------------------------------------------- */
315 : /* Create memory space to receive the data */
316 : /* -------------------------------------------------------------------- */
317 30 : col_dims[0]=nBlockYSize;
318 30 : col_dims[1]=nBlockXSize;
319 30 : memspace = H5Screate_simple( (int) rank, col_dims, NULL );
320 30 : H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
321 : status = H5Sselect_hyperslab(memspace,
322 : H5S_SELECT_SET,
323 : mem_offset, NULL,
324 30 : count, NULL);
325 :
326 : status = H5Dread ( hDatasetID,
327 : native,
328 : memspace,
329 : dataspace,
330 : H5P_DEFAULT,
331 30 : pImage );
332 :
333 30 : H5Sclose( memspace );
334 :
335 : /* -------------------------------------------------------------------- */
336 : /* Y flip the data. */
337 : /* -------------------------------------------------------------------- */
338 30 : int nLinesToFlip = count[0];
339 30 : int nLineSize = nSizeOfData * nBlockXSize;
340 30 : GByte *pabyTemp = (GByte *) CPLMalloc(nLineSize);
341 :
342 30 : for( int iY = 0; iY < nLinesToFlip/2; iY++ )
343 : {
344 : memcpy( pabyTemp,
345 : ((GByte *)pImage) + iY * nLineSize,
346 0 : nLineSize );
347 : memcpy( ((GByte *)pImage) + iY * nLineSize,
348 : ((GByte *)pImage) + (nLinesToFlip-iY-1) * nLineSize,
349 0 : nLineSize );
350 : memcpy( ((GByte *)pImage) + (nLinesToFlip-iY-1) * nLineSize,
351 : pabyTemp,
352 0 : nLineSize );
353 : }
354 :
355 30 : CPLFree( pabyTemp );
356 :
357 : /* -------------------------------------------------------------------- */
358 : /* Return success or failure. */
359 : /* -------------------------------------------------------------------- */
360 30 : if( status < 0 )
361 : {
362 : CPLError( CE_Failure, CPLE_AppDefined,
363 0 : "H5Dread() failed for block." );
364 0 : return CE_Failure;
365 : }
366 : else
367 30 : return CE_None;
368 : }
369 :
370 : /************************************************************************/
371 : /* ==================================================================== */
372 : /* BAGDataset */
373 : /* ==================================================================== */
374 : /************************************************************************/
375 :
376 : /************************************************************************/
377 : /* BAGDataset() */
378 : /************************************************************************/
379 :
380 1 : BAGDataset::BAGDataset()
381 : {
382 1 : hHDF5 = -1;
383 1 : pszXMLMetadata = NULL;
384 1 : pszProjection = NULL;
385 :
386 1 : adfGeoTransform[0] = 0.0;
387 1 : adfGeoTransform[1] = 1.0;
388 1 : adfGeoTransform[2] = 0.0;
389 1 : adfGeoTransform[3] = 0.0;
390 1 : adfGeoTransform[4] = 0.0;
391 1 : adfGeoTransform[5] = 1.0;
392 1 : }
393 :
394 : /************************************************************************/
395 : /* ~BAGDataset() */
396 : /************************************************************************/
397 1 : BAGDataset::~BAGDataset( )
398 : {
399 1 : FlushCache();
400 :
401 1 : if( hHDF5 >= 0 )
402 1 : H5Fclose( hHDF5 );
403 :
404 1 : CPLFree( pszXMLMetadata );
405 1 : CPLFree( pszProjection );
406 1 : }
407 :
408 : /************************************************************************/
409 : /* Identify() */
410 : /************************************************************************/
411 :
412 11824 : int BAGDataset::Identify( GDALOpenInfo * poOpenInfo )
413 :
414 : {
415 : /* -------------------------------------------------------------------- */
416 : /* Is it an HDF5 file? */
417 : /* -------------------------------------------------------------------- */
418 : static const char achSignature[] = "\211HDF\r\n\032\n";
419 :
420 11824 : if( poOpenInfo->pabyHeader == NULL
421 : || memcmp(poOpenInfo->pabyHeader,achSignature,8) != 0 )
422 11819 : return FALSE;
423 :
424 : /* -------------------------------------------------------------------- */
425 : /* Does it have the extension .bag? */
426 : /* -------------------------------------------------------------------- */
427 5 : if( !EQUAL(CPLGetExtension(poOpenInfo->pszFilename),"bag") )
428 4 : return FALSE;
429 :
430 1 : return TRUE;
431 : }
432 :
433 : /************************************************************************/
434 : /* Open() */
435 : /************************************************************************/
436 :
437 1733 : GDALDataset *BAGDataset::Open( GDALOpenInfo * poOpenInfo )
438 :
439 : {
440 : /* -------------------------------------------------------------------- */
441 : /* Confirm that this appears to be a BAG file. */
442 : /* -------------------------------------------------------------------- */
443 1733 : if( !Identify( poOpenInfo ) )
444 1732 : return NULL;
445 :
446 : /* -------------------------------------------------------------------- */
447 : /* Confirm the requested access is supported. */
448 : /* -------------------------------------------------------------------- */
449 1 : if( poOpenInfo->eAccess == GA_Update )
450 : {
451 : CPLError( CE_Failure, CPLE_NotSupported,
452 0 : "The BAG driver does not support update access." );
453 0 : return NULL;
454 : }
455 :
456 : /* -------------------------------------------------------------------- */
457 : /* Open the file as an HDF5 file. */
458 : /* -------------------------------------------------------------------- */
459 : hid_t hHDF5 = H5Fopen( poOpenInfo->pszFilename,
460 1 : H5F_ACC_RDONLY, H5P_DEFAULT );
461 :
462 1 : if( hHDF5 < 0 )
463 0 : return NULL;
464 :
465 : /* -------------------------------------------------------------------- */
466 : /* Confirm it is a BAG dataset by checking for the */
467 : /* BAG_Root/Bag Version attribute. */
468 : /* -------------------------------------------------------------------- */
469 1 : hid_t hBagRoot = H5Gopen( hHDF5, "/BAG_root" );
470 1 : hid_t hVersion = -1;
471 :
472 1 : if( hBagRoot >= 0 )
473 1 : hVersion = H5Aopen_name( hBagRoot, "Bag Version" );
474 :
475 1 : if( hVersion < 0 )
476 : {
477 0 : H5Fclose( hHDF5 );
478 0 : return NULL;
479 : }
480 1 : H5Aclose( hVersion );
481 :
482 : /* -------------------------------------------------------------------- */
483 : /* Create a corresponding dataset. */
484 : /* -------------------------------------------------------------------- */
485 1 : BAGDataset *poDS = new BAGDataset();
486 :
487 1 : poDS->hHDF5 = hHDF5;
488 :
489 : /* -------------------------------------------------------------------- */
490 : /* Extract version as metadata. */
491 : /* -------------------------------------------------------------------- */
492 1 : CPLString osVersion;
493 :
494 2 : if( GH5_FetchAttribute( hBagRoot, "Bag Version", osVersion ) )
495 1 : poDS->SetMetadataItem( "BagVersion", osVersion );
496 :
497 1 : H5Gclose( hBagRoot );
498 :
499 : /* -------------------------------------------------------------------- */
500 : /* Fetch the elevation dataset and attach as a band. */
501 : /* -------------------------------------------------------------------- */
502 1 : int nNextBand = 1;
503 1 : hid_t hElevation = H5Dopen( hHDF5, "/BAG_root/elevation" );
504 1 : if( hElevation < 0 )
505 : {
506 0 : delete poDS;
507 0 : return NULL;
508 : }
509 :
510 1 : BAGRasterBand *poElevBand = new BAGRasterBand( poDS, nNextBand );
511 :
512 2 : if( !poElevBand->Initialize( hElevation, "elevation" ) )
513 : {
514 0 : delete poElevBand;
515 0 : delete poDS;
516 0 : return NULL;
517 : }
518 :
519 1 : poDS->nRasterXSize = poElevBand->nRasterXSize;
520 1 : poDS->nRasterYSize = poElevBand->nRasterYSize;
521 :
522 1 : poDS->SetBand( nNextBand++, poElevBand );
523 :
524 : /* -------------------------------------------------------------------- */
525 : /* Try to do the same for the uncertainty band. */
526 : /* -------------------------------------------------------------------- */
527 1 : hid_t hUncertainty = H5Dopen( hHDF5, "/BAG_root/uncertainty" );
528 1 : BAGRasterBand *poUBand = new BAGRasterBand( poDS, nNextBand );
529 :
530 2 : if( hUncertainty >= 0 && poUBand->Initialize( hUncertainty, "uncertainty") )
531 : {
532 1 : poDS->SetBand( nNextBand++, poUBand );
533 : }
534 : else
535 0 : delete poUBand;
536 :
537 : /* -------------------------------------------------------------------- */
538 : /* Try to do the same for the uncertainty band. */
539 : /* -------------------------------------------------------------------- */
540 1 : hid_t hNominal = -1;
541 :
542 1 : H5E_BEGIN_TRY {
543 1 : hNominal = H5Dopen( hHDF5, "/BAG_root/nominal_elevation" );
544 1 : } H5E_END_TRY;
545 :
546 1 : BAGRasterBand *poNBand = new BAGRasterBand( poDS, nNextBand );
547 2 : if( hNominal >= 0 && poNBand->Initialize( hNominal,
548 : "nominal_elevation" ) )
549 : {
550 1 : poDS->SetBand( nNextBand++, poNBand );
551 : }
552 : else
553 0 : delete poNBand;
554 :
555 : /* -------------------------------------------------------------------- */
556 : /* Load the XML metadata. */
557 : /* -------------------------------------------------------------------- */
558 1 : poDS->LoadMetadata();
559 :
560 : /* -------------------------------------------------------------------- */
561 : /* Setup/check for pam .aux.xml. */
562 : /* -------------------------------------------------------------------- */
563 1 : poDS->SetDescription( poOpenInfo->pszFilename );
564 1 : poDS->TryLoadXML();
565 :
566 : /* -------------------------------------------------------------------- */
567 : /* Setup overviews. */
568 : /* -------------------------------------------------------------------- */
569 1 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
570 :
571 1 : return( poDS );
572 : }
573 :
574 : /************************************************************************/
575 : /* LoadMetadata() */
576 : /************************************************************************/
577 :
578 1 : void BAGDataset::LoadMetadata()
579 :
580 : {
581 : /* -------------------------------------------------------------------- */
582 : /* Load the metadata from the file. */
583 : /* -------------------------------------------------------------------- */
584 1 : hid_t hMDDS = H5Dopen( hHDF5, "/BAG_root/metadata" );
585 1 : hid_t datatype = H5Dget_type( hMDDS );
586 1 : hid_t dataspace = H5Dget_space( hMDDS );
587 1 : hid_t native = H5Tget_native_type( datatype, H5T_DIR_ASCEND );
588 : hsize_t dims[3], maxdims[3];
589 :
590 1 : H5Sget_simple_extent_dims( dataspace, dims, maxdims );
591 :
592 1 : pszXMLMetadata = (char *) CPLCalloc((int) (dims[0]+1),1);
593 :
594 1 : H5Dread( hMDDS, native, H5S_ALL, dataspace, H5P_DEFAULT, pszXMLMetadata );
595 :
596 1 : H5Tclose( native );
597 1 : H5Sclose( dataspace );
598 1 : H5Tclose( datatype );
599 1 : H5Dclose( hMDDS );
600 :
601 1 : if( strlen(pszXMLMetadata) == 0 )
602 0 : return;
603 :
604 : /* -------------------------------------------------------------------- */
605 : /* Try to get the geotransform. */
606 : /* -------------------------------------------------------------------- */
607 1 : CPLXMLNode *psRoot = CPLParseXMLString( pszXMLMetadata );
608 :
609 1 : if( psRoot == NULL )
610 0 : return;
611 :
612 1 : CPLStripXMLNamespace( psRoot, NULL, TRUE );
613 :
614 1 : CPLXMLNode *psGeo = CPLSearchXMLNode( psRoot, "=MD_Georectified" );
615 :
616 1 : if( psGeo != NULL )
617 : {
618 : char **papszCornerTokens =
619 : CSLTokenizeStringComplex(
620 : CPLGetXMLValue( psGeo, "cornerPoints.Point.coordinates", "" ),
621 1 : " ,", FALSE, FALSE );
622 :
623 1 : if( CSLCount(papszCornerTokens ) == 4 )
624 : {
625 1 : double dfLLX = atof( papszCornerTokens[0] );
626 1 : double dfLLY = atof( papszCornerTokens[1] );
627 1 : double dfURX = atof( papszCornerTokens[2] );
628 1 : double dfURY = atof( papszCornerTokens[3] );
629 :
630 1 : adfGeoTransform[0] = dfLLX;
631 1 : adfGeoTransform[1] = (dfURX - dfLLX) / (GetRasterXSize()-1);
632 1 : adfGeoTransform[3] = dfURY;
633 1 : adfGeoTransform[5] = (dfLLY - dfURY) / (GetRasterYSize()-1);
634 :
635 1 : adfGeoTransform[0] -= adfGeoTransform[1] * 0.5;
636 1 : adfGeoTransform[3] -= adfGeoTransform[5] * 0.5;
637 : }
638 1 : CSLDestroy( papszCornerTokens );
639 : }
640 :
641 : /* -------------------------------------------------------------------- */
642 : /* Try to get the coordinate system. */
643 : /* -------------------------------------------------------------------- */
644 1 : OGRSpatialReference oSRS;
645 :
646 1 : if( OGR_SRS_ImportFromISO19115( &oSRS, pszXMLMetadata )
647 : == OGRERR_NONE )
648 : {
649 0 : oSRS.exportToWkt( &pszProjection );
650 : }
651 : else
652 : {
653 1 : ParseWKTFromXML( pszXMLMetadata );
654 : }
655 :
656 : /* -------------------------------------------------------------------- */
657 : /* Fetch acquisition date. */
658 : /* -------------------------------------------------------------------- */
659 1 : CPLXMLNode *psDateTime = CPLSearchXMLNode( psRoot, "=dateTime" );
660 1 : if( psDateTime != NULL )
661 : {
662 1 : const char *pszDateTimeValue = CPLGetXMLValue( psDateTime, NULL, "" );
663 1 : if( pszDateTimeValue )
664 1 : SetMetadataItem( "BAG_DATETIME", pszDateTimeValue );
665 : }
666 :
667 1 : CPLDestroyXMLNode( psRoot );
668 : }
669 :
670 : /************************************************************************/
671 : /* ParseWKTFromXML() */
672 : /************************************************************************/
673 1 : OGRErr BAGDataset::ParseWKTFromXML( const char *pszISOXML )
674 : {
675 1 : OGRSpatialReference oSRS;
676 1 : CPLXMLNode *psRoot = CPLParseXMLString( pszISOXML );
677 1 : OGRErr eOGRErr = OGRERR_FAILURE;
678 :
679 1 : if( psRoot == NULL )
680 0 : return eOGRErr;
681 :
682 1 : CPLStripXMLNamespace( psRoot, NULL, TRUE );
683 :
684 1 : CPLXMLNode *psRSI = CPLSearchXMLNode( psRoot, "=referenceSystemInfo" );
685 1 : if( psRSI == NULL )
686 : {
687 : CPLError( CE_Failure, CPLE_AppDefined,
688 0 : "Unable to find <referenceSystemInfo> in metadata." );
689 0 : CPLDestroyXMLNode( psRoot );
690 0 : return eOGRErr;
691 : }
692 :
693 1 : oSRS.Clear();
694 :
695 : const char *pszSRCodeString =
696 1 : CPLGetXMLValue( psRSI, "MD_ReferenceSystem.referenceSystemIdentifier.RS_Identifier.code.CharacterString", NULL );
697 1 : if( pszSRCodeString == NULL )
698 : {
699 : CPLDebug("BAG",
700 1 : "Unable to find /MI_Metadata/referenceSystemInfo[1]/MD_ReferenceSystem[1]/referenceSystemIdentifier[1]/RS_Identifier[1]/code[1]/CharacterString[1] in metadata." );
701 1 : CPLDestroyXMLNode( psRoot );
702 1 : return eOGRErr;
703 : }
704 :
705 : const char *pszSRCodeSpace =
706 0 : CPLGetXMLValue( psRSI, "MD_ReferenceSystem.referenceSystemIdentifier.RS_Identifier.codeSpace.CharacterString", "" );
707 0 : if( !EQUAL( pszSRCodeSpace, "WKT" ) )
708 : {
709 : CPLError( CE_Failure, CPLE_AppDefined,
710 0 : "Spatial reference string is not in WKT." );
711 0 : CPLDestroyXMLNode( psRoot );
712 0 : return eOGRErr;
713 : }
714 :
715 0 : char* pszWKT = const_cast< char* >( pszSRCodeString );
716 0 : if( oSRS.importFromWkt( &pszWKT ) != OGRERR_NONE )
717 : {
718 : CPLError( CE_Failure, CPLE_AppDefined,
719 0 : "Failed parsing WKT string \"%s\".", pszSRCodeString );
720 0 : CPLDestroyXMLNode( psRoot );
721 0 : return eOGRErr;
722 : }
723 :
724 0 : oSRS.exportToWkt( &pszProjection );
725 0 : eOGRErr = OGRERR_NONE;
726 :
727 0 : psRSI = CPLSearchXMLNode( psRSI->psNext, "=referenceSystemInfo" );
728 0 : if( psRSI == NULL )
729 : {
730 : CPLError( CE_Failure, CPLE_AppDefined,
731 0 : "Unable to find second instance of <referenceSystemInfo> in metadata." );
732 0 : CPLDestroyXMLNode( psRoot );
733 0 : return eOGRErr;
734 : }
735 :
736 : pszSRCodeString =
737 0 : CPLGetXMLValue( psRSI, "MD_ReferenceSystem.referenceSystemIdentifier.RS_Identifier.code.CharacterString", NULL );
738 0 : if( pszSRCodeString == NULL )
739 : {
740 : CPLError( CE_Failure, CPLE_AppDefined,
741 0 : "Unable to find /MI_Metadata/referenceSystemInfo[2]/MD_ReferenceSystem[1]/referenceSystemIdentifier[1]/RS_Identifier[1]/code[1]/CharacterString[1] in metadata." );
742 0 : CPLDestroyXMLNode( psRoot );
743 0 : return eOGRErr;
744 : }
745 :
746 : pszSRCodeSpace =
747 0 : CPLGetXMLValue( psRSI, "MD_ReferenceSystem.referenceSystemIdentifier.RS_Identifier.codeSpace.CharacterString", "" );
748 0 : if( !EQUAL( pszSRCodeSpace, "WKT" ) )
749 : {
750 : CPLError( CE_Failure, CPLE_AppDefined,
751 0 : "Spatial reference string is not in WKT." );
752 0 : CPLDestroyXMLNode( psRoot );
753 0 : return eOGRErr;
754 : }
755 :
756 0 : if( EQUALN(pszSRCodeString, "VERTCS", 6 ) )
757 : {
758 0 : CPLString oString( pszProjection );
759 0 : oString += ",";
760 0 : oString += pszSRCodeString;
761 0 : if ( pszProjection )
762 0 : CPLFree( pszProjection );
763 0 : pszProjection = CPLStrdup( oString );
764 : }
765 :
766 0 : CPLDestroyXMLNode( psRoot );
767 :
768 0 : return eOGRErr;
769 : }
770 :
771 : /************************************************************************/
772 : /* GetGeoTransform() */
773 : /************************************************************************/
774 :
775 0 : CPLErr BAGDataset::GetGeoTransform( double *padfGeoTransform )
776 :
777 : {
778 0 : if( adfGeoTransform[0] != 0.0 || adfGeoTransform[3] != 0.0 )
779 : {
780 0 : memcpy( padfGeoTransform, adfGeoTransform, sizeof(double)*6 );
781 0 : return CE_None;
782 : }
783 : else
784 0 : return GDALPamDataset::GetGeoTransform( padfGeoTransform );
785 : }
786 :
787 : /************************************************************************/
788 : /* GetProjectionRef() */
789 : /************************************************************************/
790 :
791 0 : const char *BAGDataset::GetProjectionRef()
792 :
793 : {
794 0 : if( pszProjection )
795 0 : return pszProjection;
796 : else
797 0 : return GDALPamDataset::GetProjectionRef();
798 : }
799 :
800 : /************************************************************************/
801 : /* GetMetadata() */
802 : /************************************************************************/
803 :
804 1 : char **BAGDataset::GetMetadata( const char *pszDomain )
805 :
806 : {
807 1 : if( pszDomain != NULL && EQUAL(pszDomain,"xml:BAG") )
808 : {
809 1 : apszMDList[0] = pszXMLMetadata;
810 1 : apszMDList[1] = NULL;
811 :
812 1 : return apszMDList;
813 : }
814 : else
815 0 : return GDALPamDataset::GetMetadata( pszDomain );
816 : }
817 :
818 : /************************************************************************/
819 : /* GDALRegister_BAG() */
820 : /************************************************************************/
821 582 : void GDALRegister_BAG( )
822 :
823 : {
824 : GDALDriver *poDriver;
825 :
826 582 : if (! GDAL_CHECK_VERSION("BAG"))
827 0 : return;
828 :
829 582 : if( GDALGetDriverByName( "BAG" ) == NULL )
830 : {
831 561 : poDriver = new GDALDriver();
832 :
833 561 : poDriver->SetDescription( "BAG" );
834 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
835 561 : "Bathymetry Attributed Grid" );
836 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
837 561 : "frmt_bag.html" );
838 561 : poDriver->pfnOpen = BAGDataset::Open;
839 561 : poDriver->pfnIdentify = BAGDataset::Identify;
840 :
841 561 : GetGDALDriverManager( )->RegisterDriver( poDriver );
842 : }
843 : }
|
