1 : /******************************************************************************
2 : * $Id: hdf5imagedataset.cpp 25801 2013-03-25 17:26:49Z antonio $
3 : *
4 : * Project: Hierarchical Data Format Release 5 (HDF5)
5 : * Purpose: Read subdatasets of HDF5 file.
6 : * Author: Denis Nadeau <denis.nadeau@gmail.com>
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2005, 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 : #define H5_USE_16_API
31 :
32 : #include "hdf5.h"
33 :
34 : #include "gdal_pam.h"
35 : #include "gdal_priv.h"
36 : #include "cpl_string.h"
37 : #include "hdf5dataset.h"
38 : #include "ogr_spatialref.h"
39 :
40 : CPL_CVSID("$Id: hdf5imagedataset.cpp 25801 2013-03-25 17:26:49Z antonio $");
41 :
42 : CPL_C_START
43 : void GDALRegister_HDF5Image(void);
44 : CPL_C_END
45 :
46 : /* release 1.6.3 or 1.6.4 changed the type of count in some api functions */
47 :
48 : #if H5_VERS_MAJOR == 1 && H5_VERS_MINOR <= 6 \
49 : && (H5_VERS_MINOR < 6 || H5_VERS_RELEASE < 3)
50 : # define H5OFFSET_TYPE hssize_t
51 : #else
52 : # define H5OFFSET_TYPE hsize_t
53 : #endif
54 :
55 : class HDF5ImageDataset : public HDF5Dataset
56 : {
57 : typedef enum
58 : {
59 : UNKNOWN_PRODUCT=0,
60 : CSK_PRODUCT
61 : } Hdf5ProductType;
62 :
63 : typedef enum
64 : {
65 : PROD_UNKNOWN=0,
66 : PROD_CSK_L0,
67 : PROD_CSK_L1A,
68 : PROD_CSK_L1B,
69 : PROD_CSK_L1C,
70 : PROD_CSK_L1D
71 : } HDF5CSKProductEnum;
72 :
73 : friend class HDF5ImageRasterBand;
74 :
75 : char *pszProjection;
76 : char *pszGCPProjection;
77 : GDAL_GCP *pasGCPList;
78 : int nGCPCount;
79 : OGRSpatialReference oSRS;
80 :
81 : hsize_t *dims,*maxdims;
82 : HDF5GroupObjects *poH5Objects;
83 : int ndims,dimensions;
84 : hid_t dataset_id;
85 : hid_t dataspace_id;
86 : hsize_t size;
87 : haddr_t address;
88 : hid_t datatype;
89 : hid_t native;
90 : H5T_class_t clas;
91 : int iSubdatasetType;
92 : double adfGeoTransform[6];
93 : bool bHasGeoTransform;
94 :
95 : public:
96 : HDF5ImageDataset();
97 : ~HDF5ImageDataset();
98 :
99 : CPLErr CreateProjections( );
100 : static GDALDataset *Open( GDALOpenInfo * );
101 : static int Identify( GDALOpenInfo * );
102 :
103 : const char *GetProjectionRef();
104 : virtual int GetGCPCount( );
105 : virtual const char *GetGCPProjection();
106 : virtual const GDAL_GCP *GetGCPs( );
107 : virtual CPLErr GetGeoTransform( double * padfTransform );
108 :
109 : /**
110 : * Identify if the subdataset has a known product format
111 : * It stores a product identifier in iSubdatasetType,
112 : * UNKNOWN_PRODUCT, if it isn't a recognizable format.
113 : */
114 : void IdentifyProductType();
115 :
116 : /**
117 : * Captures Geolocation information from a COSMO-SKYMED
118 : * file.
119 : * The geoid will allways be WGS84
120 : * The projection type may be UTM or UPS, depending on the
121 : * latitude from the center of the image.
122 : * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
123 : */
124 : void CaptureCSKGeolocation(int iProductType);
125 :
126 : /**
127 : * Get Geotransform information for COSMO-SKYMED files
128 : * In case of sucess it stores the transformation
129 : * in adfGeoTransform. In case of failure it doesn't
130 : * modify adfGeoTransform
131 : * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
132 : */
133 : void CaptureCSKGeoTransform(int iProductType);
134 :
135 : /**
136 : * @param iProductType type of HDF5 subproduct, see HDF5CSKProduct
137 : */
138 : void CaptureCSKGCPs(int iProductType);
139 : };
140 :
141 : /************************************************************************/
142 : /* ==================================================================== */
143 : /* HDF5ImageDataset */
144 : /* ==================================================================== */
145 : /************************************************************************/
146 :
147 : /************************************************************************/
148 : /* HDF5ImageDataset() */
149 : /************************************************************************/
150 14 : HDF5ImageDataset::HDF5ImageDataset()
151 : {
152 14 : nGCPCount = 0;
153 14 : pszProjection = NULL;
154 14 : pszGCPProjection= NULL;
155 14 : pasGCPList = NULL;
156 14 : poH5Objects = NULL;
157 14 : poH5RootGroup = NULL;
158 14 : dims = NULL;
159 14 : maxdims = NULL;
160 14 : papszMetadata = NULL;
161 14 : adfGeoTransform[0] = 0.0;
162 14 : adfGeoTransform[1] = 1.0;
163 14 : adfGeoTransform[2] = 0.0;
164 14 : adfGeoTransform[3] = 0.0;
165 14 : adfGeoTransform[4] = 0.0;
166 14 : adfGeoTransform[5] = 1.0;
167 14 : iSubdatasetType = UNKNOWN_PRODUCT;
168 14 : bHasGeoTransform = false;
169 14 : }
170 :
171 : /************************************************************************/
172 : /* ~HDF5ImageDataset() */
173 : /************************************************************************/
174 14 : HDF5ImageDataset::~HDF5ImageDataset( )
175 : {
176 14 : FlushCache();
177 :
178 14 : CPLFree(pszProjection);
179 14 : CPLFree(pszGCPProjection);
180 :
181 14 : if( dims )
182 14 : CPLFree( dims );
183 :
184 14 : if( maxdims )
185 14 : CPLFree( maxdims );
186 :
187 14 : if( nGCPCount > 0 )
188 : {
189 5 : for( int i = 0; i < nGCPCount; i++ )
190 : {
191 4 : if( pasGCPList[i].pszId )
192 4 : CPLFree( pasGCPList[i].pszId );
193 4 : if( pasGCPList[i].pszInfo )
194 4 : CPLFree( pasGCPList[i].pszInfo );
195 : }
196 :
197 1 : CPLFree( pasGCPList );
198 : }
199 14 : }
200 :
201 : /************************************************************************/
202 : /* ==================================================================== */
203 : /* Hdf5imagerasterband */
204 : /* ==================================================================== */
205 : /************************************************************************/
206 : class HDF5ImageRasterBand : public GDALPamRasterBand
207 : {
208 : friend class HDF5ImageDataset;
209 :
210 : int bNoDataSet;
211 : double dfNoDataValue;
212 :
213 : public:
214 :
215 : HDF5ImageRasterBand( HDF5ImageDataset *, int, GDALDataType );
216 : ~HDF5ImageRasterBand();
217 :
218 : virtual CPLErr IReadBlock( int, int, void * );
219 : virtual double GetNoDataValue( int * );
220 : virtual CPLErr SetNoDataValue( double );
221 : /* virtual CPLErr IWriteBlock( int, int, void * ); */
222 : };
223 :
224 : /************************************************************************/
225 : /* ~HDF5ImageRasterBand() */
226 : /************************************************************************/
227 :
228 14 : HDF5ImageRasterBand::~HDF5ImageRasterBand()
229 : {
230 :
231 14 : }
232 : /************************************************************************/
233 : /* HDF5ImageRasterBand() */
234 : /************************************************************************/
235 14 : HDF5ImageRasterBand::HDF5ImageRasterBand( HDF5ImageDataset *poDS, int nBand,
236 14 : GDALDataType eType )
237 :
238 : {
239 : char **papszMetaGlobal;
240 14 : this->poDS = poDS;
241 14 : this->nBand = nBand;
242 14 : eDataType = eType;
243 14 : bNoDataSet = FALSE;
244 14 : dfNoDataValue = -9999;
245 14 : nBlockXSize = poDS->GetRasterXSize( );
246 14 : nBlockYSize = 1;
247 :
248 : /* -------------------------------------------------------------------- */
249 : /* Take a copy of Global Metadata since I can't pass Raster */
250 : /* variable to Iterate function. */
251 : /* -------------------------------------------------------------------- */
252 14 : papszMetaGlobal = CSLDuplicate( poDS->papszMetadata );
253 14 : CSLDestroy( poDS->papszMetadata );
254 14 : poDS->papszMetadata = NULL;
255 :
256 14 : if( poDS->poH5Objects->nType == H5G_DATASET ) {
257 14 : poDS->CreateMetadata( poDS->poH5Objects, H5G_DATASET );
258 : }
259 :
260 : /* -------------------------------------------------------------------- */
261 : /* Recover Global Metadat and set Band Metadata */
262 : /* -------------------------------------------------------------------- */
263 :
264 14 : SetMetadata( poDS->papszMetadata );
265 :
266 14 : CSLDestroy( poDS->papszMetadata );
267 14 : poDS->papszMetadata = CSLDuplicate( papszMetaGlobal );
268 14 : CSLDestroy( papszMetaGlobal );
269 :
270 : /* check for chunksize and set it as the blocksize (optimizes read) */
271 14 : hid_t listid = H5Dget_create_plist(((HDF5ImageDataset * )poDS)->dataset_id);
272 14 : if (listid>0)
273 : {
274 14 : if(H5Pget_layout(listid) == H5D_CHUNKED)
275 : {
276 : hsize_t panChunkDims[3];
277 4 : int nDimSize = H5Pget_chunk(listid, 3, panChunkDims);
278 4 : nBlockXSize = (int) panChunkDims[nDimSize-1];
279 4 : nBlockYSize = (int) panChunkDims[nDimSize-2];
280 : }
281 14 : H5Pclose(listid);
282 : }
283 :
284 14 : }
285 :
286 : /************************************************************************/
287 : /* GetNoDataValue() */
288 : /************************************************************************/
289 3 : double HDF5ImageRasterBand::GetNoDataValue( int * pbSuccess )
290 :
291 : {
292 3 : if( bNoDataSet )
293 : {
294 0 : if( pbSuccess )
295 0 : *pbSuccess = bNoDataSet;
296 :
297 0 : return dfNoDataValue;
298 : }
299 : else
300 3 : return GDALPamRasterBand::GetNoDataValue( pbSuccess );
301 : }
302 :
303 : /************************************************************************/
304 : /* SetNoDataValue() */
305 : /************************************************************************/
306 0 : CPLErr HDF5ImageRasterBand::SetNoDataValue( double dfNoData )
307 :
308 : {
309 0 : bNoDataSet = TRUE;
310 0 : dfNoDataValue = dfNoData;
311 :
312 0 : return CE_None;
313 : }
314 :
315 : /************************************************************************/
316 : /* IReadBlock() */
317 : /************************************************************************/
318 919 : CPLErr HDF5ImageRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
319 : void * pImage )
320 : {
321 : herr_t status;
322 : hsize_t count[3];
323 : H5OFFSET_TYPE offset[3];
324 : int nSizeOfData;
325 : hid_t memspace;
326 : hsize_t col_dims[3];
327 : hsize_t rank;
328 :
329 919 : HDF5ImageDataset *poGDS = ( HDF5ImageDataset * ) poDS;
330 :
331 919 : if( poGDS->eAccess == GA_Update ) {
332 : memset( pImage, 0,
333 : nBlockXSize * nBlockYSize *
334 0 : GDALGetDataTypeSize( eDataType )/8 );
335 0 : return CE_None;
336 : }
337 :
338 919 : rank=2;
339 :
340 919 : if( poGDS->ndims == 3 ){
341 0 : rank=3;
342 0 : offset[0] = nBand-1;
343 0 : count[0] = 1;
344 0 : col_dims[0] = 1;
345 : }
346 :
347 919 : offset[poGDS->ndims - 2] = nBlockYOff*nBlockYSize;
348 919 : offset[poGDS->ndims - 1] = nBlockXOff*nBlockXSize;
349 919 : count[poGDS->ndims - 2] = nBlockYSize;
350 919 : count[poGDS->ndims - 1] = nBlockXSize;
351 :
352 919 : nSizeOfData = H5Tget_size( poGDS->native );
353 919 : memset( pImage,0,nBlockXSize*nBlockYSize*nSizeOfData );
354 :
355 : /* blocksize may not be a multiple of imagesize */
356 1845 : count[poGDS->ndims - 2] = MIN( size_t(nBlockYSize),
357 : poDS->GetRasterYSize() -
358 1845 : offset[poGDS->ndims - 2]);
359 2757 : count[poGDS->ndims - 1] = MIN( size_t(nBlockXSize),
360 : poDS->GetRasterXSize()-
361 2757 : offset[poGDS->ndims - 1]);
362 :
363 : /* -------------------------------------------------------------------- */
364 : /* Select block from file space */
365 : /* -------------------------------------------------------------------- */
366 : status = H5Sselect_hyperslab( poGDS->dataspace_id,
367 : H5S_SELECT_SET,
368 : offset, NULL,
369 919 : count, NULL );
370 :
371 : /* -------------------------------------------------------------------- */
372 : /* Create memory space to receive the data */
373 : /* -------------------------------------------------------------------- */
374 919 : col_dims[poGDS->ndims-2]=nBlockYSize;
375 919 : col_dims[poGDS->ndims-1]=nBlockXSize;
376 919 : memspace = H5Screate_simple( (int) rank, col_dims, NULL );
377 919 : H5OFFSET_TYPE mem_offset[3] = {0, 0, 0};
378 : status = H5Sselect_hyperslab(memspace,
379 : H5S_SELECT_SET,
380 : mem_offset, NULL,
381 919 : count, NULL);
382 :
383 : status = H5Dread ( poGDS->dataset_id,
384 : poGDS->native,
385 : memspace,
386 : poGDS->dataspace_id,
387 : H5P_DEFAULT,
388 919 : pImage );
389 :
390 919 : H5Sclose( memspace );
391 :
392 919 : if( status < 0 )
393 : {
394 : CPLError( CE_Failure, CPLE_AppDefined,
395 0 : "H5Dread() failed for block." );
396 0 : return CE_Failure;
397 : }
398 : else
399 919 : return CE_None;
400 : }
401 :
402 : /************************************************************************/
403 : /* Identify() */
404 : /************************************************************************/
405 :
406 10129 : int HDF5ImageDataset::Identify( GDALOpenInfo *poOpenInfo )
407 :
408 : {
409 10129 : if(!EQUALN( poOpenInfo->pszFilename, "HDF5:", 5 ) )
410 10129 : return FALSE;
411 : else
412 0 : return TRUE;
413 : }
414 :
415 : /************************************************************************/
416 : /* Open() */
417 : /************************************************************************/
418 1733 : GDALDataset *HDF5ImageDataset::Open( GDALOpenInfo * poOpenInfo )
419 : {
420 : int i;
421 : HDF5ImageDataset *poDS;
422 : char szFilename[2048];
423 :
424 1733 : if(!EQUALN( poOpenInfo->pszFilename, "HDF5:", 5 ) ||
425 : strlen(poOpenInfo->pszFilename) > sizeof(szFilename) - 3 )
426 1719 : return NULL;
427 :
428 : /* -------------------------------------------------------------------- */
429 : /* Confirm the requested access is supported. */
430 : /* -------------------------------------------------------------------- */
431 14 : if( poOpenInfo->eAccess == GA_Update )
432 : {
433 : CPLError( CE_Failure, CPLE_NotSupported,
434 : "The HDF5ImageDataset driver does not support update access to existing"
435 0 : " datasets.\n" );
436 0 : return NULL;
437 : }
438 :
439 14 : poDS = new HDF5ImageDataset();
440 :
441 : /* -------------------------------------------------------------------- */
442 : /* Create a corresponding GDALDataset. */
443 : /* -------------------------------------------------------------------- */
444 : /* printf("poOpenInfo->pszFilename %s\n",poOpenInfo->pszFilename); */
445 : char **papszName =
446 : CSLTokenizeString2( poOpenInfo->pszFilename,
447 14 : ":", CSLT_HONOURSTRINGS|CSLT_PRESERVEESCAPES );
448 :
449 14 : if( !((CSLCount(papszName) == 3) || (CSLCount(papszName) == 4)) )
450 : {
451 0 : CSLDestroy(papszName);
452 0 : delete poDS;
453 0 : return NULL;
454 : }
455 :
456 14 : poDS->SetDescription( poOpenInfo->pszFilename );
457 :
458 : /* -------------------------------------------------------------------- */
459 : /* Check for drive name in windows HDF5:"D:\... */
460 : /* -------------------------------------------------------------------- */
461 14 : strcpy(szFilename, papszName[1]);
462 :
463 14 : if( strlen(papszName[1]) == 1 && papszName[3] != NULL )
464 : {
465 0 : strcat(szFilename, ":");
466 0 : strcat(szFilename, papszName[2]);
467 :
468 0 : poDS->SetSubdatasetName( papszName[3] );
469 : }
470 : else
471 14 : poDS->SetSubdatasetName( papszName[2] );
472 :
473 14 : CSLDestroy(papszName);
474 14 : papszName = NULL;
475 :
476 14 : if( !H5Fis_hdf5(szFilename) ) {
477 0 : delete poDS;
478 0 : return NULL;
479 : }
480 :
481 14 : poDS->SetPhysicalFilename( szFilename );
482 :
483 : /* -------------------------------------------------------------------- */
484 : /* Try opening the dataset. */
485 : /* -------------------------------------------------------------------- */
486 : poDS->hHDF5 = H5Fopen(szFilename,
487 : H5F_ACC_RDONLY,
488 14 : H5P_DEFAULT );
489 :
490 14 : if( poDS->hHDF5 < 0 )
491 : {
492 0 : delete poDS;
493 0 : return NULL;
494 : }
495 :
496 14 : poDS->hGroupID = H5Gopen( poDS->hHDF5, "/" );
497 14 : if( poDS->hGroupID < 0 )
498 : {
499 0 : poDS->bIsHDFEOS=false;
500 0 : delete poDS;
501 0 : return NULL;
502 : }
503 :
504 : /* -------------------------------------------------------------------- */
505 : /* THIS IS AN HDF5 FILE */
506 : /* -------------------------------------------------------------------- */
507 14 : poDS->bIsHDFEOS=TRUE;
508 14 : poDS->ReadGlobalAttributes( FALSE );
509 :
510 : /* -------------------------------------------------------------------- */
511 : /* Create HDF5 Data Hierarchy in a link list */
512 : /* -------------------------------------------------------------------- */
513 : poDS->poH5Objects =
514 : poDS->HDF5FindDatasetObjectsbyPath( poDS->poH5RootGroup,
515 14 : (char *)poDS->GetSubdatasetName() );
516 :
517 14 : if( poDS->poH5Objects == NULL ) {
518 0 : delete poDS;
519 0 : return NULL;
520 : }
521 : /* -------------------------------------------------------------------- */
522 : /* Retrieve HDF5 data information */
523 : /* -------------------------------------------------------------------- */
524 14 : poDS->dataset_id = H5Dopen( poDS->hHDF5,poDS->poH5Objects->pszPath );
525 14 : poDS->dataspace_id = H5Dget_space( poDS->dataset_id );
526 14 : poDS->ndims = H5Sget_simple_extent_ndims( poDS->dataspace_id );
527 14 : poDS->dims = (hsize_t*)CPLCalloc( poDS->ndims, sizeof(hsize_t) );
528 14 : poDS->maxdims = (hsize_t*)CPLCalloc( poDS->ndims, sizeof(hsize_t) );
529 : poDS->dimensions = H5Sget_simple_extent_dims( poDS->dataspace_id,
530 : poDS->dims,
531 14 : poDS->maxdims );
532 14 : poDS->datatype = H5Dget_type( poDS->dataset_id );
533 14 : poDS->clas = H5Tget_class( poDS->datatype );
534 14 : poDS->size = H5Tget_size( poDS->datatype );
535 14 : poDS->address = H5Dget_offset( poDS->dataset_id );
536 14 : poDS->native = H5Tget_native_type( poDS->datatype, H5T_DIR_ASCEND );
537 :
538 14 : poDS->nRasterYSize=(int)poDS->dims[poDS->ndims-2]; // Y
539 14 : poDS->nRasterXSize=(int)poDS->dims[poDS->ndims-1]; // X alway last
540 :
541 14 : poDS->nBands=1;
542 :
543 14 : if( poDS->ndims == 3 ) poDS->nBands=(int) poDS->dims[0];
544 :
545 :
546 28 : for( i = 1; i <= poDS->nBands; i++ ) {
547 : HDF5ImageRasterBand *poBand =
548 : new HDF5ImageRasterBand( poDS, i,
549 14 : poDS->GetDataType( poDS->native ) );
550 :
551 14 : poDS->SetBand( i, poBand );
552 14 : if( poBand->bNoDataSet )
553 0 : poBand->SetNoDataValue( 255 );
554 : }
555 :
556 : // CSK code in IdentifyProductType() and CreateProjections()
557 : // uses dataset metadata.
558 14 : poDS->SetMetadata( poDS->papszMetadata );
559 :
560 : // Check if the hdf5 is a well known product type
561 14 : poDS->IdentifyProductType();
562 :
563 14 : poDS->CreateProjections( );
564 :
565 : /* -------------------------------------------------------------------- */
566 : /* Setup/check for pam .aux.xml. */
567 : /* -------------------------------------------------------------------- */
568 14 : poDS->TryLoadXML();
569 :
570 : /* -------------------------------------------------------------------- */
571 : /* Setup overviews. */
572 : /* -------------------------------------------------------------------- */
573 14 : poDS->oOvManager.Initialize( poDS, ":::VIRTUAL:::" );
574 :
575 14 : return( poDS );
576 : }
577 :
578 :
579 : /************************************************************************/
580 : /* GDALRegister_HDF5Image() */
581 : /************************************************************************/
582 610 : void GDALRegister_HDF5Image( )
583 :
584 : {
585 : GDALDriver *poDriver;
586 :
587 610 : if (! GDAL_CHECK_VERSION("HDF5Image driver"))
588 0 : return;
589 :
590 610 : if( GDALGetDriverByName( "HDF5Image" ) == NULL )
591 : {
592 588 : poDriver = new GDALDriver( );
593 :
594 588 : poDriver->SetDescription( "HDF5Image" );
595 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
596 588 : "HDF5 Dataset" );
597 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
598 588 : "frmt_hdf5.html" );
599 588 : poDriver->pfnOpen = HDF5ImageDataset::Open;
600 588 : poDriver->pfnIdentify = HDF5ImageDataset::Identify;
601 :
602 588 : GetGDALDriverManager( )->RegisterDriver( poDriver );
603 : }
604 : }
605 :
606 : /************************************************************************/
607 : /* CreateProjections() */
608 : /************************************************************************/
609 14 : CPLErr HDF5ImageDataset::CreateProjections()
610 : {
611 14 : switch(iSubdatasetType)
612 : {
613 : case CSK_PRODUCT:
614 : {
615 2 : const char *osMissionLevel = NULL;
616 2 : int productType = PROD_UNKNOWN;
617 :
618 2 : if(GetMetadataItem("Product_Type")!=NULL)
619 : {
620 : //Get the format's level
621 2 : osMissionLevel = HDF5Dataset::GetMetadataItem("Product_Type");
622 :
623 2 : if(EQUALN(osMissionLevel,"RAW",3))
624 0 : productType = PROD_CSK_L0;
625 :
626 2 : if(EQUALN(osMissionLevel,"SSC",3))
627 0 : productType = PROD_CSK_L1A;
628 :
629 2 : if(EQUALN(osMissionLevel,"DGM",3))
630 1 : productType = PROD_CSK_L1B;
631 :
632 2 : if(EQUALN(osMissionLevel,"GEC",3))
633 1 : productType = PROD_CSK_L1C;
634 :
635 2 : if(EQUALN(osMissionLevel,"GTC",3))
636 0 : productType = PROD_CSK_L1D;
637 : }
638 :
639 2 : CaptureCSKGeoTransform(productType);
640 2 : CaptureCSKGeolocation(productType);
641 2 : CaptureCSKGCPs(productType);
642 :
643 2 : break;
644 : }
645 : case UNKNOWN_PRODUCT:
646 : {
647 : #define NBGCPLAT 100
648 : #define NBGCPLON 30
649 :
650 12 : hid_t LatitudeDatasetID = -1;
651 12 : hid_t LongitudeDatasetID = -1;
652 : hid_t LatitudeDataspaceID;
653 : hid_t LongitudeDataspaceID;
654 : float* Latitude;
655 : float* Longitude;
656 : int i,j;
657 : int nDeltaLat;
658 : int nDeltaLon;
659 :
660 12 : nDeltaLat = nRasterYSize / NBGCPLAT;
661 12 : nDeltaLon = nRasterXSize / NBGCPLON;
662 :
663 12 : if( nDeltaLat == 0 || nDeltaLon == 0 )
664 10 : return CE_None;
665 :
666 : /* -------------------------------------------------------------------- */
667 : /* Create HDF5 Data Hierarchy in a link list */
668 : /* -------------------------------------------------------------------- */
669 2 : poH5Objects=HDF5FindDatasetObjects( poH5RootGroup, "Latitude" );
670 2 : if( !poH5Objects ) {
671 2 : return CE_None;
672 : }
673 : /* -------------------------------------------------------------------- */
674 : /* The Lattitude and Longitude arrays must have a rank of 2 to */
675 : /* retrieve GCPs. */
676 : /* -------------------------------------------------------------------- */
677 0 : if( poH5Objects->nRank != 2 ) {
678 0 : return CE_None;
679 : }
680 :
681 : /* -------------------------------------------------------------------- */
682 : /* Retrieve HDF5 data information */
683 : /* -------------------------------------------------------------------- */
684 0 : LatitudeDatasetID = H5Dopen( hHDF5,poH5Objects->pszPath );
685 0 : LatitudeDataspaceID = H5Dget_space( dataset_id );
686 :
687 0 : poH5Objects=HDF5FindDatasetObjects( poH5RootGroup, "Longitude" );
688 0 : LongitudeDatasetID = H5Dopen( hHDF5,poH5Objects->pszPath );
689 0 : LongitudeDataspaceID = H5Dget_space( dataset_id );
690 :
691 0 : if( ( LatitudeDatasetID > 0 ) && ( LongitudeDatasetID > 0) ) {
692 :
693 : Latitude = ( float * ) CPLCalloc( nRasterYSize*nRasterXSize,
694 0 : sizeof( float ) );
695 : Longitude = ( float * ) CPLCalloc( nRasterYSize*nRasterXSize,
696 0 : sizeof( float ) );
697 0 : memset( Latitude, 0, nRasterXSize*nRasterYSize*sizeof( float ) );
698 0 : memset( Longitude, 0, nRasterXSize*nRasterYSize*sizeof( float ) );
699 :
700 : H5Dread ( LatitudeDatasetID,
701 : H5T_NATIVE_FLOAT,
702 : H5S_ALL,
703 : H5S_ALL,
704 : H5P_DEFAULT,
705 0 : Latitude );
706 :
707 : H5Dread ( LongitudeDatasetID,
708 : H5T_NATIVE_FLOAT,
709 : H5S_ALL,
710 : H5S_ALL,
711 : H5P_DEFAULT,
712 0 : Longitude );
713 :
714 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
715 0 : CPLFree(pszProjection);
716 0 : CPLFree(pszGCPProjection);
717 0 : oSRS.exportToWkt( &pszProjection );
718 0 : oSRS.exportToWkt( &pszGCPProjection );
719 :
720 : /* -------------------------------------------------------------------- */
721 : /* Fill the GCPs list. */
722 : /* -------------------------------------------------------------------- */
723 0 : nGCPCount = nRasterYSize/nDeltaLat * nRasterXSize/nDeltaLon;
724 :
725 : pasGCPList = ( GDAL_GCP * )
726 0 : CPLCalloc( nGCPCount, sizeof( GDAL_GCP ) );
727 :
728 0 : GDALInitGCPs( nGCPCount, pasGCPList );
729 0 : int k=0;
730 :
731 0 : int nYLimit = ((int)nRasterYSize/nDeltaLat) * nDeltaLat;
732 0 : int nXLimit = ((int)nRasterXSize/nDeltaLon) * nDeltaLon;
733 0 : for( j = 0; j < nYLimit; j+=nDeltaLat )
734 : {
735 0 : for( i = 0; i < nXLimit; i+=nDeltaLon )
736 : {
737 0 : int iGCP = j * nRasterXSize + i;
738 0 : pasGCPList[k].dfGCPX = ( double ) Longitude[iGCP]+180.0;
739 0 : pasGCPList[k].dfGCPY = ( double ) Latitude[iGCP];
740 :
741 0 : pasGCPList[k].dfGCPPixel = i + 0.5;
742 0 : pasGCPList[k++].dfGCPLine = j + 0.5;
743 :
744 : }
745 : }
746 :
747 0 : CPLFree( Latitude );
748 0 : CPLFree( Longitude );
749 : }
750 : break;
751 : }
752 : }
753 2 : return CE_None;
754 :
755 : }
756 : /************************************************************************/
757 : /* GetProjectionRef() */
758 : /************************************************************************/
759 :
760 1 : const char *HDF5ImageDataset::GetProjectionRef( )
761 :
762 : {
763 1 : if( pszProjection )
764 1 : return pszProjection;
765 : else
766 0 : return GDALPamDataset::GetProjectionRef();
767 : }
768 :
769 : /************************************************************************/
770 : /* GetGCPCount() */
771 : /************************************************************************/
772 :
773 1 : int HDF5ImageDataset::GetGCPCount( )
774 :
775 : {
776 1 : if( nGCPCount > 0 )
777 1 : return nGCPCount;
778 : else
779 0 : return GDALPamDataset::GetGCPCount();
780 : }
781 :
782 : /************************************************************************/
783 : /* GetGCPProjection() */
784 : /************************************************************************/
785 :
786 1 : const char *HDF5ImageDataset::GetGCPProjection( )
787 :
788 : {
789 1 : if( nGCPCount > 0 )
790 1 : return pszGCPProjection;
791 : else
792 0 : return GDALPamDataset::GetGCPProjection();
793 : }
794 :
795 : /************************************************************************/
796 : /* GetGCPs() */
797 : /************************************************************************/
798 :
799 1 : const GDAL_GCP *HDF5ImageDataset::GetGCPs( )
800 : {
801 1 : if( nGCPCount > 0 )
802 1 : return pasGCPList;
803 : else
804 0 : return GDALPamDataset::GetGCPs();
805 : }
806 :
807 : /************************************************************************/
808 : /* GetGeoTransform() */
809 : /************************************************************************/
810 :
811 1 : CPLErr HDF5ImageDataset::GetGeoTransform( double * padfTransform )
812 : {
813 1 : if ( bHasGeoTransform )
814 : {
815 1 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
816 1 : return CE_None;
817 : }
818 : else
819 0 : return GDALPamDataset::GetGeoTransform(padfTransform);
820 : }
821 :
822 : /************************************************************************/
823 : /* IdentifyProductType() */
824 : /************************************************************************/
825 :
826 : /**
827 : * Identify if the subdataset has a known product format
828 : * It stores a product identifier in iSubdatasetType,
829 : * UNKNOWN_PRODUCT, if it isn't a recognizable format.
830 : */
831 14 : void HDF5ImageDataset::IdentifyProductType()
832 : {
833 14 : iSubdatasetType = UNKNOWN_PRODUCT;
834 :
835 : /************************************************************************/
836 : /* COSMO-SKYMED */
837 : /************************************************************************/
838 : const char *pszMissionId;
839 :
840 : //Get the Mission Id as a char *, because the
841 : //field may not exist
842 14 : pszMissionId = HDF5Dataset::GetMetadataItem("Mission_ID");
843 :
844 : //If there is a Mission_ID field
845 14 : if(pszMissionId != NULL && strstr(GetDescription(), "QLK") == NULL)
846 : //Check if the mission type is CSK
847 2 : if(EQUAL(pszMissionId,"CSK"))
848 2 : iSubdatasetType = CSK_PRODUCT;
849 14 : }
850 :
851 : /************************************************************************/
852 : /* CaptureCSKGeolocation() */
853 : /************************************************************************/
854 :
855 : /**
856 : * Captures Geolocation information from a COSMO-SKYMED
857 : * file.
858 : * The geoid will allways be WGS84
859 : * The projection type may be UTM or UPS, depending on the
860 : * latitude from the center of the image.
861 : * @param iProductType type of CSK subproduct, see HDF5CSKProduct
862 : */
863 2 : void HDF5ImageDataset::CaptureCSKGeolocation(int iProductType)
864 : {
865 : double *dfProjFalseEastNorth;
866 : double *dfProjScaleFactor;
867 : double *dfCenterCoord;
868 :
869 3 : if(iProductType == PROD_CSK_L1C||iProductType == PROD_CSK_L1D)
870 : {
871 : //Check if all the metadata attributes are present
872 2 : if(HDF5ReadDoubleAttr("Map Projection False East-North", &dfProjFalseEastNorth) == CE_Failure||
873 : HDF5ReadDoubleAttr("Map Projection Scale Factor", &dfProjScaleFactor) == CE_Failure||
874 : HDF5ReadDoubleAttr("Map Projection Centre", &dfCenterCoord) == CE_Failure||
875 1 : GetMetadataItem("Projection_ID") == NULL)
876 : {
877 :
878 0 : pszProjection = CPLStrdup("");
879 0 : pszGCPProjection = CPLStrdup("");
880 : CPLError( CE_Failure, CPLE_OpenFailed,
881 0 : "The CSK hdf5 file geolocation information is malformed\n" );
882 : }
883 : else
884 : {
885 :
886 : //Fetch projection Type
887 1 : CPLString osProjectionID = GetMetadataItem("Projection_ID");
888 :
889 : //If the projection is UTM
890 1 : if(EQUAL(osProjectionID,"UTM"))
891 : {
892 : // @TODO: use SetUTM
893 1 : oSRS.SetProjCS(SRS_PT_TRANSVERSE_MERCATOR);
894 : oSRS.SetTM(dfCenterCoord[0],
895 : dfCenterCoord[1],
896 : dfProjScaleFactor[0],
897 : dfProjFalseEastNorth[0],
898 1 : dfProjFalseEastNorth[1]);
899 : }
900 : else
901 : {
902 : //TODO Test! I didn't had any UPS projected files to test!
903 : //If the projection is UPS
904 0 : if(EQUAL(osProjectionID,"UPS"))
905 : {
906 0 : oSRS.SetProjCS(SRS_PT_POLAR_STEREOGRAPHIC);
907 : oSRS.SetPS(dfCenterCoord[0],
908 : dfCenterCoord[1],
909 : dfProjScaleFactor[0],
910 : dfProjFalseEastNorth[0],
911 0 : dfProjFalseEastNorth[1]);
912 : }
913 : }
914 :
915 : //Export Projection to Wkt.
916 : //In case of error then clean the projection
917 1 : if (oSRS.exportToWkt(&pszProjection) != OGRERR_NONE)
918 0 : pszProjection = CPLStrdup("");
919 :
920 1 : CPLFree(dfCenterCoord);
921 1 : CPLFree(dfProjScaleFactor);
922 1 : CPLFree(dfProjFalseEastNorth);
923 : }
924 : }
925 : else
926 : {
927 : //Set the ellipsoid to WGS84
928 1 : oSRS.SetWellKnownGeogCS( "WGS84" );
929 :
930 : //Export GCPProjection to Wkt.
931 : //In case of error then clean the projection
932 1 : if(oSRS.exportToWkt(&pszGCPProjection) != OGRERR_NONE)
933 0 : pszGCPProjection = CPLStrdup("");
934 : }
935 2 : }
936 :
937 : /************************************************************************/
938 : /* CaptureCSKGeoTransform() */
939 : /************************************************************************/
940 :
941 : /**
942 : * Get Geotransform information for COSMO-SKYMED files
943 : * In case of sucess it stores the transformation
944 : * in adfGeoTransform. In case of failure it doesn't
945 : * modify adfGeoTransform
946 : * @param iProductType type of CSK subproduct, see HDF5CSKProduct
947 : */
948 2 : void HDF5ImageDataset::CaptureCSKGeoTransform(int iProductType)
949 : {
950 : double *pdOutUL;
951 : double *pdLineSpacing;
952 : double *pdColumnSpacing;
953 :
954 2 : CPLString osULCoord;
955 2 : CPLString osULPath;
956 2 : CPLString osLineSpacingPath, osColumnSpacingPath;
957 :
958 2 : const char *pszSubdatasetName = GetSubdatasetName();
959 :
960 2 : bHasGeoTransform = FALSE;
961 : //If the product level is not L1C or L1D then
962 : //it doesn't have a valid projection
963 2 : if(iProductType == PROD_CSK_L1C||iProductType == PROD_CSK_L1D)
964 : {
965 : //If there is a subdataset
966 1 : if(pszSubdatasetName != NULL)
967 : {
968 :
969 1 : osULPath = pszSubdatasetName ;
970 1 : osULPath += "/Top Left East-North";
971 :
972 1 : osLineSpacingPath = pszSubdatasetName;
973 1 : osLineSpacingPath += "/Line Spacing";
974 :
975 1 : osColumnSpacingPath = pszSubdatasetName;
976 1 : osColumnSpacingPath += "/Column Spacing";
977 :
978 :
979 : //If it could find the attributes on the metadata
980 1 : if(HDF5ReadDoubleAttr(osULPath.c_str(), &pdOutUL) == CE_Failure ||
981 : HDF5ReadDoubleAttr(osLineSpacingPath.c_str(), &pdLineSpacing) == CE_Failure ||
982 : HDF5ReadDoubleAttr(osColumnSpacingPath.c_str(), &pdColumnSpacing) == CE_Failure)
983 : {
984 0 : bHasGeoTransform = FALSE;
985 : }
986 : else
987 : {
988 : // geotransform[1] : width of pixel
989 : // geotransform[4] : rotational coefficient, zero for north up images.
990 : // geotransform[2] : rotational coefficient, zero for north up images.
991 : // geotransform[5] : height of pixel (but negative)
992 :
993 1 : adfGeoTransform[0] = pdOutUL[0];
994 1 : adfGeoTransform[1] = pdLineSpacing[0];
995 1 : adfGeoTransform[2] = 0;
996 1 : adfGeoTransform[3] = pdOutUL[1];
997 1 : adfGeoTransform[4] = 0;
998 1 : adfGeoTransform[5] = -pdColumnSpacing[0];
999 :
1000 :
1001 1 : CPLFree(pdOutUL);
1002 1 : CPLFree(pdLineSpacing);
1003 1 : CPLFree(pdColumnSpacing);
1004 :
1005 1 : bHasGeoTransform = TRUE;
1006 : }
1007 : }
1008 2 : }
1009 2 : }
1010 :
1011 :
1012 : /************************************************************************/
1013 : /* CaptureCSKGCPs() */
1014 : /************************************************************************/
1015 :
1016 : /**
1017 : * Retrieves and stores the GCPs from a COSMO-SKYMED dataset
1018 : * It only retrieves the GCPs for L0, L1A and L1B products
1019 : * for L1C and L1D products, geotransform is provided.
1020 : * The GCPs provided will be the Image's corners.
1021 : * @param iProductType type of CSK product @see HDF5CSKProductEnum
1022 : */
1023 2 : void HDF5ImageDataset::CaptureCSKGCPs(int iProductType)
1024 : {
1025 : //Only retrieve GCPs for L0,L1A and L1B products
1026 2 : if(iProductType == PROD_CSK_L0||iProductType == PROD_CSK_L1A||
1027 : iProductType == PROD_CSK_L1B)
1028 : {
1029 : int i;
1030 : double *pdCornerCoordinates;
1031 :
1032 1 : nGCPCount=4;
1033 1 : pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),4);
1034 1 : CPLString osCornerName[4];
1035 : double pdCornerPixel[4];
1036 : double pdCornerLine[4];
1037 :
1038 2 : const char *pszSubdatasetName = GetSubdatasetName();
1039 :
1040 : //Load the subdataset name first
1041 5 : for(i=0;i <4;i++)
1042 4 : osCornerName[i] = pszSubdatasetName;
1043 :
1044 : //Load the attribute name, and raster coordinates for
1045 : //all the corners
1046 1 : osCornerName[0] += "/Top Left Geodetic Coordinates";
1047 1 : pdCornerPixel[0] = 0;
1048 1 : pdCornerLine[0] = 0;
1049 :
1050 1 : osCornerName[1] += "/Top Right Geodetic Coordinates";
1051 1 : pdCornerPixel[1] = GetRasterXSize();
1052 1 : pdCornerLine[1] = 0;
1053 :
1054 1 : osCornerName[2] += "/Bottom Left Geodetic Coordinates";
1055 1 : pdCornerPixel[2] = 0;
1056 1 : pdCornerLine[2] = GetRasterYSize();
1057 :
1058 1 : osCornerName[3] += "/Bottom Right Geodetic Coordinates";
1059 1 : pdCornerPixel[3] = GetRasterXSize();
1060 1 : pdCornerLine[3] = GetRasterYSize();
1061 :
1062 : //For all the image's corners
1063 5 : for(i=0;i<4;i++)
1064 : {
1065 4 : GDALInitGCPs( 1, pasGCPList + i );
1066 :
1067 4 : CPLFree( pasGCPList[i].pszId );
1068 4 : pasGCPList[i].pszId = NULL;
1069 :
1070 : //Retrieve the attributes
1071 4 : if(HDF5ReadDoubleAttr(osCornerName[i].c_str(),
1072 : &pdCornerCoordinates) == CE_Failure)
1073 : {
1074 : CPLError( CE_Failure, CPLE_OpenFailed,
1075 0 : "Error retrieving CSK GCPs\n" );
1076 : // Free on failure, e.g. in case of QLK subdataset.
1077 0 : for( int i = 0; i < 4; i++ )
1078 : {
1079 0 : if( pasGCPList[i].pszId )
1080 0 : CPLFree( pasGCPList[i].pszId );
1081 0 : if( pasGCPList[i].pszInfo )
1082 0 : CPLFree( pasGCPList[i].pszInfo );
1083 : }
1084 0 : CPLFree( pasGCPList );
1085 0 : pasGCPList = NULL;
1086 0 : nGCPCount = 0;
1087 0 : break;
1088 : }
1089 :
1090 : //Fill the GCPs name
1091 4 : pasGCPList[i].pszId = CPLStrdup( osCornerName[i].c_str() );
1092 :
1093 : //Fill the coordinates
1094 4 : pasGCPList[i].dfGCPX = pdCornerCoordinates[1];
1095 4 : pasGCPList[i].dfGCPY = pdCornerCoordinates[0];
1096 4 : pasGCPList[i].dfGCPZ = pdCornerCoordinates[2];
1097 4 : pasGCPList[i].dfGCPPixel = pdCornerPixel[i];
1098 4 : pasGCPList[i].dfGCPLine = pdCornerLine[i];
1099 :
1100 : //Free the returned coordinates
1101 4 : CPLFree(pdCornerCoordinates);
1102 1 : }
1103 : }
1104 2 : }
1105 :
|
