1 : /******************************************************************************
2 : * $Id: hdf4imagedataset.cpp 25607 2013-02-06 20:00:43Z warmerdam $
3 : *
4 : * Project: Hierarchical Data Format Release 4 (HDF4)
5 : * Purpose: Read subdatasets of HDF4 file.
6 : * This driver initially based on code supplied by Markus Neteler
7 : * Author: Andrey Kiselev, dron@ak4719.spb.edu
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2002, Andrey Kiselev <dron@ak4719.spb.edu>
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ****************************************************************************/
30 :
31 : #include <string.h>
32 : #include <math.h>
33 : #include "cpl_multiproc.h"
34 :
35 : #include "hdf.h"
36 : #include "mfhdf.h"
37 :
38 : #include "HdfEosDef.h"
39 :
40 : #include "gdal_priv.h"
41 : #include "cpl_string.h"
42 : #include "ogr_spatialref.h"
43 :
44 : #include "hdf4compat.h"
45 : #include "hdf4dataset.h"
46 :
47 : #include "nasakeywordhandler.h"
48 :
49 : CPL_CVSID("$Id: hdf4imagedataset.cpp 25607 2013-02-06 20:00:43Z warmerdam $");
50 :
51 : CPL_C_START
52 : void GDALRegister_HDF4(void);
53 : CPL_C_END
54 :
55 : #define HDF4_SDS_MAXNAMELEN 65
56 :
57 : // Signature to recognize files written by GDAL
58 : const char *pszGDALSignature =
59 : "Created with GDAL (http://www.remotesensing.org/gdal/)";
60 :
61 : extern void *hHDF4Mutex;
62 :
63 : /************************************************************************/
64 : /* ==================================================================== */
65 : /* List of HDF-EOS Swath product types. */
66 : /* ==================================================================== */
67 : /************************************************************************/
68 :
69 : enum HDF4EOSProduct
70 : {
71 : PROD_UNKNOWN,
72 : PROD_ASTER_L1A,
73 : PROD_ASTER_L1B,
74 : PROD_ASTER_L2,
75 : PROD_ASTER_L3,
76 : PROD_AST14DEM,
77 : PROD_MODIS_L1B,
78 : PROD_MODIS_L2
79 : };
80 :
81 : /************************************************************************/
82 : /* ==================================================================== */
83 : /* HDF4ImageDataset */
84 : /* ==================================================================== */
85 : /************************************************************************/
86 :
87 : class HDF4ImageDataset : public HDF4Dataset
88 : {
89 : friend class HDF4ImageRasterBand;
90 :
91 : char *pszFilename;
92 : int32 hHDF4, iGR, iPal, iDataset;
93 : int32 iRank, iNumType, nAttrs,
94 : iInterlaceMode, iPalInterlaceMode, iPalDataType;
95 : int32 nComps, nPalEntries;
96 : int32 aiDimSizes[H4_MAX_VAR_DIMS];
97 : int iXDim, iYDim, iBandDim, i4Dim;
98 : int nBandCount;
99 : char **papszLocalMetadata;
100 : #define N_COLOR_ENTRIES 256
101 : uint8 aiPaletteData[N_COLOR_ENTRIES][3]; // XXX: Static array for now
102 : char szName[HDF4_SDS_MAXNAMELEN];
103 : char *pszSubdatasetName;
104 : char *pszFieldName;
105 :
106 : GDALColorTable *poColorTable;
107 :
108 : OGRSpatialReference oSRS;
109 : int bHasGeoTransform;
110 : double adfGeoTransform[6];
111 : char *pszProjection;
112 : char *pszGCPProjection;
113 : GDAL_GCP *pasGCPList;
114 : int nGCPCount;
115 :
116 : HDF4DatasetType iDatasetType;
117 :
118 : int32 iSDS;
119 :
120 : int nBlockPreferredXSize;
121 : int nBlockPreferredYSize;
122 : bool bReadTile;
123 :
124 : void ToGeoref( double *, double * );
125 : void GetImageDimensions( char * );
126 : void GetSwatAttrs( int32 );
127 : void GetGridAttrs( int32 hGD );
128 : void CaptureNRLGeoTransform(void);
129 : void CaptureL1GMTLInfo(void);
130 : void CaptureCoastwatchGCTPInfo(void);
131 : void ProcessModisSDSGeolocation(void);
132 : int ProcessSwathGeolocation( int32, char ** );
133 :
134 : static long USGSMnemonicToCode( const char* );
135 : static void ReadCoordinates( const char*, double*, double* );
136 :
137 : public:
138 : HDF4ImageDataset();
139 : ~HDF4ImageDataset();
140 :
141 : static GDALDataset *Open( GDALOpenInfo * );
142 : static GDALDataset *Create( const char * pszFilename,
143 : int nXSize, int nYSize, int nBands,
144 : GDALDataType eType, char ** papszParmList );
145 : virtual void FlushCache( void );
146 : CPLErr GetGeoTransform( double * padfTransform );
147 : virtual CPLErr SetGeoTransform( double * );
148 : const char *GetProjectionRef();
149 : virtual CPLErr SetProjection( const char * );
150 : virtual int GetGCPCount();
151 : virtual const char *GetGCPProjection();
152 : virtual const GDAL_GCP *GetGCPs();
153 : };
154 :
155 : /************************************************************************/
156 : /* ==================================================================== */
157 : /* HDF4ImageRasterBand */
158 : /* ==================================================================== */
159 : /************************************************************************/
160 :
161 : class HDF4ImageRasterBand : public GDALPamRasterBand
162 493 : {
163 : friend class HDF4ImageDataset;
164 :
165 : int bNoDataSet;
166 : double dfNoDataValue;
167 :
168 : int bHaveScale, bHaveOffset;
169 : double dfScale;
170 : double dfOffset;
171 :
172 : CPLString osUnitType;
173 :
174 : public:
175 :
176 : HDF4ImageRasterBand( HDF4ImageDataset *, int, GDALDataType );
177 :
178 : virtual CPLErr IReadBlock( int, int, void * );
179 : virtual CPLErr IWriteBlock( int, int, void * );
180 : virtual GDALColorInterp GetColorInterpretation();
181 : virtual GDALColorTable *GetColorTable();
182 : virtual double GetNoDataValue( int * );
183 : virtual CPLErr SetNoDataValue( double );
184 : virtual double GetOffset( int *pbSuccess );
185 : virtual double GetScale( int *pbSuccess );
186 : virtual const char *GetUnitType();
187 : };
188 :
189 : /************************************************************************/
190 : /* HDF4ImageRasterBand() */
191 : /************************************************************************/
192 :
193 493 : HDF4ImageRasterBand::HDF4ImageRasterBand( HDF4ImageDataset *poDS, int nBand,
194 493 : GDALDataType eType )
195 :
196 : {
197 493 : this->poDS = poDS;
198 493 : this->nBand = nBand;
199 493 : eDataType = eType;
200 493 : bNoDataSet = FALSE;
201 493 : dfNoDataValue = -9999.0;
202 :
203 493 : bHaveScale = bHaveOffset = FALSE;
204 493 : dfScale = 1.0;
205 493 : dfOffset = 0.0;
206 :
207 493 : nBlockXSize = poDS->GetRasterXSize();
208 :
209 : // Aim for a block of about 1000000 pixels. Chunking up substantially
210 : // improves performance in some situations. For now we only chunk up for
211 : // SDS and EOS based datasets since other variations haven't been tested. #2208
212 986 : if( poDS->iDatasetType == HDF4_SDS ||
213 : poDS->iDatasetType == HDF4_EOS)
214 : {
215 : int nChunkSize =
216 493 : atoi( CPLGetConfigOption("HDF4_BLOCK_PIXELS", "1000000") );
217 :
218 493 : nBlockYSize = nChunkSize / poDS->GetRasterXSize();
219 493 : nBlockYSize = MAX(1,MIN(nBlockYSize,poDS->GetRasterYSize()));
220 : }
221 : else
222 : {
223 0 : nBlockYSize = 1;
224 : }
225 :
226 : /* HDF4_EOS:EOS_GRID case. We ensure that the block size matches */
227 : /* the raster width, as the IReadBlock() code can only handle multiple */
228 : /* blocks per row */
229 493 : if ( poDS->nBlockPreferredXSize == nBlockXSize &&
230 : poDS->nBlockPreferredYSize > 0 )
231 : {
232 6 : if (poDS->nBlockPreferredYSize == 1)
233 : {
234 : /* Avoid defaulting to tile reading when the preferred height is 1 */
235 : /* as it leads to very poor performance with : */
236 : /* ftp://e4ftl01u.ecs.nasa.gov/MODIS_Composites/MOLT/MOD13Q1.005/2006.06.10/MOD13Q1.A2006161.h21v13.005.2008234103220.hd */
237 2 : poDS->bReadTile = FALSE;
238 : }
239 : else
240 : {
241 4 : nBlockYSize = poDS->nBlockPreferredYSize;
242 : }
243 : }
244 :
245 : /* -------------------------------------------------------------------- */
246 : /* We need to avoid using the tile based api if we aren't */
247 : /* matching the tile size. (#4672) */
248 : /* -------------------------------------------------------------------- */
249 493 : if( nBlockXSize != poDS->nBlockPreferredXSize
250 : || nBlockYSize != poDS->nBlockPreferredYSize )
251 : {
252 489 : poDS->bReadTile = FALSE;
253 : }
254 493 : }
255 :
256 : /************************************************************************/
257 : /* IReadBlock() */
258 : /************************************************************************/
259 :
260 281 : CPLErr HDF4ImageRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
261 : void * pImage )
262 : {
263 281 : HDF4ImageDataset *poGDS = (HDF4ImageDataset *) poDS;
264 : int32 aiStart[H4_MAX_NC_DIMS], aiEdges[H4_MAX_NC_DIMS];
265 281 : CPLErr eErr = CE_None;
266 :
267 281 : CPLMutexHolderD(&hHDF4Mutex);
268 :
269 281 : if( poGDS->eAccess == GA_Update )
270 : {
271 : memset( pImage, 0,
272 0 : nBlockXSize * nBlockYSize * GDALGetDataTypeSize(eDataType) / 8 );
273 0 : return CE_None;
274 : }
275 :
276 : /* -------------------------------------------------------------------- */
277 : /* Work out some block oriented details. */
278 : /* -------------------------------------------------------------------- */
279 281 : int nYOff = nBlockYOff * nBlockYSize;
280 281 : int nYSize = MIN(nYOff + nBlockYSize, poDS->GetRasterYSize()) - nYOff;
281 281 : CPLAssert( nBlockXOff == 0 );
282 :
283 : /* -------------------------------------------------------------------- */
284 : /* HDF files with external data files, such as some landsat */
285 : /* products (eg. data/hdf/L1G) need to be told what directory */
286 : /* to look in to find the external files. Normally this is the */
287 : /* directory holding the hdf file. */
288 : /* -------------------------------------------------------------------- */
289 281 : HXsetdir(CPLGetPath(poGDS->pszFilename));
290 :
291 : /* -------------------------------------------------------------------- */
292 : /* Handle different configurations. */
293 : /* -------------------------------------------------------------------- */
294 281 : switch ( poGDS->iDatasetType )
295 : {
296 : case HDF4_SDS:
297 : {
298 : /* We avoid doing SDselect() / SDendaccess() for each block access */
299 : /* as this is very slow when zlib compression is used */
300 :
301 78 : if (poGDS->iSDS == FAIL)
302 44 : poGDS->iSDS = SDselect( poGDS->hSD, poGDS->iDataset );
303 :
304 : /* HDF rank:
305 : A rank 2 dataset is an image read in scan-line order (2D).
306 : A rank 3 dataset is a series of images which are read in
307 : an image at a time to form a volume.
308 : A rank 4 dataset may be thought of as a series of volumes.
309 :
310 : The "aiStart" array specifies the multi-dimensional index of the
311 : starting corner of the hyperslab to read. The values are zero
312 : based.
313 :
314 : The "edge" array specifies the number of values to read along
315 : each dimension of the hyperslab.
316 :
317 : The "iStride" array allows for sub-sampling along each
318 : dimension. If a iStride value is specified for a dimension,
319 : that many values will be skipped over when reading along that
320 : dimension. Specifying iStride = NULL in the C interface or
321 : iStride = 1 in either interface specifies contiguous reading
322 : of data. If the iStride values are set to 0, SDreaddata
323 : returns FAIL (or -1). No matter what iStride value is
324 : provided, data is always placed contiguously in buffer.
325 : */
326 78 : switch ( poGDS->iRank )
327 : {
328 : case 4: // 4Dim: volume-time
329 : // FIXME: needs sample file. Does not work currently.
330 0 : aiStart[3] = 0/* range: 0--aiDimSizes[3]-1 */;
331 0 : aiEdges[3] = 1;
332 0 : aiStart[2] = 0/* range: 0--aiDimSizes[2]-1 */;
333 0 : aiEdges[2] = 1;
334 0 : aiStart[1] = nYOff; aiEdges[1] = nYSize;
335 0 : aiStart[0] = nBlockXOff; aiEdges[0] = nBlockXSize;
336 0 : break;
337 : case 3: // 3Dim: volume
338 24 : aiStart[poGDS->iBandDim] = nBand - 1;
339 24 : aiEdges[poGDS->iBandDim] = 1;
340 :
341 24 : aiStart[poGDS->iYDim] = nYOff;
342 24 : aiEdges[poGDS->iYDim] = nYSize;
343 :
344 24 : aiStart[poGDS->iXDim] = nBlockXOff;
345 24 : aiEdges[poGDS->iXDim] = nBlockXSize;
346 24 : break;
347 : case 2: // 2Dim: rows/cols
348 54 : aiStart[poGDS->iYDim] = nYOff;
349 54 : aiEdges[poGDS->iYDim] = nYSize;
350 :
351 54 : aiStart[poGDS->iXDim] = nBlockXOff;
352 54 : aiEdges[poGDS->iXDim] = nBlockXSize;
353 54 : break;
354 : case 1: //1Dim:
355 0 : aiStart[poGDS->iXDim] = nBlockXOff;
356 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
357 : break;
358 : }
359 :
360 : // Read HDF SDS array
361 78 : if( SDreaddata( poGDS->iSDS, aiStart, NULL, aiEdges, pImage ) < 0 )
362 : {
363 : CPLError( CE_Failure, CPLE_AppDefined,
364 0 : "SDreaddata() failed for block." );
365 0 : eErr = CE_Failure;
366 : }
367 :
368 : //SDendaccess( iSDS );
369 : }
370 78 : break;
371 :
372 : case HDF4_GR:
373 : {
374 : int nDataTypeSize =
375 0 : GDALGetDataTypeSize(poGDS->GetDataType(poGDS->iNumType)) / 8;
376 : GByte *pbBuffer = (GByte *)
377 0 : CPLMalloc(nBlockXSize*nBlockYSize*poGDS->iRank*nBlockYSize);
378 : int i, j;
379 :
380 0 : aiStart[poGDS->iYDim] = nYOff;
381 0 : aiEdges[poGDS->iYDim] = nYSize;
382 :
383 0 : aiStart[poGDS->iXDim] = nBlockXOff;
384 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
385 :
386 0 : if( GRreadimage(poGDS->iGR, aiStart, NULL, aiEdges, pbBuffer) < 0 )
387 : {
388 : CPLError( CE_Failure, CPLE_AppDefined,
389 0 : "GRreaddata() failed for block." );
390 0 : eErr = CE_Failure;
391 : }
392 : else
393 : {
394 0 : for ( i = 0, j = (nBand - 1) * nDataTypeSize;
395 : i < nBlockXSize * nDataTypeSize;
396 : i += nDataTypeSize, j += poGDS->nBands * nDataTypeSize )
397 0 : memcpy( (GByte *)pImage + i, pbBuffer + j, nDataTypeSize );
398 : }
399 :
400 0 : CPLFree( pbBuffer );
401 : }
402 0 : break;
403 :
404 : case HDF4_EOS:
405 : {
406 203 : switch ( poGDS->iSubdatasetType )
407 : {
408 : case H4ST_EOS_GRID:
409 : {
410 : int32 hGD;
411 :
412 : hGD = GDattach( poGDS->hHDF4,
413 203 : poGDS->pszSubdatasetName );
414 203 : switch ( poGDS->iRank )
415 : {
416 : case 4: // 4Dim: volume
417 0 : aiStart[poGDS->i4Dim] =
418 : (nBand - 1)
419 0 : / poGDS->aiDimSizes[poGDS->iBandDim];
420 0 : aiEdges[poGDS->i4Dim] = 1;
421 :
422 0 : aiStart[poGDS->iBandDim] =
423 : (nBand - 1)
424 0 : % poGDS->aiDimSizes[poGDS->iBandDim];
425 0 : aiEdges[poGDS->iBandDim] = 1;
426 :
427 0 : aiStart[poGDS->iYDim] = nYOff;
428 0 : aiEdges[poGDS->iYDim] = nYSize;
429 :
430 0 : aiStart[poGDS->iXDim] = nBlockXOff;
431 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
432 0 : break;
433 : case 3: // 3Dim: volume
434 0 : aiStart[poGDS->iBandDim] = nBand - 1;
435 0 : aiEdges[poGDS->iBandDim] = 1;
436 :
437 0 : aiStart[poGDS->iYDim] = nYOff;
438 0 : aiEdges[poGDS->iYDim] = nYSize;
439 :
440 0 : aiStart[poGDS->iXDim] = nBlockXOff;
441 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
442 0 : break;
443 : case 2: // 2Dim: rows/cols
444 203 : aiStart[poGDS->iYDim] = nYOff;
445 203 : aiEdges[poGDS->iYDim] = nYSize;
446 :
447 203 : aiStart[poGDS->iXDim] = nBlockXOff;
448 203 : aiEdges[poGDS->iXDim] = nBlockXSize;
449 : break;
450 : }
451 :
452 : /* Ensure that we don't overlap the bottom or right edges */
453 : /* of the dataset in order to use the GDreadtile() API */
454 355 : if ( poGDS->bReadTile &&
455 : (nBlockXOff + 1) * nBlockXSize <= nRasterXSize &&
456 : (nBlockYOff + 1) * nBlockYSize <= nRasterYSize )
457 : {
458 152 : int32 tilecoords[] = { nBlockYOff , nBlockXOff };
459 152 : if( GDreadtile( hGD, poGDS->pszFieldName , tilecoords , pImage ) != 0 )
460 : {
461 0 : CPLError( CE_Failure, CPLE_AppDefined, "GDreadtile() failed for block." );
462 0 : eErr = CE_Failure;
463 : }
464 : }
465 51 : else if( GDreadfield( hGD, poGDS->pszFieldName,
466 : aiStart, NULL, aiEdges, pImage ) < 0 )
467 : {
468 : CPLError( CE_Failure, CPLE_AppDefined,
469 0 : "GDreadfield() failed for block." );
470 0 : eErr = CE_Failure;
471 : }
472 203 : GDdetach( hGD );
473 : }
474 203 : break;
475 :
476 : case H4ST_EOS_SWATH:
477 : case H4ST_EOS_SWATH_GEOL:
478 : {
479 : int32 hSW;
480 :
481 : hSW = SWattach( poGDS->hHDF4,
482 0 : poGDS->pszSubdatasetName );
483 0 : switch ( poGDS->iRank )
484 : {
485 : case 3: // 3Dim: volume
486 0 : aiStart[poGDS->iBandDim] = nBand - 1;
487 0 : aiEdges[poGDS->iBandDim] = 1;
488 :
489 0 : aiStart[poGDS->iYDim] = nYOff;
490 0 : aiEdges[poGDS->iYDim] = nYSize;
491 :
492 0 : aiStart[poGDS->iXDim] = nBlockXOff;
493 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
494 0 : break;
495 : case 2: // 2Dim: rows/cols
496 0 : aiStart[poGDS->iYDim] = nYOff;
497 0 : aiEdges[poGDS->iYDim] = nYSize;
498 :
499 0 : aiStart[poGDS->iXDim] = nBlockXOff;
500 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
501 : break;
502 : }
503 0 : if( SWreadfield( hSW, poGDS->pszFieldName,
504 : aiStart, NULL, aiEdges, pImage ) < 0 )
505 : {
506 : CPLError( CE_Failure, CPLE_AppDefined,
507 0 : "SWreadfield() failed for block." );
508 0 : eErr = CE_Failure;
509 : }
510 0 : SWdetach( hSW );
511 : }
512 : break;
513 : default:
514 : break;
515 : }
516 : }
517 203 : break;
518 :
519 : default:
520 0 : eErr = CE_Failure;
521 : break;
522 : }
523 :
524 281 : return eErr;
525 : }
526 :
527 : /************************************************************************/
528 : /* IWriteBlock() */
529 : /************************************************************************/
530 :
531 49 : CPLErr HDF4ImageRasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff,
532 : void * pImage )
533 : {
534 49 : HDF4ImageDataset *poGDS = (HDF4ImageDataset *)poDS;
535 : int32 aiStart[H4_MAX_NC_DIMS], aiEdges[H4_MAX_NC_DIMS];
536 49 : CPLErr eErr = CE_None;
537 :
538 49 : CPLMutexHolderD(&hHDF4Mutex);
539 :
540 : CPLAssert( poGDS != NULL
541 : && nBlockXOff == 0
542 : && nBlockYOff >= 0
543 49 : && pImage != NULL );
544 :
545 : /* -------------------------------------------------------------------- */
546 : /* Work out some block oriented details. */
547 : /* -------------------------------------------------------------------- */
548 49 : int nYOff = nBlockYOff * nBlockYSize;
549 49 : int nYSize = MIN(nYOff + nBlockYSize, poDS->GetRasterYSize()) - nYOff;
550 49 : CPLAssert( nBlockXOff == 0 );
551 :
552 : /* -------------------------------------------------------------------- */
553 : /* Process based on rank. */
554 : /* -------------------------------------------------------------------- */
555 49 : switch ( poGDS->iRank )
556 : {
557 : case 3:
558 : {
559 41 : int32 iSDS = SDselect( poGDS->hSD, poGDS->iDataset );
560 :
561 41 : aiStart[poGDS->iBandDim] = nBand - 1;
562 41 : aiEdges[poGDS->iBandDim] = 1;
563 :
564 41 : aiStart[poGDS->iYDim] = nYOff;
565 41 : aiEdges[poGDS->iYDim] = nYSize;
566 :
567 41 : aiStart[poGDS->iXDim] = nBlockXOff;
568 41 : aiEdges[poGDS->iXDim] = nBlockXSize;
569 :
570 41 : if ( (SDwritedata( iSDS, aiStart, NULL,
571 : aiEdges, (VOIDP)pImage )) < 0 )
572 0 : eErr = CE_Failure;
573 :
574 41 : SDendaccess( iSDS );
575 : }
576 41 : break;
577 :
578 : case 2:
579 : {
580 8 : int32 iSDS = SDselect( poGDS->hSD, nBand - 1 );
581 8 : aiStart[poGDS->iYDim] = nYOff;
582 8 : aiEdges[poGDS->iYDim] = nYSize;
583 :
584 8 : aiStart[poGDS->iXDim] = nBlockXOff;
585 8 : aiEdges[poGDS->iXDim] = nBlockXSize;
586 :
587 8 : if ( (SDwritedata( iSDS, aiStart, NULL,
588 : aiEdges, (VOIDP)pImage )) < 0 )
589 0 : eErr = CE_Failure;
590 :
591 8 : SDendaccess( iSDS );
592 : }
593 8 : break;
594 :
595 : default:
596 0 : eErr = CE_Failure;
597 : break;
598 : }
599 :
600 49 : return eErr;
601 : }
602 :
603 : /************************************************************************/
604 : /* GetColorTable() */
605 : /************************************************************************/
606 :
607 0 : GDALColorTable *HDF4ImageRasterBand::GetColorTable()
608 : {
609 0 : HDF4ImageDataset *poGDS = (HDF4ImageDataset *) poDS;
610 :
611 0 : return poGDS->poColorTable;
612 : }
613 :
614 : /************************************************************************/
615 : /* GetColorInterpretation() */
616 : /************************************************************************/
617 :
618 36 : GDALColorInterp HDF4ImageRasterBand::GetColorInterpretation()
619 : {
620 36 : HDF4ImageDataset *poGDS = (HDF4ImageDataset *) poDS;
621 :
622 36 : if ( poGDS->iDatasetType == HDF4_SDS )
623 36 : return GCI_GrayIndex;
624 0 : else if ( poGDS->iDatasetType == HDF4_GR )
625 : {
626 0 : if ( poGDS->poColorTable != NULL )
627 0 : return GCI_PaletteIndex;
628 0 : else if ( poGDS->nBands != 1 )
629 : {
630 0 : if ( nBand == 1 )
631 0 : return GCI_RedBand;
632 0 : else if ( nBand == 2 )
633 0 : return GCI_GreenBand;
634 0 : else if ( nBand == 3 )
635 0 : return GCI_BlueBand;
636 0 : else if ( nBand == 4 )
637 0 : return GCI_AlphaBand;
638 : else
639 0 : return GCI_Undefined;
640 : }
641 : else
642 0 : return GCI_GrayIndex;
643 : }
644 : else
645 0 : return GCI_GrayIndex;
646 : }
647 :
648 : /************************************************************************/
649 : /* GetNoDataValue() */
650 : /************************************************************************/
651 :
652 64 : double HDF4ImageRasterBand::GetNoDataValue( int * pbSuccess )
653 :
654 : {
655 64 : if( pbSuccess )
656 64 : *pbSuccess = bNoDataSet;
657 :
658 64 : return dfNoDataValue;
659 : }
660 :
661 : /************************************************************************/
662 : /* SetNoDataValue() */
663 : /************************************************************************/
664 :
665 55 : CPLErr HDF4ImageRasterBand::SetNoDataValue( double dfNoData )
666 :
667 : {
668 55 : bNoDataSet = TRUE;
669 55 : dfNoDataValue = dfNoData;
670 :
671 55 : return CE_None;
672 : }
673 :
674 : /************************************************************************/
675 : /* GetUnitType() */
676 : /************************************************************************/
677 :
678 0 : const char *HDF4ImageRasterBand::GetUnitType()
679 :
680 : {
681 0 : if( osUnitType.size() > 0 )
682 0 : return osUnitType;
683 : else
684 0 : return GDALRasterBand::GetUnitType();
685 : }
686 :
687 : /************************************************************************/
688 : /* GetOffset() */
689 : /************************************************************************/
690 :
691 0 : double HDF4ImageRasterBand::GetOffset( int *pbSuccess )
692 :
693 : {
694 0 : if( bHaveOffset )
695 : {
696 0 : if( pbSuccess != NULL )
697 0 : *pbSuccess = TRUE;
698 0 : return dfOffset;
699 : }
700 : else
701 0 : return GDALRasterBand::GetOffset( pbSuccess );
702 : }
703 :
704 : /************************************************************************/
705 : /* GetScale() */
706 : /************************************************************************/
707 :
708 0 : double HDF4ImageRasterBand::GetScale( int *pbSuccess )
709 :
710 : {
711 0 : if( bHaveScale )
712 : {
713 0 : if( pbSuccess != NULL )
714 0 : *pbSuccess = TRUE;
715 0 : return dfScale;
716 : }
717 : else
718 0 : return GDALRasterBand::GetScale( pbSuccess );
719 : }
720 :
721 : /************************************************************************/
722 : /* ==================================================================== */
723 : /* HDF4ImageDataset */
724 : /* ==================================================================== */
725 : /************************************************************************/
726 :
727 : /************************************************************************/
728 : /* HDF4ImageDataset() */
729 : /************************************************************************/
730 :
731 426 : HDF4ImageDataset::HDF4ImageDataset()
732 : {
733 426 : pszFilename = NULL;
734 426 : hHDF4 = 0;
735 426 : iGR = 0;
736 426 : iPal = 0;
737 426 : iDataset = 0;
738 426 : iRank = 0;
739 426 : iNumType = 0;
740 426 : nAttrs = 0;
741 426 : iInterlaceMode = 0;
742 426 : iPalInterlaceMode = 0;
743 426 : iPalDataType = 0;
744 426 : nComps = 0;
745 426 : nPalEntries = 0;
746 426 : memset(aiDimSizes, 0, sizeof(aiDimSizes));
747 426 : iXDim = 0;
748 426 : iYDim = 0;
749 426 : iBandDim = -1;
750 426 : i4Dim = 0;
751 426 : nBandCount = 0;
752 426 : papszLocalMetadata = NULL;
753 426 : memset(aiPaletteData, 0, sizeof(aiPaletteData));
754 426 : memset(szName, 0, sizeof(szName));
755 426 : pszSubdatasetName = NULL;
756 426 : pszFieldName = NULL;
757 426 : poColorTable = NULL;
758 426 : bHasGeoTransform = FALSE;
759 426 : adfGeoTransform[0] = 0.0;
760 426 : adfGeoTransform[1] = 1.0;
761 426 : adfGeoTransform[2] = 0.0;
762 426 : adfGeoTransform[3] = 0.0;
763 426 : adfGeoTransform[4] = 0.0;
764 426 : adfGeoTransform[5] = 1.0;
765 426 : pszProjection = CPLStrdup( "" );
766 426 : pszGCPProjection = CPLStrdup( "" );
767 426 : pasGCPList = NULL;
768 426 : nGCPCount = 0;
769 :
770 426 : iDatasetType = HDF4_UNKNOWN;
771 426 : iSDS = FAIL;
772 :
773 426 : nBlockPreferredXSize = -1;
774 426 : nBlockPreferredYSize = -1;
775 426 : bReadTile = false;
776 :
777 426 : }
778 :
779 : /************************************************************************/
780 : /* ~HDF4ImageDataset() */
781 : /************************************************************************/
782 :
783 426 : HDF4ImageDataset::~HDF4ImageDataset()
784 : {
785 426 : CPLMutexHolderD(&hHDF4Mutex);
786 :
787 426 : FlushCache();
788 :
789 426 : if ( pszFilename )
790 265 : CPLFree( pszFilename );
791 426 : if ( iSDS != FAIL )
792 44 : SDendaccess( iSDS );
793 426 : if ( hSD > 0 )
794 409 : SDend( hSD );
795 426 : hSD = 0;
796 426 : if ( iGR > 0 )
797 0 : GRendaccess( iGR );
798 426 : if ( hGR > 0 )
799 0 : GRend( hGR );
800 426 : hGR = 0;
801 426 : if ( pszSubdatasetName )
802 9 : CPLFree( pszSubdatasetName );
803 426 : if ( pszFieldName )
804 9 : CPLFree( pszFieldName );
805 426 : if ( papszLocalMetadata )
806 264 : CSLDestroy( papszLocalMetadata );
807 426 : if ( poColorTable != NULL )
808 0 : delete poColorTable;
809 426 : if ( pszProjection )
810 426 : CPLFree( pszProjection );
811 426 : if ( pszGCPProjection )
812 426 : CPLFree( pszGCPProjection );
813 426 : if( nGCPCount > 0 )
814 : {
815 5 : for( int i = 0; i < nGCPCount; i++ )
816 : {
817 4 : if ( pasGCPList[i].pszId )
818 4 : CPLFree( pasGCPList[i].pszId );
819 4 : if ( pasGCPList[i].pszInfo )
820 4 : CPLFree( pasGCPList[i].pszInfo );
821 : }
822 :
823 1 : CPLFree( pasGCPList );
824 : }
825 426 : if ( hHDF4 > 0 )
826 : {
827 265 : switch ( iDatasetType )
828 : {
829 : case HDF4_EOS:
830 9 : switch ( iSubdatasetType )
831 : {
832 : case H4ST_EOS_SWATH:
833 : case H4ST_EOS_SWATH_GEOL:
834 0 : SWclose( hHDF4 );
835 0 : break;
836 : case H4ST_EOS_GRID:
837 9 : GDclose( hHDF4 );
838 : default:
839 : break;
840 : }
841 9 : break;
842 : case HDF4_SDS:
843 : case HDF4_GR:
844 256 : hHDF4 = Hclose( hHDF4 );
845 : break;
846 : default:
847 : break;
848 : }
849 426 : }
850 426 : }
851 :
852 : /************************************************************************/
853 : /* GetGeoTransform() */
854 : /************************************************************************/
855 :
856 42 : CPLErr HDF4ImageDataset::GetGeoTransform( double * padfTransform )
857 : {
858 42 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
859 :
860 42 : if ( !bHasGeoTransform )
861 0 : return CE_Failure;
862 :
863 42 : return CE_None;
864 : }
865 :
866 : /************************************************************************/
867 : /* SetGeoTransform() */
868 : /************************************************************************/
869 :
870 88 : CPLErr HDF4ImageDataset::SetGeoTransform( double * padfTransform )
871 : {
872 88 : bHasGeoTransform = TRUE;
873 88 : memcpy( adfGeoTransform, padfTransform, sizeof(double) * 6 );
874 :
875 88 : return CE_None;
876 : }
877 :
878 : /************************************************************************/
879 : /* GetProjectionRef() */
880 : /************************************************************************/
881 :
882 17 : const char *HDF4ImageDataset::GetProjectionRef()
883 :
884 : {
885 17 : return pszProjection;
886 : }
887 :
888 : /************************************************************************/
889 : /* SetProjection() */
890 : /************************************************************************/
891 :
892 72 : CPLErr HDF4ImageDataset::SetProjection( const char *pszNewProjection )
893 :
894 : {
895 72 : if ( pszProjection )
896 72 : CPLFree( pszProjection );
897 72 : pszProjection = CPLStrdup( pszNewProjection );
898 :
899 72 : return CE_None;
900 : }
901 :
902 : /************************************************************************/
903 : /* GetGCPCount() */
904 : /************************************************************************/
905 :
906 1 : int HDF4ImageDataset::GetGCPCount()
907 :
908 : {
909 1 : return nGCPCount;
910 : }
911 :
912 : /************************************************************************/
913 : /* GetGCPProjection() */
914 : /************************************************************************/
915 :
916 0 : const char *HDF4ImageDataset::GetGCPProjection()
917 :
918 : {
919 0 : if( nGCPCount > 0 )
920 0 : return pszGCPProjection;
921 : else
922 0 : return "";
923 : }
924 :
925 : /************************************************************************/
926 : /* GetGCPs() */
927 : /************************************************************************/
928 :
929 0 : const GDAL_GCP *HDF4ImageDataset::GetGCPs()
930 : {
931 0 : return pasGCPList;
932 : }
933 :
934 : /************************************************************************/
935 : /* FlushCache() */
936 : /************************************************************************/
937 :
938 426 : void HDF4ImageDataset::FlushCache()
939 :
940 : {
941 : int iBand;
942 : char *pszName;
943 : const char *pszValue;
944 :
945 426 : CPLMutexHolderD(&hHDF4Mutex);
946 :
947 426 : GDALDataset::FlushCache();
948 :
949 426 : if( eAccess == GA_ReadOnly )
950 : return;
951 :
952 : // Write out transformation matrix
953 : pszValue = CPLSPrintf( "%f, %f, %f, %f, %f, %f",
954 : adfGeoTransform[0], adfGeoTransform[1],
955 : adfGeoTransform[2], adfGeoTransform[3],
956 144 : adfGeoTransform[4], adfGeoTransform[5] );
957 144 : if ( (SDsetattr( hSD, "TransformationMatrix", DFNT_CHAR8,
958 : strlen(pszValue) + 1, pszValue )) < 0 )
959 : {
960 : CPLDebug( "HDF4Image",
961 0 : "Cannot write transformation matrix to output file" );
962 : }
963 :
964 : // Write out projection
965 144 : if ( pszProjection != NULL && !EQUAL( pszProjection, "" ) )
966 : {
967 72 : if ( (SDsetattr( hSD, "Projection", DFNT_CHAR8,
968 : strlen(pszProjection) + 1, pszProjection )) < 0 )
969 : {
970 : CPLDebug( "HDF4Image",
971 0 : "Cannot write projection information to output file");
972 : }
973 : }
974 :
975 : // Store all metadata from source dataset as HDF attributes
976 144 : if( GetMetadata() )
977 : {
978 47 : char **papszMeta = GetMetadata();
979 :
980 141 : while ( *papszMeta )
981 : {
982 47 : pszName = NULL;
983 47 : pszValue = CPLParseNameValue( *papszMeta++, &pszName );
984 47 : if ( pszName != NULL && (SDsetattr( hSD, pszName, DFNT_CHAR8,
985 : strlen(pszValue) + 1, pszValue )) < 0 )
986 : {
987 : CPLDebug( "HDF4Image",
988 0 : "Cannot write metadata information to output file");
989 : }
990 :
991 47 : CPLFree( pszName );
992 : }
993 : }
994 :
995 : // Write out NoData values
996 328 : for ( iBand = 1; iBand <= nBands; iBand++ )
997 : {
998 : HDF4ImageRasterBand *poBand =
999 184 : (HDF4ImageRasterBand *)GetRasterBand(iBand);
1000 :
1001 184 : if ( poBand->bNoDataSet )
1002 : {
1003 16 : pszName = CPLStrdup( CPLSPrintf( "NoDataValue%d", iBand ) );
1004 16 : pszValue = CPLSPrintf( "%f", poBand->dfNoDataValue );
1005 16 : if ( (SDsetattr( hSD, pszName, DFNT_CHAR8,
1006 : strlen(pszValue) + 1, pszValue )) < 0 )
1007 : {
1008 : CPLDebug( "HDF4Image",
1009 : "Cannot write NoData value for band %d "
1010 0 : "to output file", iBand);
1011 : }
1012 :
1013 16 : CPLFree( pszName );
1014 : }
1015 : }
1016 :
1017 : // Write out band descriptions
1018 328 : for ( iBand = 1; iBand <= nBands; iBand++ )
1019 : {
1020 : HDF4ImageRasterBand *poBand =
1021 184 : (HDF4ImageRasterBand *)GetRasterBand(iBand);
1022 :
1023 184 : pszName = CPLStrdup( CPLSPrintf( "BandDesc%d", iBand ) );
1024 184 : pszValue = poBand->GetDescription();
1025 184 : if ( pszValue != NULL && !EQUAL( pszValue, "" ) )
1026 : {
1027 16 : if ( (SDsetattr( hSD, pszName, DFNT_CHAR8,
1028 : strlen(pszValue) + 1, pszValue )) < 0 )
1029 : {
1030 : CPLDebug( "HDF4Image",
1031 : "Cannot write band's %d description to output file",
1032 0 : iBand);
1033 : }
1034 : }
1035 :
1036 184 : CPLFree( pszName );
1037 0 : }
1038 : }
1039 :
1040 : /************************************************************************/
1041 : /* USGSMnemonicToCode() */
1042 : /************************************************************************/
1043 :
1044 0 : long HDF4ImageDataset::USGSMnemonicToCode( const char* pszMnemonic )
1045 : {
1046 0 : if ( EQUAL(pszMnemonic, "UTM") )
1047 0 : return 1L;
1048 0 : else if ( EQUAL(pszMnemonic, "LAMCC") )
1049 0 : return 4L;
1050 0 : else if ( EQUAL(pszMnemonic, "PS") )
1051 0 : return 6L;
1052 0 : else if ( EQUAL(pszMnemonic, "PC") )
1053 0 : return 7L;
1054 0 : else if ( EQUAL(pszMnemonic, "TM") )
1055 0 : return 9L;
1056 0 : else if ( EQUAL(pszMnemonic, "EQRECT") )
1057 0 : return 17L;
1058 0 : else if ( EQUAL(pszMnemonic, "OM") )
1059 0 : return 20L;
1060 0 : else if ( EQUAL(pszMnemonic, "SOM") )
1061 0 : return 22L;
1062 : else
1063 0 : return 1L; // UTM by default
1064 : }
1065 :
1066 : /************************************************************************/
1067 : /* ToGeoref() */
1068 : /************************************************************************/
1069 :
1070 0 : void HDF4ImageDataset::ToGeoref( double *pdfGeoX, double *pdfGeoY )
1071 : {
1072 0 : OGRCoordinateTransformation *poTransform = NULL;
1073 0 : OGRSpatialReference *poLatLong = NULL;
1074 0 : poLatLong = oSRS.CloneGeogCS();
1075 0 : poTransform = OGRCreateCoordinateTransformation( poLatLong, &oSRS );
1076 :
1077 0 : if( poTransform != NULL )
1078 0 : poTransform->Transform( 1, pdfGeoX, pdfGeoY, NULL );
1079 :
1080 0 : if( poTransform != NULL )
1081 0 : delete poTransform;
1082 :
1083 0 : if( poLatLong != NULL )
1084 0 : delete poLatLong;
1085 0 : }
1086 :
1087 : /************************************************************************/
1088 : /* ReadCoordinates() */
1089 : /************************************************************************/
1090 :
1091 0 : void HDF4ImageDataset::ReadCoordinates( const char *pszString,
1092 : double *pdfX, double *pdfY )
1093 : {
1094 : char **papszStrList;
1095 0 : papszStrList = CSLTokenizeString2( pszString, ", ", 0 );
1096 0 : *pdfX = CPLAtof( papszStrList[0] );
1097 0 : *pdfY = CPLAtof( papszStrList[1] );
1098 0 : CSLDestroy( papszStrList );
1099 0 : }
1100 :
1101 : /************************************************************************/
1102 : /* CaptureL1GMTLInfo() */
1103 : /************************************************************************/
1104 :
1105 : /* FILE L71002025_02520010722_M_MTL.L1G
1106 :
1107 : GROUP = L1_METADATA_FILE
1108 : ...
1109 : GROUP = PRODUCT_METADATA
1110 : PRODUCT_TYPE = "L1G"
1111 : PROCESSING_SOFTWARE = "IAS_5.1"
1112 : EPHEMERIS_TYPE = "DEFINITIVE"
1113 : SPACECRAFT_ID = "Landsat7"
1114 : SENSOR_ID = "ETM+"
1115 : ACQUISITION_DATE = 2001-07-22
1116 : WRS_PATH = 002
1117 : STARTING_ROW = 025
1118 : ENDING_ROW = 025
1119 : BAND_COMBINATION = "12345--7-"
1120 : PRODUCT_UL_CORNER_LAT = 51.2704805
1121 : PRODUCT_UL_CORNER_LON = -53.8914311
1122 : PRODUCT_UR_CORNER_LAT = 50.8458100
1123 : PRODUCT_UR_CORNER_LON = -50.9869091
1124 : PRODUCT_LL_CORNER_LAT = 49.6960897
1125 : PRODUCT_LL_CORNER_LON = -54.4047933
1126 : PRODUCT_LR_CORNER_LAT = 49.2841436
1127 : PRODUCT_LR_CORNER_LON = -51.5900428
1128 : PRODUCT_UL_CORNER_MAPX = 298309.894
1129 : PRODUCT_UL_CORNER_MAPY = 5683875.631
1130 : PRODUCT_UR_CORNER_MAPX = 500921.624
1131 : PRODUCT_UR_CORNER_MAPY = 5632678.683
1132 : PRODUCT_LL_CORNER_MAPX = 254477.193
1133 : PRODUCT_LL_CORNER_MAPY = 5510407.880
1134 : PRODUCT_LR_CORNER_MAPX = 457088.923
1135 : PRODUCT_LR_CORNER_MAPY = 5459210.932
1136 : PRODUCT_SAMPLES_REF = 6967
1137 : PRODUCT_LINES_REF = 5965
1138 : BAND1_FILE_NAME = "L71002025_02520010722_B10.L1G"
1139 : BAND2_FILE_NAME = "L71002025_02520010722_B20.L1G"
1140 : BAND3_FILE_NAME = "L71002025_02520010722_B30.L1G"
1141 : BAND4_FILE_NAME = "L71002025_02520010722_B40.L1G"
1142 : BAND5_FILE_NAME = "L71002025_02520010722_B50.L1G"
1143 : BAND7_FILE_NAME = "L72002025_02520010722_B70.L1G"
1144 : METADATA_L1_FILE_NAME = "L71002025_02520010722_MTL.L1G"
1145 : CPF_FILE_NAME = "L7CPF20010701_20010930_06"
1146 : HDF_DIR_FILE_NAME = "L71002025_02520010722_HDF.L1G"
1147 : END_GROUP = PRODUCT_METADATA
1148 : ...
1149 : GROUP = PROJECTION_PARAMETERS
1150 : REFERENCE_DATUM = "NAD83"
1151 : REFERENCE_ELLIPSOID = "GRS80"
1152 : GRID_CELL_SIZE_PAN = 15.000000
1153 : GRID_CELL_SIZE_THM = 60.000000
1154 : GRID_CELL_SIZE_REF = 30.000000
1155 : ORIENTATION = "NOM"
1156 : RESAMPLING_OPTION = "CC"
1157 : MAP_PROJECTION = "UTM"
1158 : END_GROUP = PROJECTION_PARAMETERS
1159 : GROUP = UTM_PARAMETERS
1160 : ZONE_NUMBER = 22
1161 : END_GROUP = UTM_PARAMETERS
1162 : END_GROUP = L1_METADATA_FILE
1163 : END
1164 : */
1165 :
1166 256 : void HDF4ImageDataset::CaptureL1GMTLInfo()
1167 :
1168 : {
1169 : /* -------------------------------------------------------------------- */
1170 : /* Does the physical file look like it matches our expected */
1171 : /* name pattern? */
1172 : /* -------------------------------------------------------------------- */
1173 256 : if( strlen(pszFilename) < 8
1174 : || !EQUAL(pszFilename+strlen(pszFilename)-8,"_HDF.L1G") )
1175 255 : return;
1176 :
1177 : /* -------------------------------------------------------------------- */
1178 : /* Construct the name of the corresponding MTL file. We should */
1179 : /* likely be able to extract that from the HDF itself but I'm */
1180 : /* not sure where to find it. */
1181 : /* -------------------------------------------------------------------- */
1182 1 : CPLString osMTLFilename = pszFilename;
1183 1 : osMTLFilename.resize(osMTLFilename.length() - 8);
1184 1 : osMTLFilename += "_MTL.L1G";
1185 :
1186 : /* -------------------------------------------------------------------- */
1187 : /* Ingest the MTL using the NASAKeywordHandler written for the */
1188 : /* PDS driver. */
1189 : /* -------------------------------------------------------------------- */
1190 1 : NASAKeywordHandler oMTL;
1191 :
1192 1 : VSILFILE *fp = VSIFOpenL( osMTLFilename, "r" );
1193 1 : if( fp == NULL )
1194 : return;
1195 :
1196 1 : if( !oMTL.Ingest( fp, 0 ) )
1197 : {
1198 0 : VSIFCloseL( fp );
1199 : return;
1200 : }
1201 :
1202 1 : VSIFCloseL( fp );
1203 :
1204 : /* -------------------------------------------------------------------- */
1205 : /* Note: Different variation of MTL files use different group names. */
1206 : /* Check for LPGS_METADATA_FILE and L1_METADATA_FILE. */
1207 : /* -------------------------------------------------------------------- */
1208 : double dfULX, dfULY, dfLRX, dfLRY;
1209 : double dfLLX, dfLLY, dfURX, dfURY;
1210 1 : CPLString osPrefix;
1211 :
1212 1 : if( oMTL.GetKeyword( "LPGS_METADATA_FILE.PRODUCT_METADATA"
1213 : ".PRODUCT_UL_CORNER_LON", NULL ) )
1214 0 : osPrefix = "LPGS_METADATA_FILE.PRODUCT_METADATA.PRODUCT_";
1215 1 : else if( oMTL.GetKeyword( "L1_METADATA_FILE.PRODUCT_METADATA"
1216 : ".PRODUCT_UL_CORNER_LON", NULL ) )
1217 1 : osPrefix = "L1_METADATA_FILE.PRODUCT_METADATA.PRODUCT_";
1218 : else
1219 : return;
1220 :
1221 1 : dfULX = atof(oMTL.GetKeyword((osPrefix+"UL_CORNER_LON").c_str(), "0" ));
1222 1 : dfULY = atof(oMTL.GetKeyword((osPrefix+"UL_CORNER_LAT").c_str(), "0" ));
1223 1 : dfLRX = atof(oMTL.GetKeyword((osPrefix+"LR_CORNER_LON").c_str(), "0" ));
1224 1 : dfLRY = atof(oMTL.GetKeyword((osPrefix+"LR_CORNER_LAT").c_str(), "0" ));
1225 1 : dfLLX = atof(oMTL.GetKeyword((osPrefix+"LL_CORNER_LON").c_str(), "0" ));
1226 1 : dfLLY = atof(oMTL.GetKeyword((osPrefix+"LL_CORNER_LAT").c_str(), "0" ));
1227 1 : dfURX = atof(oMTL.GetKeyword((osPrefix+"UR_CORNER_LON").c_str(), "0" ));
1228 1 : dfURY = atof(oMTL.GetKeyword((osPrefix+"UR_CORNER_LAT").c_str(), "0" ));
1229 :
1230 1 : CPLFree( pszGCPProjection );
1231 1 : pszGCPProjection = CPLStrdup( "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9108\"]],AXIS[\"Lat\",NORTH],AXIS[\"Long\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]" );
1232 :
1233 1 : nGCPCount = 4;
1234 1 : pasGCPList = (GDAL_GCP *) CPLCalloc( nGCPCount, sizeof( GDAL_GCP ) );
1235 1 : GDALInitGCPs( nGCPCount, pasGCPList );
1236 :
1237 1 : pasGCPList[0].dfGCPX = dfULX;
1238 1 : pasGCPList[0].dfGCPY = dfULY;
1239 1 : pasGCPList[0].dfGCPPixel = 0.0;
1240 1 : pasGCPList[0].dfGCPLine = 0.0;
1241 :
1242 1 : pasGCPList[1].dfGCPX = dfURX;
1243 1 : pasGCPList[1].dfGCPY = dfURY;
1244 1 : pasGCPList[1].dfGCPPixel = GetRasterXSize();
1245 1 : pasGCPList[1].dfGCPLine = 0.0;
1246 :
1247 1 : pasGCPList[2].dfGCPX = dfLLX;
1248 1 : pasGCPList[2].dfGCPY = dfLLY;
1249 1 : pasGCPList[2].dfGCPPixel = 0.0;
1250 1 : pasGCPList[2].dfGCPLine = GetRasterYSize();
1251 :
1252 1 : pasGCPList[3].dfGCPX = dfLRX;
1253 1 : pasGCPList[3].dfGCPY = dfLRY;
1254 1 : pasGCPList[3].dfGCPPixel = GetRasterXSize();
1255 1 : pasGCPList[3].dfGCPLine = GetRasterYSize();
1256 : }
1257 :
1258 : /************************************************************************/
1259 : /* CaptureNRLGeoTransform() */
1260 : /* */
1261 : /* Capture geotransform and coordinate system from NRL (Naval */
1262 : /* Research Laboratory, Stennis Space Center) metadata. */
1263 : /************************************************************************/
1264 :
1265 : /* Example metadata:
1266 : Metadata:
1267 : createTime=Fri Oct 1 18:00:07 2004
1268 : createSoftware=APS v2.8.4
1269 : createPlatform=i686-pc-linux-gnu
1270 : createAgency=Naval Research Laboratory, Stennis Space Center
1271 : sensor=MODIS
1272 : sensorPlatform=TERRA-AM
1273 : sensorAgency=NASA
1274 : sensorType=whiskbroom
1275 : sensorSpectrum=Visible/Thermal
1276 : sensorNumberOfBands=36
1277 : sensorBandUnits=nano meters
1278 : sensorBands=645, 858.5, 469, 555, 1240, 1640, 2130, 412.5, 443, 488, 531, 551,
1279 : 667, 678, 748, 869.5, 905, 936, 940, 3750, 3959, 3959, 4050, 4465.5, 4515.5, 13
1280 : 75, 6715, 7325, 8550, 9730, 11130, 12020, 13335, 13635, 13935, 14235
1281 : sensorBandWidths=50, 35, 20, 20, 20, 24, 50, 15, 10, 10, 10, 10, 10, 10, 10, 1
1282 : 5, 30, 10, 50, 90, 60, 60, 60, 65, 67, 30, 360, 300, 300, 300, 500, 500, 300, 30
1283 : 0, 300, 300
1284 : sensorNominalAltitudeInKM=705
1285 : sensorScanWidthInKM=2330
1286 : sensorResolutionInKM=1
1287 : sensorPlatformType=Polar-orbiting Satellite
1288 : timeStartYear=2004
1289 : timeStartDay=275
1290 : timeStartTime=56400000
1291 : timeStart=Fri Oct 1 15:40:00 2004
1292 : timeDayNight=Day
1293 : timeEndYear=2004
1294 : timeEndDay=275
1295 : timeEndTime=56700000
1296 : timeEnd=Fri Oct 1 15:45:00 2004
1297 : inputMasks=HIGLINT,CLDICE,LAND,ATMFAIL
1298 : inputMasksInt=523
1299 : processedVersion=1.2
1300 : file=MODAM2004275.L3_Mosaic_NOAA_GMX
1301 : fileTitle=NRL Level-3 Mosaic
1302 : fileVersion=3.0
1303 : fileClassification=UNCLASSIFIED
1304 : fileStatus=EXPERIMENTAL
1305 : navType=mapped
1306 : mapProjectionSystem=NRL(USGS)
1307 : mapProjection=Gomex
1308 : mapUpperLeft=31, -99
1309 : mapUpperRight=31, -79
1310 : mapLowerLeft=14.9844128048645, -99
1311 : mapLowerRight=14.9844128048645, -79
1312 : inputFiles=MODAM2004275154000.L3_NOAA_GMX
1313 : ...
1314 : */
1315 :
1316 0 : void HDF4ImageDataset::CaptureNRLGeoTransform()
1317 :
1318 : {
1319 : /* -------------------------------------------------------------------- */
1320 : /* Collect the four corners. */
1321 : /* -------------------------------------------------------------------- */
1322 : double adfXY[8];
1323 : static const char *apszItems[] = {
1324 : "mapUpperLeft", "mapUpperRight", "mapLowerLeft", "mapLowerRight" };
1325 : int iCorner;
1326 0 : int bLLPossible = TRUE;
1327 :
1328 0 : for( iCorner = 0; iCorner < 4; iCorner++ )
1329 : {
1330 : const char *pszCornerLoc =
1331 0 : CSLFetchNameValue( papszGlobalMetadata, apszItems[iCorner] );
1332 :
1333 0 : if( pszCornerLoc == NULL )
1334 0 : return;
1335 :
1336 : char **papszTokens = CSLTokenizeStringComplex( pszCornerLoc, ",",
1337 0 : FALSE, FALSE );
1338 0 : if( CSLCount( papszTokens ) != 2 )
1339 0 : return;
1340 :
1341 0 : adfXY[iCorner*2+0] = CPLAtof( papszTokens[1] );
1342 0 : adfXY[iCorner*2+1] = CPLAtof( papszTokens[0] );
1343 :
1344 0 : if( adfXY[iCorner*2+0] < -360 || adfXY[iCorner*2+0] > 360
1345 0 : || adfXY[iCorner*2+1] < -90 || adfXY[iCorner*2+1] > 90 )
1346 0 : bLLPossible = FALSE;
1347 :
1348 0 : CSLDestroy( papszTokens );
1349 : }
1350 :
1351 : /* -------------------------------------------------------------------- */
1352 : /* Does this look like nice clean "northup" lat/long data? */
1353 : /* -------------------------------------------------------------------- */
1354 0 : if( adfXY[0*2+0] == adfXY[2*2+0] && adfXY[0*2+1] == adfXY[1*2+1]
1355 : && bLLPossible )
1356 : {
1357 0 : bHasGeoTransform = TRUE;
1358 0 : adfGeoTransform[0] = adfXY[0*2+0];
1359 0 : adfGeoTransform[1] = (adfXY[1*2+0] - adfXY[0*2+0]) / nRasterXSize;
1360 0 : adfGeoTransform[2] = 0.0;
1361 0 : adfGeoTransform[3] = adfXY[0*2+1];
1362 0 : adfGeoTransform[4] = 0.0;
1363 0 : adfGeoTransform[5] = (adfXY[2*2+1] - adfXY[0*2+1]) / nRasterYSize;
1364 :
1365 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
1366 0 : CPLFree( pszProjection );
1367 0 : oSRS.exportToWkt( &pszProjection );
1368 : }
1369 :
1370 : /* -------------------------------------------------------------------- */
1371 : /* Can we find the USGS Projection Parameters? */
1372 : /* -------------------------------------------------------------------- */
1373 0 : int bGotGCTPProjection = FALSE;
1374 0 : int iSDSIndex = FAIL, iSDS = FAIL;
1375 : const char *mapProjection = CSLFetchNameValue( papszGlobalMetadata,
1376 0 : "mapProjection" );
1377 :
1378 0 : if( mapProjection )
1379 0 : iSDSIndex = SDnametoindex( hSD, mapProjection );
1380 :
1381 0 : if( iSDSIndex != FAIL )
1382 0 : iSDS = SDselect( hSD, iSDSIndex );
1383 :
1384 0 : if( iSDS != FAIL )
1385 : {
1386 : char szName[HDF4_SDS_MAXNAMELEN];
1387 : int32 iRank, iNumType, nAttrs;
1388 : int32 aiDimSizes[H4_MAX_VAR_DIMS];
1389 :
1390 : double adfGCTP[29];
1391 : int32 aiStart[H4_MAX_NC_DIMS], aiEdges[H4_MAX_NC_DIMS];
1392 :
1393 0 : aiStart[0] = 0;
1394 0 : aiEdges[0] = 29;
1395 :
1396 0 : if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType,
1397 : &nAttrs) == 0
1398 : && iNumType == DFNT_FLOAT64
1399 : && iRank == 1
1400 0 : && aiDimSizes[0] >= 29
1401 : && SDreaddata( iSDS, aiStart, NULL, aiEdges, adfGCTP ) == 0
1402 0 : && oSRS.importFromUSGS( (long) adfGCTP[1], (long) adfGCTP[2],
1403 : adfGCTP+4,
1404 0 : (long) adfGCTP[3] ) == OGRERR_NONE )
1405 : {
1406 : CPLDebug( "HDF4Image", "GCTP Parms = %g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g",
1407 : adfGCTP[0],
1408 : adfGCTP[1],
1409 : adfGCTP[2],
1410 : adfGCTP[3],
1411 : adfGCTP[4],
1412 : adfGCTP[5],
1413 : adfGCTP[6],
1414 : adfGCTP[7],
1415 : adfGCTP[8],
1416 : adfGCTP[9],
1417 : adfGCTP[10],
1418 : adfGCTP[11],
1419 : adfGCTP[12],
1420 : adfGCTP[13],
1421 : adfGCTP[14],
1422 : adfGCTP[15],
1423 : adfGCTP[16],
1424 : adfGCTP[17],
1425 : adfGCTP[18],
1426 : adfGCTP[19],
1427 : adfGCTP[20],
1428 : adfGCTP[21],
1429 : adfGCTP[22],
1430 : adfGCTP[23],
1431 : adfGCTP[24],
1432 : adfGCTP[25],
1433 : adfGCTP[26],
1434 : adfGCTP[27],
1435 0 : adfGCTP[28] );
1436 :
1437 0 : CPLFree( pszProjection );
1438 0 : oSRS.exportToWkt( &pszProjection );
1439 0 : bGotGCTPProjection = TRUE;
1440 : }
1441 :
1442 0 : SDendaccess(iSDS);
1443 : }
1444 :
1445 : /* -------------------------------------------------------------------- */
1446 : /* If we derived a GCTP based projection, then we need to */
1447 : /* transform the lat/long corners into this projection and use */
1448 : /* them to establish the geotransform. */
1449 : /* -------------------------------------------------------------------- */
1450 0 : if( bLLPossible && bGotGCTPProjection )
1451 : {
1452 : double dfULX, dfULY, dfLRX, dfLRY;
1453 0 : OGRSpatialReference oWGS84;
1454 :
1455 0 : oWGS84.SetWellKnownGeogCS( "WGS84" );
1456 :
1457 : OGRCoordinateTransformation *poCT =
1458 0 : OGRCreateCoordinateTransformation( &oWGS84, &oSRS );
1459 :
1460 0 : dfULX = adfXY[0*2+0];
1461 0 : dfULY = adfXY[0*2+1];
1462 :
1463 0 : dfLRX = adfXY[3*2+0];
1464 0 : dfLRY = adfXY[3*2+1];
1465 :
1466 0 : if( poCT->Transform( 1, &dfULX, &dfULY )
1467 0 : && poCT->Transform( 1, &dfLRX, &dfLRY ) )
1468 : {
1469 0 : bHasGeoTransform = TRUE;
1470 0 : adfGeoTransform[0] = dfULX;
1471 0 : adfGeoTransform[1] = (dfLRX - dfULX) / nRasterXSize;
1472 0 : adfGeoTransform[2] = 0.0;
1473 0 : adfGeoTransform[3] = dfULY;
1474 0 : adfGeoTransform[4] = 0.0;
1475 0 : adfGeoTransform[5] = (dfLRY - dfULY) / nRasterYSize;
1476 : }
1477 :
1478 0 : delete poCT;
1479 : }
1480 : }
1481 :
1482 : /************************************************************************/
1483 : /* CaptureCoastwatchGCTPInfo() */
1484 : /************************************************************************/
1485 :
1486 : /* Example Metadata from:
1487 :
1488 : http://coastwatch.noaa.gov/interface/most_recent.php?sensor=MODIS&product=chlorNASA
1489 :
1490 : Definitions at:
1491 : http://coastwatch.noaa.gov/cw_form_hdf.html
1492 :
1493 : Metadata:
1494 : satellite=Aqua
1495 : sensor=MODIS
1496 : origin=USDOC/NOAA/NESDIS CoastWatch
1497 : history=PGE01:4.1.12;PGE02:4.3.1.12;SeaDAS Version ?.?, MSl12 4.0.2, Linux 2.4.21-27.0.1.EL
1498 : cwregister GulfOfMexicoSinusoidal.hdf MODSCW.P2005023.1835.swath09.hdf MODSCW.P2005023.1835.GM16.mapped09.hdf
1499 : cwgraphics MODSCW.P2005023.1835.GM16.closest.hdf
1500 : cwmath --template chlor_a --expr chlor_a=select(and(l2_flags,514)!=0,nan,chlor_a) /data/aps/browse/lvl3/seadas/coastwatch/hdf/MODSCW_P2005023_1835_GM16_closest.hdf /data/aps/browse/lvl3/seadas/coastwatch/maskhdf/MODSCW_P2005023_1835_GM16_closest_chlora.hdf
1501 : cwmath --template latitude --expr latitude=latitude /data/aps/browse/lvl3/seadas/coastwatch/hdf/MODSCW_P2005023_1835_GM16_closest.hdf /data/aps/browse/lvl3/seadas/coastwatch/maskhdf/MODSCW_P2005023_1835_GM16_closest_chlora.hdf
1502 : cwmath --template longitude --expr longitude=longitude /data/aps/browse/lvl3/seadas/coastwatch/hdf/MODSCW_P2005023_1835_GM16_closest.hdf /data/aps/browse/lvl3/seadas/coastwatch/maskhdf/MODSCW_P2005023_1835_GM16_closest_chlora.hdf
1503 : cwhdf_version=3.2
1504 : pass_type=day
1505 : pass_date=12806
1506 : start_time=66906
1507 : temporal_extent=298
1508 : projection_type=mapped
1509 : projection=Sinusoidal
1510 : gctp_sys=16
1511 : gctp_zone=62
1512 : gctp_parm=6378137, 0, 0, 0, -89000000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
1513 : gctp_datum=12
1514 : et_affine=0, -1008.74836097881, 1008.74836097881, 0, -953126.102425113, 3447041.10282512
1515 : rows=1540
1516 : cols=2000
1517 : polygon_latitude=31, 31, 31, 31, 31, 27.5095879249529, 24.0191758499058, 20.5287637748587, 17.0383516998116, 17.0383516998116, 17.0383516998116, 17.0383516998116, 17.0383516998116, 20.5287637748587, 24.0191758499058, 27.5095879249529, 31
1518 : polygon_longitude=-99, -93.7108573344442, -88.4217146688883, -83.1325720033325, -77.8434293377767, -78.217853417453, -78.5303805448579, -78.7884829057512, -78.9979508907244, -83.7397542896832, -88.481557688642, -93.2233610876007, -97.9651644865595, -98.1529175079091, -98.3842631146439, -98.664391423662, -99
1519 : orbit_type=ascending
1520 : raster_type=RasterPixelIsArea
1521 : swath_sync_lines=1
1522 :
1523 : */
1524 :
1525 0 : void HDF4ImageDataset::CaptureCoastwatchGCTPInfo()
1526 :
1527 : {
1528 0 : if( CSLFetchNameValue( papszGlobalMetadata, "gctp_sys" ) == NULL
1529 : || CSLFetchNameValue( papszGlobalMetadata, "gctp_zone" ) == NULL
1530 : || CSLFetchNameValue( papszGlobalMetadata, "gctp_parm" ) == NULL
1531 : || CSLFetchNameValue( papszGlobalMetadata, "gctp_datum" ) == NULL
1532 : || CSLFetchNameValue( papszGlobalMetadata, "et_affine" ) == NULL )
1533 0 : return;
1534 :
1535 : /* -------------------------------------------------------------------- */
1536 : /* Grab USGS/GCTP Parameters. */
1537 : /* -------------------------------------------------------------------- */
1538 : int nSys, nZone, nDatum, iParm;
1539 : double adfParms[15];
1540 : char **papszTokens;
1541 :
1542 0 : nSys = atoi( CSLFetchNameValue( papszGlobalMetadata, "gctp_sys" ) );
1543 0 : nZone = atoi( CSLFetchNameValue( papszGlobalMetadata, "gctp_zone" ) );
1544 0 : nDatum = atoi( CSLFetchNameValue( papszGlobalMetadata, "gctp_datum" ) );
1545 :
1546 : papszTokens = CSLTokenizeStringComplex(
1547 : CSLFetchNameValue( papszGlobalMetadata, "gctp_parm" ), ",",
1548 0 : FALSE, FALSE );
1549 0 : if( CSLCount(papszTokens) < 15 )
1550 0 : return;
1551 :
1552 0 : for( iParm = 0; iParm < 15; iParm++ )
1553 0 : adfParms[iParm] = CPLAtof( papszTokens[iParm] );
1554 0 : CSLDestroy( papszTokens );
1555 :
1556 : /* -------------------------------------------------------------------- */
1557 : /* Convert into an SRS. */
1558 : /* -------------------------------------------------------------------- */
1559 :
1560 0 : if( oSRS.importFromUSGS( nSys, nZone, adfParms, nDatum ) != OGRERR_NONE )
1561 0 : return;
1562 :
1563 0 : CPLFree( pszProjection );
1564 0 : oSRS.exportToWkt( &pszProjection );
1565 :
1566 : /* -------------------------------------------------------------------- */
1567 : /* Capture the affine transform info. */
1568 : /* -------------------------------------------------------------------- */
1569 :
1570 : papszTokens = CSLTokenizeStringComplex(
1571 : CSLFetchNameValue( papszGlobalMetadata, "et_affine" ), ",",
1572 0 : FALSE, FALSE );
1573 0 : if( CSLCount(papszTokens) != 6 )
1574 0 : return;
1575 :
1576 : // We don't seem to have proper ef_affine docs so I don't
1577 : // know which of these two coefficients goes where.
1578 0 : if( CPLAtof(papszTokens[0]) != 0.0 || CPLAtof(papszTokens[3]) != 0.0 )
1579 0 : return;
1580 :
1581 0 : bHasGeoTransform = TRUE;
1582 0 : adfGeoTransform[0] = CPLAtof( papszTokens[4] );
1583 0 : adfGeoTransform[1] = CPLAtof( papszTokens[2] );
1584 0 : adfGeoTransform[2] = 0.0;
1585 0 : adfGeoTransform[3] = CPLAtof( papszTokens[5] );
1586 0 : adfGeoTransform[4] = 0.0;
1587 0 : adfGeoTransform[5] = CPLAtof( papszTokens[1] );
1588 :
1589 : // Middle of pixel adjustment.
1590 0 : adfGeoTransform[0] -= adfGeoTransform[1] * 0.5;
1591 0 : adfGeoTransform[3] -= adfGeoTransform[5] * 0.5;
1592 : }
1593 :
1594 : /************************************************************************/
1595 : /* GetImageDimensions() */
1596 : /************************************************************************/
1597 :
1598 9 : void HDF4ImageDataset::GetImageDimensions( char *pszDimList )
1599 : {
1600 : char **papszDimList = CSLTokenizeString2( pszDimList, ",",
1601 9 : CSLT_HONOURSTRINGS );
1602 9 : int i, nDimCount = CSLCount( papszDimList );
1603 :
1604 : // TODO: check whether nDimCount is > 1 and do something if it isn't.
1605 :
1606 : // Search for the "Band" word in the name of dimension
1607 : // or take the first one as a number of bands
1608 9 : if ( iRank == 2 )
1609 9 : nBandCount = 1;
1610 : else
1611 : {
1612 0 : for ( i = 0; i < nDimCount; i++ )
1613 : {
1614 0 : if ( strstr( papszDimList[i], "band" ) )
1615 : {
1616 0 : iBandDim = i;
1617 0 : nBandCount = aiDimSizes[i];
1618 : // Handle 4D datasets
1619 0 : if ( iRank > 3 && i < nDimCount - 1 )
1620 : {
1621 : // FIXME: is there a better way to search for
1622 : // the 4th dimension?
1623 0 : i4Dim = i + 1;
1624 0 : nBandCount *= aiDimSizes[i4Dim];
1625 : }
1626 0 : break;
1627 : }
1628 : }
1629 : }
1630 :
1631 : // Search for the starting "X" and "Y" in the names or take
1632 : // the last two dimensions as X and Y sizes
1633 9 : iXDim = nDimCount - 1;
1634 9 : iYDim = nDimCount - 2;
1635 :
1636 27 : for ( i = 0; i < nDimCount; i++ )
1637 : {
1638 27 : if ( EQUALN( papszDimList[i], "X", 1 ) && iBandDim != i )
1639 9 : iXDim = i;
1640 9 : else if ( EQUALN( papszDimList[i], "Y", 1 ) && iBandDim != i )
1641 9 : iYDim = i;
1642 : }
1643 :
1644 : // If didn't get a band dimension yet, but have an extra
1645 : // dimension, use it as the band dimension.
1646 9 : if ( iRank > 2 && iBandDim == -1 )
1647 : {
1648 0 : if( iXDim != 0 && iYDim != 0 )
1649 0 : iBandDim = 0;
1650 0 : else if( iXDim != 1 && iYDim != 1 )
1651 0 : iBandDim = 1;
1652 0 : else if( iXDim != 2 && iYDim != 2 )
1653 0 : iBandDim = 2;
1654 :
1655 0 : nBandCount = aiDimSizes[iBandDim];
1656 : }
1657 :
1658 9 : CSLDestroy( papszDimList );
1659 9 : }
1660 :
1661 : /************************************************************************/
1662 : /* GetSwatAttrs() */
1663 : /************************************************************************/
1664 :
1665 0 : void HDF4ImageDataset::GetSwatAttrs( int32 hSW )
1666 : {
1667 : /* -------------------------------------------------------------------- */
1668 : /* At the start we will fetch the global HDF attributes. */
1669 : /* -------------------------------------------------------------------- */
1670 : int32 hDummy;
1671 :
1672 0 : EHidinfo( hHDF4, &hDummy, &hSD );
1673 0 : ReadGlobalAttributes( hSD );
1674 0 : papszLocalMetadata = CSLDuplicate( papszGlobalMetadata );
1675 :
1676 : /* -------------------------------------------------------------------- */
1677 : /* Fetch the esoteric HDF-EOS attributes then. */
1678 : /* -------------------------------------------------------------------- */
1679 0 : int32 nStrBufSize = 0;
1680 :
1681 0 : if ( SWinqattrs( hSW, NULL, &nStrBufSize ) > 0 && nStrBufSize > 0 )
1682 : {
1683 : char *pszAttrList;
1684 : char **papszAttributes;
1685 : int i, nAttrs;
1686 :
1687 0 : pszAttrList = (char *)CPLMalloc( nStrBufSize + 1 );
1688 0 : SWinqattrs( hSW, pszAttrList, &nStrBufSize );
1689 :
1690 : #ifdef DEBUG
1691 : CPLDebug( "HDF4Image", "List of attributes in swath \"%s\": %s",
1692 0 : pszFieldName, pszAttrList );
1693 : #endif
1694 :
1695 : papszAttributes = CSLTokenizeString2( pszAttrList, ",",
1696 0 : CSLT_HONOURSTRINGS );
1697 0 : nAttrs = CSLCount( papszAttributes );
1698 0 : for ( i = 0; i < nAttrs; i++ )
1699 : {
1700 : int32 iNumType, nValues;
1701 0 : void *pData = NULL;
1702 0 : char *pszTemp = NULL;
1703 :
1704 0 : SWattrinfo( hSW, papszAttributes[i], &iNumType, &nValues );
1705 :
1706 0 : if ( iNumType == DFNT_CHAR8 || iNumType == DFNT_UCHAR8 )
1707 0 : pData = CPLMalloc( (nValues + 1) * GetDataTypeSize(iNumType) );
1708 : else
1709 0 : pData = CPLMalloc( nValues * GetDataTypeSize(iNumType) );
1710 :
1711 0 : SWreadattr( hSW, papszAttributes[i], pData );
1712 :
1713 0 : if ( iNumType == DFNT_CHAR8 || iNumType == DFNT_UCHAR8 )
1714 : {
1715 0 : ((char *)pData)[nValues] = '\0';
1716 : papszLocalMetadata = CSLAddNameValue( papszLocalMetadata,
1717 0 : papszAttributes[i],
1718 0 : (const char *) pData );
1719 : }
1720 : else
1721 : {
1722 : pszTemp = SPrintArray( GetDataType(iNumType), pData,
1723 0 : nValues, ", " );
1724 : papszLocalMetadata = CSLAddNameValue( papszLocalMetadata,
1725 0 : papszAttributes[i],
1726 0 : pszTemp );
1727 0 : if ( pszTemp )
1728 0 : CPLFree( pszTemp );
1729 : }
1730 :
1731 0 : if ( pData )
1732 0 : CPLFree( pData );
1733 :
1734 : }
1735 :
1736 0 : CSLDestroy( papszAttributes );
1737 0 : CPLFree( pszAttrList );
1738 : }
1739 :
1740 : /* -------------------------------------------------------------------- */
1741 : /* After fetching HDF-EOS specific stuff we will read the generic */
1742 : /* HDF attributes and append them to the list of metadata. */
1743 : /* -------------------------------------------------------------------- */
1744 : int32 iSDS;
1745 0 : if ( SWsdid(hSW, pszFieldName, &iSDS) != -1 )
1746 : {
1747 : int32 iRank, iNumType, iAttribute, nAttrs, nValues;
1748 : char szName[HDF4_SDS_MAXNAMELEN];
1749 : int32 aiDimSizes[H4_MAX_VAR_DIMS];
1750 :
1751 0 : if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType,
1752 : &nAttrs) == 0 )
1753 : {
1754 0 : for ( iAttribute = 0; iAttribute < nAttrs; iAttribute++ )
1755 : {
1756 : char szAttrName[H4_MAX_NC_NAME];
1757 : SDattrinfo( iSDS, iAttribute, szAttrName,
1758 0 : &iNumType, &nValues );
1759 : papszLocalMetadata =
1760 : TranslateHDF4Attributes( iSDS, iAttribute,
1761 : szAttrName, iNumType,
1762 0 : nValues, papszLocalMetadata );
1763 : }
1764 : }
1765 : }
1766 :
1767 : /* -------------------------------------------------------------------- */
1768 : /* Finally make the whole list visible. */
1769 : /* -------------------------------------------------------------------- */
1770 0 : SetMetadata( papszLocalMetadata );
1771 0 : }
1772 :
1773 : /************************************************************************/
1774 : /* GetGridAttrs() */
1775 : /************************************************************************/
1776 :
1777 9 : void HDF4ImageDataset::GetGridAttrs( int32 hGD )
1778 : {
1779 : /* -------------------------------------------------------------------- */
1780 : /* At the start we will fetch the global HDF attributes. */
1781 : /* -------------------------------------------------------------------- */
1782 : int32 hDummy;
1783 :
1784 9 : EHidinfo( hHDF4, &hDummy, &hSD );
1785 9 : ReadGlobalAttributes( hSD );
1786 9 : papszLocalMetadata = CSLDuplicate( papszGlobalMetadata );
1787 :
1788 : /* -------------------------------------------------------------------- */
1789 : /* Fetch the esoteric HDF-EOS attributes then. */
1790 : /* -------------------------------------------------------------------- */
1791 9 : int32 nStrBufSize = 0;
1792 :
1793 9 : if ( GDinqattrs( hGD, NULL, &nStrBufSize ) > 0 && nStrBufSize > 0 )
1794 : {
1795 : char *pszAttrList;
1796 : char **papszAttributes;
1797 : int i, nAttrs;
1798 :
1799 0 : pszAttrList = (char *)CPLMalloc( nStrBufSize + 1 );
1800 0 : GDinqattrs( hGD, pszAttrList, &nStrBufSize );
1801 : #ifdef DEBUG
1802 : CPLDebug( "HDF4Image", "List of attributes in grid %s: %s",
1803 0 : pszFieldName, pszAttrList );
1804 : #endif
1805 : papszAttributes = CSLTokenizeString2( pszAttrList, ",",
1806 0 : CSLT_HONOURSTRINGS );
1807 0 : nAttrs = CSLCount( papszAttributes );
1808 0 : for ( i = 0; i < nAttrs; i++ )
1809 : {
1810 : int32 iNumType, nValues;
1811 0 : void *pData = NULL;
1812 0 : char *pszTemp = NULL;
1813 :
1814 0 : GDattrinfo( hGD, papszAttributes[i], &iNumType, &nValues );
1815 :
1816 0 : if ( iNumType == DFNT_CHAR8 || iNumType == DFNT_UCHAR8 )
1817 0 : pData = CPLMalloc( (nValues + 1) * GetDataTypeSize(iNumType) );
1818 : else
1819 0 : pData = CPLMalloc( nValues * GetDataTypeSize(iNumType) );
1820 :
1821 0 : GDreadattr( hGD, papszAttributes[i], pData );
1822 :
1823 0 : if ( iNumType == DFNT_CHAR8 || iNumType == DFNT_UCHAR8 )
1824 : {
1825 0 : ((char *)pData)[nValues] = '\0';
1826 : papszLocalMetadata = CSLAddNameValue( papszLocalMetadata,
1827 0 : papszAttributes[i],
1828 0 : (const char *) pData );
1829 : }
1830 : else
1831 : {
1832 : pszTemp = SPrintArray( GetDataType(iNumType), pData,
1833 0 : nValues, ", " );
1834 : papszLocalMetadata = CSLAddNameValue( papszLocalMetadata,
1835 0 : papszAttributes[i],
1836 0 : pszTemp );
1837 0 : if ( pszTemp )
1838 0 : CPLFree( pszTemp );
1839 : }
1840 :
1841 0 : if ( pData )
1842 0 : CPLFree( pData );
1843 :
1844 : }
1845 :
1846 0 : CSLDestroy( papszAttributes );
1847 0 : CPLFree( pszAttrList );
1848 : }
1849 :
1850 : /* -------------------------------------------------------------------- */
1851 : /* After fetching HDF-EOS specific stuff we will read the generic */
1852 : /* HDF attributes and append them to the list of metadata. */
1853 : /* -------------------------------------------------------------------- */
1854 : int32 iSDS;
1855 9 : if ( GDsdid(hGD, pszFieldName, &iSDS) != -1 )
1856 : {
1857 : int32 iRank, iNumType, iAttribute, nAttrs, nValues;
1858 : char szName[HDF4_SDS_MAXNAMELEN];
1859 : int32 aiDimSizes[H4_MAX_VAR_DIMS];
1860 :
1861 9 : if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType,
1862 : &nAttrs) == 0 )
1863 : {
1864 74 : for ( iAttribute = 0; iAttribute < nAttrs; iAttribute++ )
1865 : {
1866 : char szAttrName[H4_MAX_NC_NAME];
1867 : SDattrinfo( iSDS, iAttribute, szAttrName,
1868 65 : &iNumType, &nValues );
1869 : papszLocalMetadata =
1870 : TranslateHDF4Attributes( iSDS, iAttribute,
1871 : szAttrName, iNumType,
1872 65 : nValues, papszLocalMetadata );
1873 : }
1874 : }
1875 : }
1876 :
1877 : /* -------------------------------------------------------------------- */
1878 : /* Finally make the whole list visible. */
1879 : /* -------------------------------------------------------------------- */
1880 9 : SetMetadata( papszLocalMetadata );
1881 9 : }
1882 :
1883 : /************************************************************************/
1884 : /* ProcessModisSDSGeolocation() */
1885 : /* */
1886 : /* Recognise latitude and longitude geolocation arrays in */
1887 : /* simple SDS datasets like: */
1888 : /* */
1889 : /* download.osgeo.org/gdal/data/hdf4/A2006005182000.L2_LAC_SST.x.hdf */
1890 : /* */
1891 : /* As reported in ticket #1895. */
1892 : /* */
1893 : /* Note that we don't check that the dimensions of the latitude */
1894 : /* and longitude exactly match the dimensions of the basedata, */
1895 : /* though we ought to. */
1896 : /************************************************************************/
1897 :
1898 256 : void HDF4ImageDataset::ProcessModisSDSGeolocation(void)
1899 :
1900 : {
1901 256 : int iDSIndex, iXIndex=-1, iYIndex=-1;
1902 :
1903 : // No point in assigning geolocation to the geolocation SDSes themselves.
1904 256 : if( EQUAL(szName,"longitude") || EQUAL(szName,"latitude") )
1905 0 : return;
1906 :
1907 256 : if (nRasterYSize == 1)
1908 0 : return;
1909 :
1910 : /* -------------------------------------------------------------------- */
1911 : /* Scan for latitude and longitude sections. */
1912 : /* -------------------------------------------------------------------- */
1913 : int32 nDatasets, nAttributes;
1914 :
1915 256 : if ( SDfileinfo( hSD, &nDatasets, &nAttributes ) != 0 )
1916 0 : return;
1917 :
1918 672 : for( iDSIndex = 0; iDSIndex < nDatasets; iDSIndex++ )
1919 : {
1920 : int32 iRank, iNumType, nAttrs, iSDS;
1921 : char szName[HDF4_SDS_MAXNAMELEN];
1922 : int32 aiDimSizes[H4_MAX_VAR_DIMS];
1923 :
1924 416 : iSDS = SDselect( hSD, iDSIndex );
1925 :
1926 416 : if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType,
1927 : &nAttrs) == 0 )
1928 : {
1929 416 : if( EQUAL(szName,"latitude") )
1930 2 : iYIndex = iDSIndex;
1931 :
1932 416 : if( EQUAL(szName,"longitude") )
1933 2 : iXIndex = iDSIndex;
1934 : }
1935 :
1936 416 : SDendaccess(iSDS);
1937 : }
1938 :
1939 256 : if( iXIndex == -1 || iYIndex == -1 )
1940 254 : return;
1941 :
1942 : /* -------------------------------------------------------------------- */
1943 : /* We found geolocation information. Record it as metadata. */
1944 : /* -------------------------------------------------------------------- */
1945 2 : CPLString osWrk;
1946 :
1947 2 : SetMetadataItem( "SRS", SRS_WKT_WGS84, "GEOLOCATION" );
1948 :
1949 : osWrk.Printf( "HDF4_SDS:UNKNOWN:\"%s\":%d",
1950 2 : pszFilename, iXIndex );
1951 2 : SetMetadataItem( "X_DATASET", osWrk, "GEOLOCATION" );
1952 2 : SetMetadataItem( "X_BAND", "1" , "GEOLOCATION" );
1953 :
1954 : osWrk.Printf( "HDF4_SDS:UNKNOWN:\"%s\":%d",
1955 2 : pszFilename, iYIndex );
1956 2 : SetMetadataItem( "Y_DATASET", osWrk, "GEOLOCATION" );
1957 2 : SetMetadataItem( "Y_BAND", "1" , "GEOLOCATION" );
1958 :
1959 2 : SetMetadataItem( "PIXEL_OFFSET", "0", "GEOLOCATION" );
1960 2 : SetMetadataItem( "PIXEL_STEP", "1", "GEOLOCATION" );
1961 :
1962 2 : SetMetadataItem( "LINE_OFFSET", "0", "GEOLOCATION" );
1963 2 : SetMetadataItem( "LINE_STEP", "1", "GEOLOCATION" );
1964 : }
1965 :
1966 : /************************************************************************/
1967 : /* ProcessSwathGeolocation() */
1968 : /* */
1969 : /* Handle the swath geolocation data for a swath. Attach */
1970 : /* geolocation metadata corresponding to it (if there is no */
1971 : /* lattice), and also attach it as GCPs. This is only invoked */
1972 : /* for EOS_SWATH, not EOS_SWATH_GEOL datasets. */
1973 : /************************************************************************/
1974 :
1975 0 : int HDF4ImageDataset::ProcessSwathGeolocation( int32 hSW, char **papszDimList )
1976 : {
1977 0 : char szXGeo[8192] = "";
1978 0 : char szYGeo[8192] = "";
1979 0 : char szPixel[8192]= "";
1980 0 : char szLine[8192] = "";
1981 : int32 iWrkNumType;
1982 0 : void *pLat = NULL, *pLong = NULL;
1983 0 : void *pLatticeX = NULL, *pLatticeY = NULL;
1984 0 : int32 iLatticeType, iLatticeDataSize = 0, iRank;
1985 0 : int32 nLatCount = 0, nLongCount = 0;
1986 0 : int32 nXPoints=0, nYPoints=0;
1987 : int32 nStrBufSize;
1988 0 : int i, j, iDataSize = 0, iPixelDim=-1,iLineDim=-1, iLongDim=-1, iLatDim=-1;
1989 0 : int32 *paiRank = NULL, *paiNumType = NULL,
1990 0 : *paiOffset = NULL, *paiIncrement = NULL;
1991 :
1992 : /* -------------------------------------------------------------------- */
1993 : /* Determine a product name. */
1994 : /* -------------------------------------------------------------------- */
1995 : const char *pszProduct =
1996 0 : CSLFetchNameValue( papszLocalMetadata, "SHORTNAME" );
1997 :
1998 0 : HDF4EOSProduct eProduct = PROD_UNKNOWN;
1999 0 : if ( pszProduct )
2000 : {
2001 0 : if ( EQUALN(pszProduct, "ASTL1A", 6) )
2002 0 : eProduct = PROD_ASTER_L1A;
2003 0 : else if ( EQUALN(pszProduct, "ASTL1B", 6) )
2004 0 : eProduct = PROD_ASTER_L1B;
2005 0 : else if ( EQUALN(pszProduct, "AST_04", 6)
2006 : || EQUALN(pszProduct, "AST_05", 6)
2007 : || EQUALN(pszProduct, "AST_06", 6)
2008 : || EQUALN(pszProduct, "AST_07", 6)
2009 : || EQUALN(pszProduct, "AST_08", 6)
2010 : || EQUALN(pszProduct, "AST_09", 6)
2011 : || EQUALN(pszProduct, "AST13", 5)
2012 : || EQUALN(pszProduct, "AST3", 4) )
2013 0 : eProduct = PROD_ASTER_L2;
2014 0 : else if ( EQUALN(pszProduct, "AST14", 5) )
2015 0 : eProduct = PROD_ASTER_L3;
2016 0 : else if ( EQUALN(pszProduct, "MOD02", 5)
2017 : || EQUALN(pszProduct, "MYD02", 5) )
2018 0 : eProduct = PROD_MODIS_L1B;
2019 0 : else if ( EQUALN(pszProduct, "MOD07_L2", 8) )
2020 0 : eProduct = PROD_MODIS_L2;
2021 : }
2022 :
2023 : /* -------------------------------------------------------------------- */
2024 : /* Read names of geolocation fields and corresponding */
2025 : /* geolocation maps. */
2026 : /* -------------------------------------------------------------------- */
2027 0 : int32 nDataFields = SWnentries( hSW, HDFE_NENTGFLD, &nStrBufSize );
2028 0 : char *pszGeoList = (char *)CPLMalloc( nStrBufSize + 1 );
2029 0 : paiRank = (int32 *)CPLMalloc( nDataFields * sizeof(int32) );
2030 0 : paiNumType = (int32 *)CPLMalloc( nDataFields * sizeof(int32) );
2031 :
2032 0 : if ( nDataFields !=
2033 : SWinqgeofields(hSW, pszGeoList, paiRank, paiNumType) )
2034 : {
2035 : CPLDebug( "HDF4Image",
2036 : "Can't get the list of geolocation fields in swath \"%s\"",
2037 0 : pszSubdatasetName );
2038 : }
2039 :
2040 : #ifdef DEBUG
2041 : else
2042 : {
2043 : char *pszTmp;
2044 : CPLDebug( "HDF4Image",
2045 : "Number of geolocation fields in swath \"%s\": %ld",
2046 0 : pszSubdatasetName, (long)nDataFields );
2047 : CPLDebug( "HDF4Image",
2048 : "List of geolocation fields in swath \"%s\": %s",
2049 0 : pszSubdatasetName, pszGeoList );
2050 : pszTmp = SPrintArray( GDT_UInt32, paiRank,
2051 0 : nDataFields, "," );
2052 : CPLDebug( "HDF4Image",
2053 0 : "Geolocation fields ranks: %s", pszTmp );
2054 0 : CPLFree( pszTmp );
2055 : }
2056 : #endif
2057 :
2058 : // Read geolocation data
2059 0 : int32 nDimMaps = SWnentries( hSW, HDFE_NENTMAP, &nStrBufSize );
2060 0 : if ( nDimMaps <= 0 )
2061 : {
2062 :
2063 : #ifdef DEBUG
2064 : CPLDebug( "HDF4Image", "No geolocation maps in swath \"%s\"",
2065 0 : pszSubdatasetName );
2066 : CPLDebug( "HDF4Image",
2067 : "Suppose one-to-one mapping. X field is \"%s\", Y field is \"%s\"",
2068 0 : papszDimList[iXDim], papszDimList[iYDim] );
2069 : #endif
2070 :
2071 0 : strncpy( szPixel, papszDimList[iXDim], 8192 );
2072 0 : strncpy( szLine, papszDimList[iYDim], 8192 );
2073 0 : strncpy( szXGeo, papszDimList[iXDim], 8192 );
2074 0 : strncpy( szYGeo, papszDimList[iYDim], 8192 );
2075 0 : paiOffset = (int32 *)CPLCalloc( 2, sizeof(int32) );
2076 0 : paiIncrement = (int32 *)CPLCalloc( 2, sizeof(int32) );
2077 0 : paiOffset[0] = paiOffset[1] = 0;
2078 0 : paiIncrement[0] = paiIncrement[1] = 1;
2079 : }
2080 : else
2081 : {
2082 0 : char *pszDimMaps = (char *)CPLMalloc( nStrBufSize + 1 );
2083 0 : paiOffset = (int32 *)CPLCalloc( nDimMaps, sizeof(int32) );
2084 0 : paiIncrement = (int32 *)CPLCalloc( nDimMaps, sizeof(int32) );
2085 :
2086 0 : *pszDimMaps = '\0';
2087 0 : if ( nDimMaps != SWinqmaps(hSW, pszDimMaps, paiOffset, paiIncrement) )
2088 : {
2089 : CPLDebug( "HDF4Image",
2090 : "Can't get the list of geolocation maps in swath \"%s\"",
2091 0 : pszSubdatasetName );
2092 : }
2093 :
2094 : #ifdef DEBUG
2095 : else
2096 : {
2097 : char *pszTmp;
2098 :
2099 : CPLDebug( "HDF4Image",
2100 : "List of geolocation maps in swath \"%s\": %s",
2101 0 : pszSubdatasetName, pszDimMaps );
2102 : pszTmp = SPrintArray( GDT_Int32, paiOffset,
2103 0 : nDimMaps, "," );
2104 : CPLDebug( "HDF4Image",
2105 0 : "Geolocation map offsets: %s", pszTmp );
2106 0 : CPLFree( pszTmp );
2107 : pszTmp = SPrintArray( GDT_Int32, paiIncrement,
2108 0 : nDimMaps, "," );
2109 : CPLDebug( "HDF4Image",
2110 0 : "Geolocation map increments: %s", pszTmp );
2111 0 : CPLFree( pszTmp );
2112 : }
2113 : #endif
2114 :
2115 : char **papszDimMap = CSLTokenizeString2( pszDimMaps, ",",
2116 0 : CSLT_HONOURSTRINGS );
2117 0 : int nDimMapCount = CSLCount(papszDimMap);
2118 :
2119 0 : for ( i = 0; i < nDimMapCount; i++ )
2120 : {
2121 0 : if ( strstr(papszDimMap[i], papszDimList[iXDim]) )
2122 : {
2123 0 : strncpy( szPixel, papszDimList[iXDim], 8192 );
2124 0 : strncpy( szXGeo, papszDimMap[i], 8192 );
2125 0 : char *pszTemp = strchr( szXGeo, '/' );
2126 0 : if ( pszTemp )
2127 0 : *pszTemp = '\0';
2128 : }
2129 0 : else if ( strstr(papszDimMap[i], papszDimList[iYDim]) )
2130 : {
2131 0 : strncpy( szLine, papszDimList[iYDim], 8192 );
2132 0 : strncpy( szYGeo, papszDimMap[i], 8192 );
2133 0 : char *pszTemp = strchr( szYGeo, '/' );
2134 0 : if ( pszTemp )
2135 0 : *pszTemp = '\0';
2136 : }
2137 : }
2138 :
2139 0 : CSLDestroy( papszDimMap );
2140 0 : CPLFree( pszDimMaps );
2141 : }
2142 :
2143 0 : if ( *szXGeo == 0 || *szYGeo == 0 )
2144 0 : return FALSE;
2145 :
2146 : /* -------------------------------------------------------------------- */
2147 : /* Read geolocation fields. */
2148 : /* -------------------------------------------------------------------- */
2149 0 : char szGeoDimList[8192] = "";
2150 : char **papszGeolocations = CSLTokenizeString2( pszGeoList, ",",
2151 0 : CSLT_HONOURSTRINGS );
2152 0 : int nGeolocationsCount = CSLCount( papszGeolocations );
2153 : int32 aiDimSizes[H4_MAX_VAR_DIMS];
2154 :
2155 0 : for ( i = 0; i < nGeolocationsCount; i++ )
2156 : {
2157 0 : char **papszGeoDimList = NULL;
2158 :
2159 0 : if ( SWfieldinfo( hSW, papszGeolocations[i], &iRank,
2160 : aiDimSizes, &iWrkNumType, szGeoDimList ) < 0 )
2161 : {
2162 :
2163 : CPLDebug( "HDF4Image",
2164 : "Can't read attributes of geolocation field \"%s\"",
2165 0 : papszGeolocations[i] );
2166 0 : return FALSE;
2167 : }
2168 :
2169 : CPLDebug( "HDF4Image",
2170 : "List of dimensions in geolocation field \"%s\": %s",
2171 0 : papszGeolocations[i], szGeoDimList );
2172 :
2173 : papszGeoDimList = CSLTokenizeString2( szGeoDimList,
2174 0 : ",", CSLT_HONOURSTRINGS );
2175 :
2176 0 : if( CSLCount(papszGeoDimList) > H4_MAX_VAR_DIMS
2177 : || CSLFindString( papszGeoDimList, szXGeo ) == -1
2178 : || CSLFindString( papszGeoDimList, szYGeo ) == -1 )
2179 : {
2180 0 : CSLDestroy( papszGeoDimList );
2181 0 : return FALSE;
2182 : }
2183 :
2184 0 : nXPoints = aiDimSizes[CSLFindString( papszGeoDimList, szXGeo )];
2185 0 : nYPoints = aiDimSizes[CSLFindString( papszGeoDimList, szYGeo )];
2186 :
2187 0 : if ( EQUAL(szPixel, papszDimList[iXDim]) )
2188 : {
2189 0 : iPixelDim = 1;
2190 0 : iLineDim = 0;
2191 : }
2192 : else
2193 : {
2194 0 : iPixelDim = 0;
2195 0 : iLineDim = 1;
2196 : }
2197 :
2198 0 : iDataSize = GetDataTypeSize( iWrkNumType );
2199 0 : if ( strstr( papszGeolocations[i], "Latitude" ) )
2200 : {
2201 0 : iLatDim = i;
2202 0 : nLatCount = nXPoints * nYPoints;
2203 0 : pLat = CPLMalloc( nLatCount * iDataSize );
2204 0 : if (SWreadfield( hSW, papszGeolocations[i], NULL,
2205 : NULL, NULL, (VOIDP)pLat ) < 0)
2206 : {
2207 : CPLDebug( "HDF4Image",
2208 : "Can't read geolocation field %s",
2209 0 : papszGeolocations[i]);
2210 0 : CPLFree( pLat );
2211 0 : pLat = NULL;
2212 : }
2213 : }
2214 0 : else if ( strstr( papszGeolocations[i], "Longitude" ) )
2215 : {
2216 0 : iLongDim = i;
2217 0 : nLongCount = nXPoints * nYPoints;
2218 0 : pLong = CPLMalloc( nLongCount * iDataSize );
2219 0 : if (SWreadfield( hSW, papszGeolocations[i], NULL,
2220 : NULL, NULL, (VOIDP)pLong ) < 0)
2221 : {
2222 : CPLDebug( "HDF4Image",
2223 : "Can't read geolocation field %s",
2224 0 : papszGeolocations[i]);
2225 0 : CPLFree( pLong );
2226 0 : pLong = NULL;
2227 : }
2228 : }
2229 :
2230 0 : CSLDestroy( papszGeoDimList );
2231 : }
2232 :
2233 : /* -------------------------------------------------------------------- */
2234 : /* Do we have a lattice table? */
2235 : /* -------------------------------------------------------------------- */
2236 0 : if (SWfieldinfo(hSW, "LatticePoint", &iRank, aiDimSizes,
2237 : &iLatticeType, szGeoDimList) == 0
2238 : && iRank == 3
2239 0 : && nXPoints == aiDimSizes[1]
2240 0 : && nYPoints == aiDimSizes[0]
2241 0 : && aiDimSizes[2] == 2 )
2242 : {
2243 : int32 iStart[H4_MAX_NC_DIMS], iEdges[H4_MAX_NC_DIMS];
2244 :
2245 : iLatticeDataSize =
2246 0 : GetDataTypeSize( iLatticeType );
2247 :
2248 0 : iStart[1] = 0;
2249 0 : iEdges[1] = nXPoints;
2250 :
2251 0 : iStart[0] = 0;
2252 0 : iEdges[0] = nYPoints;
2253 :
2254 0 : iStart[2] = 0;
2255 0 : iEdges[2] = 1;
2256 :
2257 0 : pLatticeX = CPLMalloc( nLatCount * iLatticeDataSize );
2258 0 : if (SWreadfield( hSW, "LatticePoint", iStart, NULL,
2259 : iEdges, (VOIDP)pLatticeX ) < 0)
2260 : {
2261 0 : CPLDebug( "HDF4Image", "Can't read lattice field" );
2262 0 : CPLFree( pLatticeX );
2263 0 : pLatticeX = NULL;
2264 : }
2265 :
2266 0 : iStart[2] = 1;
2267 0 : iEdges[2] = 1;
2268 :
2269 0 : pLatticeY = CPLMalloc( nLatCount * iLatticeDataSize );
2270 0 : if (SWreadfield( hSW, "LatticePoint", iStart, NULL,
2271 : iEdges, (VOIDP)pLatticeY ) < 0)
2272 : {
2273 0 : CPLDebug( "HDF4Image", "Can't read lattice field" );
2274 0 : CPLFree( pLatticeY );
2275 0 : pLatticeY = NULL;
2276 : }
2277 :
2278 : }
2279 :
2280 : /* -------------------------------------------------------------------- */
2281 : /* Determine whether to use no, partial or full GCPs. */
2282 : /* -------------------------------------------------------------------- */
2283 : const char *pszGEOL_AS_GCPS = CPLGetConfigOption( "GEOL_AS_GCPS",
2284 0 : "PARTIAL" );
2285 : int iGCPStepX, iGCPStepY;
2286 :
2287 0 : if( EQUAL(pszGEOL_AS_GCPS,"NONE") )
2288 : {
2289 0 : iGCPStepX = iGCPStepY = 0;
2290 : }
2291 0 : else if( EQUAL(pszGEOL_AS_GCPS,"FULL") )
2292 : {
2293 0 : iGCPStepX = iGCPStepY = 1;
2294 : }
2295 : else
2296 : {
2297 : // aim for 10x10 grid or so.
2298 0 : iGCPStepX = MAX(1,((nXPoints-1) / 11));
2299 0 : iGCPStepY = MAX(1,((nYPoints-1) / 11));
2300 : }
2301 :
2302 : /* -------------------------------------------------------------------- */
2303 : /* Fetch projection information for various datasets. */
2304 : /* -------------------------------------------------------------------- */
2305 0 : if ( nLatCount && nLongCount && nLatCount == nLongCount
2306 : && pLat && pLong )
2307 : {
2308 0 : CPLFree( pszGCPProjection );
2309 0 : pszGCPProjection = NULL;
2310 :
2311 : // ASTER Level 1A
2312 0 : if ( eProduct == PROD_ASTER_L1A )
2313 : {
2314 0 : pszGCPProjection = CPLStrdup( "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9108\"]],AXIS[\"Lat\",NORTH],AXIS[\"Long\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]" );
2315 : }
2316 :
2317 : // ASTER Level 1B, Level 2
2318 0 : else if ( eProduct == PROD_ASTER_L1B
2319 : || eProduct == PROD_ASTER_L2 )
2320 : {
2321 : // Constuct the metadata keys.
2322 : // A band number is taken from the field name.
2323 0 : const char *pszBand = strpbrk( pszFieldName, "0123456789" );
2324 :
2325 0 : if ( !pszBand )
2326 0 : pszBand = "";
2327 :
2328 : char *pszProjLine =
2329 0 : CPLStrdup(CPLSPrintf("MPMETHOD%s", pszBand));
2330 : char *pszParmsLine =
2331 : CPLStrdup(CPLSPrintf("PROJECTIONPARAMETERS%s",
2332 0 : pszBand));
2333 : char *pszZoneLine =
2334 : CPLStrdup(CPLSPrintf("UTMZONECODE%s",
2335 0 : pszBand));
2336 : char *pszEllipsoidLine =
2337 : CPLStrdup(CPLSPrintf("ELLIPSOIDANDDATUM%s",
2338 0 : pszBand));
2339 :
2340 : // Fetch projection related values from the
2341 : // metadata.
2342 : const char *pszProj =
2343 : CSLFetchNameValue( papszLocalMetadata,
2344 0 : pszProjLine );
2345 : const char *pszParms =
2346 : CSLFetchNameValue( papszLocalMetadata,
2347 0 : pszParmsLine );
2348 : const char *pszZone =
2349 : CSLFetchNameValue( papszLocalMetadata,
2350 0 : pszZoneLine );
2351 : const char* pszEllipsoid =
2352 : CSLFetchNameValue( papszLocalMetadata,
2353 0 : pszEllipsoidLine );
2354 :
2355 : #ifdef DEBUG
2356 : CPLDebug( "HDF4Image",
2357 : "Projection %s=%s, parameters %s=%s, "
2358 : "zone %s=%s",
2359 : pszProjLine, pszProj, pszParmsLine,
2360 0 : pszParms, pszZoneLine, pszZone );
2361 : CPLDebug( "HDF4Image", "Ellipsoid %s=%s",
2362 0 : pszEllipsoidLine, pszEllipsoid );
2363 : #endif
2364 :
2365 : // Transform all mnemonical codes in the values.
2366 : int i, nParms;
2367 : // Projection is UTM by default
2368 : long iProjSys = (pszProj) ?
2369 0 : USGSMnemonicToCode(pszProj) : 1L;
2370 : long iZone =
2371 0 : (pszZone && iProjSys == 1L) ? atoi(pszZone): 0L;
2372 : char **papszEllipsoid = (pszEllipsoid) ?
2373 : CSLTokenizeString2( pszEllipsoid, ",",
2374 0 : CSLT_HONOURSTRINGS ) : NULL;
2375 :
2376 0 : long iEllipsoid = 8L; // WGS84 by default
2377 0 : if ( papszEllipsoid
2378 : && CSLCount(papszEllipsoid) > 0 )
2379 : {
2380 0 : if (EQUAL( papszEllipsoid[0], "WGS84"))
2381 0 : iEllipsoid = 8L;
2382 : }
2383 :
2384 : double adfProjParms[15];
2385 : char **papszParms = (pszParms) ?
2386 : CSLTokenizeString2( pszParms, ",",
2387 0 : CSLT_HONOURSTRINGS ) : NULL;
2388 0 : nParms = CSLCount(papszParms);
2389 0 : if (nParms >= 15)
2390 0 : nParms = 15;
2391 0 : for (i = 0; i < nParms; i++)
2392 0 : adfProjParms[i] = CPLAtof( papszParms[i] );
2393 0 : for (; i < 15; i++)
2394 0 : adfProjParms[i] = 0.0;
2395 :
2396 : // Create projection definition
2397 : oSRS.importFromUSGS( iProjSys, iZone,
2398 0 : adfProjParms, iEllipsoid );
2399 0 : oSRS.SetLinearUnits( SRS_UL_METER, 1.0 );
2400 0 : oSRS.exportToWkt( &pszGCPProjection );
2401 :
2402 0 : CSLDestroy( papszParms );
2403 0 : CPLFree( pszEllipsoidLine );
2404 0 : CPLFree( pszZoneLine );
2405 0 : CPLFree( pszParmsLine );
2406 0 : CPLFree( pszProjLine );
2407 : }
2408 :
2409 : // ASTER Level 3 (DEM)
2410 0 : else if ( eProduct == PROD_ASTER_L3 )
2411 : {
2412 : double dfCenterX, dfCenterY;
2413 : int iZone;
2414 :
2415 : ReadCoordinates( CSLFetchNameValue(
2416 : papszGlobalMetadata, "SCENECENTER" ),
2417 0 : &dfCenterY, &dfCenterX );
2418 :
2419 : // Calculate UTM zone from scene center coordinates
2420 0 : iZone = 30 + (int) ((dfCenterX + 6.0) / 6.0);
2421 :
2422 : // Create projection definition
2423 0 : if( dfCenterY > 0 )
2424 0 : oSRS.SetUTM( iZone, TRUE );
2425 : else
2426 0 : oSRS.SetUTM( - iZone, FALSE );
2427 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
2428 0 : oSRS.SetLinearUnits( SRS_UL_METER, 1.0 );
2429 0 : oSRS.exportToWkt( &pszGCPProjection );
2430 : }
2431 :
2432 : // MODIS L1B
2433 0 : else if ( eProduct == PROD_MODIS_L1B
2434 : || eProduct == PROD_MODIS_L2 )
2435 : {
2436 0 : pszGCPProjection = CPLStrdup( SRS_WKT_WGS84 );
2437 : }
2438 : }
2439 :
2440 : /* -------------------------------------------------------------------- */
2441 : /* Fill the GCPs list. */
2442 : /* -------------------------------------------------------------------- */
2443 0 : if( iGCPStepX > 0 )
2444 : {
2445 : nGCPCount = (((nXPoints-1) / iGCPStepX) + 1)
2446 0 : * (((nYPoints-1) / iGCPStepY) + 1);
2447 :
2448 0 : pasGCPList = (GDAL_GCP *) CPLCalloc( nGCPCount, sizeof( GDAL_GCP ) );
2449 0 : GDALInitGCPs( nGCPCount, pasGCPList );
2450 :
2451 0 : int iGCP = 0;
2452 0 : for ( i = 0; i < nYPoints; i += iGCPStepY )
2453 : {
2454 0 : for ( j = 0; j < nXPoints; j += iGCPStepX )
2455 : {
2456 0 : int iGeoOff = i * nXPoints + j;
2457 :
2458 0 : pasGCPList[iGCP].dfGCPX =
2459 : AnyTypeToDouble(iWrkNumType,
2460 0 : (void *)((char*)pLong+ iGeoOff*iDataSize));
2461 0 : pasGCPList[iGCP].dfGCPY =
2462 : AnyTypeToDouble(iWrkNumType,
2463 0 : (void *)((char*)pLat + iGeoOff*iDataSize));
2464 :
2465 : // GCPs in Level 1A/1B dataset are in geocentric
2466 : // coordinates. Convert them in geodetic (we
2467 : // will convert latitudes only, longitudes
2468 : // do not need to be converted, because
2469 : // they are the same).
2470 : // This calculation valid for WGS84 datum only.
2471 0 : if ( eProduct == PROD_ASTER_L1A
2472 : || eProduct == PROD_ASTER_L1B )
2473 : {
2474 0 : pasGCPList[iGCP].dfGCPY =
2475 0 : atan(tan(pasGCPList[iGCP].dfGCPY
2476 0 : *PI/180)/0.99330562)*180/PI;
2477 : }
2478 :
2479 0 : ToGeoref(&pasGCPList[iGCP].dfGCPX,
2480 0 : &pasGCPList[iGCP].dfGCPY);
2481 :
2482 0 : pasGCPList[iGCP].dfGCPZ = 0.0;
2483 :
2484 0 : if ( pLatticeX && pLatticeY )
2485 : {
2486 0 : pasGCPList[iGCP].dfGCPPixel =
2487 : AnyTypeToDouble(iLatticeType,
2488 : (void *)((char *)pLatticeX
2489 0 : + iGeoOff*iLatticeDataSize))+0.5;
2490 0 : pasGCPList[iGCP].dfGCPLine =
2491 : AnyTypeToDouble(iLatticeType,
2492 : (void *)((char *)pLatticeY
2493 0 : + iGeoOff*iLatticeDataSize))+0.5;
2494 : }
2495 0 : else if ( paiOffset && paiIncrement )
2496 : {
2497 0 : pasGCPList[iGCP].dfGCPPixel =
2498 0 : paiOffset[iPixelDim] +
2499 0 : j * paiIncrement[iPixelDim] + 0.5;
2500 0 : pasGCPList[iGCP].dfGCPLine =
2501 0 : paiOffset[iLineDim] +
2502 0 : i * paiIncrement[iLineDim] + 0.5;
2503 : }
2504 :
2505 0 : iGCP++;
2506 : }
2507 : }
2508 : }
2509 :
2510 : /* -------------------------------------------------------------------- */
2511 : /* Establish geolocation metadata, but only if there is no */
2512 : /* lattice. The lattice destroys the regularity of the grid. */
2513 : /* -------------------------------------------------------------------- */
2514 0 : if( pLatticeX == NULL
2515 : && iLatDim != -1 && iLongDim != -1
2516 : && iPixelDim != -1 && iLineDim != -1 )
2517 : {
2518 0 : CPLString osWrk;
2519 :
2520 0 : SetMetadataItem( "SRS", pszGCPProjection, "GEOLOCATION" );
2521 :
2522 : osWrk.Printf( "HDF4_EOS:EOS_SWATH_GEOL:\"%s\":%s:%s",
2523 : pszFilename, pszSubdatasetName,
2524 0 : papszGeolocations[iLongDim] );
2525 0 : SetMetadataItem( "X_DATASET", osWrk, "GEOLOCATION" );
2526 0 : SetMetadataItem( "X_BAND", "1" , "GEOLOCATION" );
2527 :
2528 : osWrk.Printf( "HDF4_EOS:EOS_SWATH_GEOL:\"%s\":%s:%s",
2529 : pszFilename, pszSubdatasetName,
2530 0 : papszGeolocations[iLatDim] );
2531 0 : SetMetadataItem( "Y_DATASET", osWrk, "GEOLOCATION" );
2532 0 : SetMetadataItem( "Y_BAND", "1" , "GEOLOCATION" );
2533 :
2534 0 : if ( paiOffset && paiIncrement )
2535 : {
2536 0 : osWrk.Printf( "%ld", (long)paiOffset[iPixelDim] );
2537 0 : SetMetadataItem( "PIXEL_OFFSET", osWrk, "GEOLOCATION" );
2538 0 : osWrk.Printf( "%ld", (long)paiIncrement[iPixelDim] );
2539 0 : SetMetadataItem( "PIXEL_STEP", osWrk, "GEOLOCATION" );
2540 :
2541 0 : osWrk.Printf( "%ld", (long)paiOffset[iLineDim] );
2542 0 : SetMetadataItem( "LINE_OFFSET", osWrk, "GEOLOCATION" );
2543 0 : osWrk.Printf( "%ld", (long)paiIncrement[iLineDim] );
2544 0 : SetMetadataItem( "LINE_STEP", osWrk, "GEOLOCATION" );
2545 0 : }
2546 : }
2547 :
2548 : /* -------------------------------------------------------------------- */
2549 : /* Cleanup */
2550 : /* -------------------------------------------------------------------- */
2551 0 : CPLFree( pLatticeX );
2552 0 : CPLFree( pLatticeY );
2553 0 : CPLFree( pLat );
2554 0 : CPLFree( pLong );
2555 :
2556 0 : CPLFree( paiOffset );
2557 0 : CPLFree( paiIncrement );
2558 0 : CPLFree( paiNumType );
2559 0 : CPLFree( paiRank );
2560 :
2561 0 : CSLDestroy( papszGeolocations );
2562 0 : CPLFree( pszGeoList );
2563 :
2564 0 : if( iGCPStepX == 0 )
2565 : {
2566 0 : CPLFree( pszGCPProjection );
2567 0 : pszGCPProjection = NULL;
2568 : }
2569 :
2570 0 : return TRUE;
2571 : }
2572 :
2573 : /************************************************************************/
2574 : /* Open() */
2575 : /************************************************************************/
2576 :
2577 13094 : GDALDataset *HDF4ImageDataset::Open( GDALOpenInfo * poOpenInfo )
2578 : {
2579 : int i;
2580 :
2581 13094 : if( !EQUALN( poOpenInfo->pszFilename, "HDF4_SDS:", 9 ) &&
2582 : !EQUALN( poOpenInfo->pszFilename, "HDF4_GR:", 8 ) &&
2583 : !EQUALN( poOpenInfo->pszFilename, "HDF4_GD:", 8 ) &&
2584 : !EQUALN( poOpenInfo->pszFilename, "HDF4_EOS:", 9 ) )
2585 12829 : return NULL;
2586 :
2587 : /* -------------------------------------------------------------------- */
2588 : /* Create a corresponding GDALDataset. */
2589 : /* -------------------------------------------------------------------- */
2590 : char **papszSubdatasetName;
2591 : HDF4ImageDataset *poDS;
2592 :
2593 265 : poDS = new HDF4ImageDataset( );
2594 265 : poDS->fp = poOpenInfo->fp;
2595 265 : poOpenInfo->fp = NULL;
2596 :
2597 265 : CPLMutexHolderD(&hHDF4Mutex);
2598 :
2599 : papszSubdatasetName = CSLTokenizeString2( poOpenInfo->pszFilename,
2600 530 : ":", CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES);
2601 265 : if ( CSLCount( papszSubdatasetName ) != 4
2602 : && CSLCount( papszSubdatasetName ) != 5
2603 : && CSLCount( papszSubdatasetName ) != 6 )
2604 : {
2605 0 : CSLDestroy( papszSubdatasetName );
2606 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
2607 0 : delete poDS;
2608 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2609 0 : return NULL;
2610 : }
2611 :
2612 : /* -------------------------------------------------------------------- */
2613 : /* Check for drive name in windows HDF4:"D:\... */
2614 : /* -------------------------------------------------------------------- */
2615 265 : if (strlen(papszSubdatasetName[2]) == 1)
2616 : {
2617 0 : char* pszFilename = (char*) CPLMalloc( 2 + strlen(papszSubdatasetName[3]) + 1);
2618 0 : sprintf(pszFilename, "%s:%s", papszSubdatasetName[2], papszSubdatasetName[3]);
2619 0 : CPLFree(papszSubdatasetName[2]);
2620 0 : CPLFree(papszSubdatasetName[3]);
2621 0 : papszSubdatasetName[2] = pszFilename;
2622 :
2623 : /* Move following arguments one rank upper */
2624 0 : papszSubdatasetName[3] = papszSubdatasetName[4];
2625 0 : if (papszSubdatasetName[4] != NULL)
2626 : {
2627 0 : papszSubdatasetName[4] = papszSubdatasetName[5];
2628 0 : papszSubdatasetName[5] = NULL;
2629 : }
2630 : }
2631 :
2632 265 : poDS->pszFilename = CPLStrdup( papszSubdatasetName[2] );
2633 :
2634 265 : if( EQUAL( papszSubdatasetName[0], "HDF4_SDS" ) )
2635 256 : poDS->iDatasetType = HDF4_SDS;
2636 9 : else if ( EQUAL( papszSubdatasetName[0], "HDF4_GR" ) )
2637 0 : poDS->iDatasetType = HDF4_GR;
2638 9 : else if ( EQUAL( papszSubdatasetName[0], "HDF4_EOS" ) )
2639 9 : poDS->iDatasetType = HDF4_EOS;
2640 : else
2641 0 : poDS->iDatasetType = HDF4_UNKNOWN;
2642 :
2643 265 : if( EQUAL( papszSubdatasetName[1], "GDAL_HDF4" ) )
2644 249 : poDS->iSubdatasetType = H4ST_GDAL;
2645 16 : else if( EQUAL( papszSubdatasetName[1], "EOS_GRID" ) )
2646 9 : poDS->iSubdatasetType = H4ST_EOS_GRID;
2647 7 : else if( EQUAL( papszSubdatasetName[1], "EOS_SWATH" ) )
2648 0 : poDS->iSubdatasetType = H4ST_EOS_SWATH;
2649 7 : else if( EQUAL( papszSubdatasetName[1], "EOS_SWATH_GEOL" ) )
2650 0 : poDS->iSubdatasetType = H4ST_EOS_SWATH_GEOL;
2651 7 : else if( EQUAL( papszSubdatasetName[1], "SEAWIFS_L3" ) )
2652 0 : poDS->iSubdatasetType= H4ST_SEAWIFS_L3;
2653 7 : else if( EQUAL( papszSubdatasetName[1], "HYPERION_L1" ) )
2654 0 : poDS->iSubdatasetType= H4ST_HYPERION_L1;
2655 : else
2656 7 : poDS->iSubdatasetType = H4ST_UNKNOWN;
2657 :
2658 : // Is our file still here?
2659 265 : if ( !Hishdf( poDS->pszFilename ) )
2660 : {
2661 0 : CSLDestroy( papszSubdatasetName );
2662 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
2663 0 : delete poDS;
2664 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2665 0 : return NULL;
2666 : }
2667 :
2668 : /* -------------------------------------------------------------------- */
2669 : /* Collect the remain (post filename) components to treat as */
2670 : /* the subdataset name. */
2671 : /* -------------------------------------------------------------------- */
2672 265 : CPLString osSubdatasetName;
2673 :
2674 265 : osSubdatasetName = papszSubdatasetName[3];
2675 265 : if( papszSubdatasetName[4] != NULL )
2676 : {
2677 9 : osSubdatasetName += ":";
2678 9 : osSubdatasetName += papszSubdatasetName[4];
2679 : }
2680 :
2681 : /* -------------------------------------------------------------------- */
2682 : /* Try opening the dataset. */
2683 : /* -------------------------------------------------------------------- */
2684 : int32 iAttribute, nValues, iAttrNumType;
2685 : double dfNoData, dfScale, dfOffset;
2686 265 : int bNoDataSet = FALSE, bHaveScale = FALSE, bHaveOffset = FALSE;
2687 265 : const char *pszUnits = NULL, *pszDescription = NULL;
2688 :
2689 : /* -------------------------------------------------------------------- */
2690 : /* Select SDS or GR to read from. */
2691 : /* -------------------------------------------------------------------- */
2692 265 : if ( poDS->iDatasetType == HDF4_EOS )
2693 : {
2694 9 : poDS->pszSubdatasetName = CPLStrdup( papszSubdatasetName[3] );
2695 9 : if (papszSubdatasetName[4] == NULL)
2696 : {
2697 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
2698 0 : delete poDS;
2699 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2700 0 : return NULL;
2701 : }
2702 9 : poDS->pszFieldName = CPLStrdup( papszSubdatasetName[4] );
2703 : }
2704 : else
2705 256 : poDS->iDataset = atoi( papszSubdatasetName[3] );
2706 265 : CSLDestroy( papszSubdatasetName );
2707 :
2708 265 : switch ( poDS->iDatasetType )
2709 : {
2710 : case HDF4_EOS:
2711 : {
2712 9 : void *pNoDataValue = NULL;
2713 :
2714 9 : switch ( poDS->iSubdatasetType )
2715 : {
2716 :
2717 : /* -------------------------------------------------------------------- */
2718 : /* HDF-EOS Swath. */
2719 : /* -------------------------------------------------------------------- */
2720 : case H4ST_EOS_SWATH:
2721 : case H4ST_EOS_SWATH_GEOL:
2722 : {
2723 : int32 hSW;
2724 0 : char *pszDimList = NULL;
2725 :
2726 0 : if( poOpenInfo->eAccess == GA_ReadOnly )
2727 0 : poDS->hHDF4 = SWopen( poDS->pszFilename, DFACC_READ );
2728 : else
2729 0 : poDS->hHDF4 = SWopen( poDS->pszFilename, DFACC_WRITE );
2730 :
2731 0 : if( poDS->hHDF4 <= 0 )
2732 : {
2733 : CPLDebug( "HDF4Image", "Can't open HDF4 file %s",
2734 0 : poDS->pszFilename );
2735 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
2736 0 : delete poDS;
2737 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2738 0 : return( NULL );
2739 : }
2740 :
2741 0 : hSW = SWattach( poDS->hHDF4, poDS->pszSubdatasetName );
2742 0 : if( hSW < 0 )
2743 : {
2744 : CPLDebug( "HDF4Image", "Can't attach to subdataset %s",
2745 0 : poDS->pszSubdatasetName );
2746 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
2747 0 : delete poDS;
2748 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2749 0 : return( NULL );
2750 : }
2751 :
2752 : /* -------------------------------------------------------------------- */
2753 : /* Decode the dimension map. */
2754 : /* -------------------------------------------------------------------- */
2755 0 : int32 nStrBufSize = 0;
2756 :
2757 0 : if ( SWnentries( hSW, HDFE_NENTDIM, &nStrBufSize ) < 0
2758 : || nStrBufSize <= 0 )
2759 : {
2760 : CPLDebug( "HDF4Image",
2761 0 : "Can't read a number of dimension maps." );
2762 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
2763 0 : delete poDS;
2764 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2765 0 : return NULL;
2766 : }
2767 :
2768 0 : pszDimList = (char *)CPLMalloc( nStrBufSize + 1 );
2769 0 : if ( SWfieldinfo( hSW, poDS->pszFieldName, &poDS->iRank,
2770 : poDS->aiDimSizes, &poDS->iNumType,
2771 : pszDimList ) < 0 )
2772 : {
2773 0 : CPLDebug( "HDF4Image", "Can't read dimension maps." );
2774 0 : CPLFree( pszDimList );
2775 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
2776 0 : delete poDS;
2777 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2778 0 : return NULL;
2779 : }
2780 0 : pszDimList[nStrBufSize] = '\0';
2781 : #ifdef DEBUG
2782 : CPLDebug( "HDF4Image",
2783 : "List of dimensions in swath \"%s\": %s",
2784 0 : poDS->pszFieldName, pszDimList );
2785 : #endif
2786 0 : poDS->GetImageDimensions( pszDimList );
2787 :
2788 : #ifdef DEBUG
2789 : CPLDebug( "HDF4Image",
2790 : "X dimension is %d, Y dimension is %d",
2791 0 : poDS->iXDim, poDS->iYDim );
2792 : #endif
2793 :
2794 : /* -------------------------------------------------------------------- */
2795 : /* Fetch metadata. */
2796 : /* -------------------------------------------------------------------- */
2797 0 : if( poDS->iSubdatasetType == H4ST_EOS_SWATH ) /* Not H4ST_EOS_SWATH_GEOL */
2798 0 : poDS->GetSwatAttrs( hSW );
2799 :
2800 : /* -------------------------------------------------------------------- */
2801 : /* Fetch NODATA value. */
2802 : /* -------------------------------------------------------------------- */
2803 : pNoDataValue =
2804 0 : CPLMalloc( poDS->GetDataTypeSize(poDS->iNumType) );
2805 0 : if ( SWgetfillvalue( hSW, poDS->pszFieldName,
2806 : pNoDataValue ) != -1 )
2807 : {
2808 : dfNoData = poDS->AnyTypeToDouble( poDS->iNumType,
2809 0 : pNoDataValue );
2810 0 : bNoDataSet = TRUE;
2811 : }
2812 : else
2813 : {
2814 : const char *pszNoData =
2815 : CSLFetchNameValue( poDS->papszLocalMetadata,
2816 0 : "_FillValue" );
2817 0 : if ( pszNoData )
2818 : {
2819 0 : dfNoData = CPLAtof( pszNoData );
2820 0 : bNoDataSet = TRUE;
2821 : }
2822 : }
2823 0 : CPLFree( pNoDataValue );
2824 :
2825 : /* -------------------------------------------------------------------- */
2826 : /* Handle Geolocation processing. */
2827 : /* -------------------------------------------------------------------- */
2828 0 : if( poDS->iSubdatasetType == H4ST_EOS_SWATH ) /* Not H4ST_SWATH_GEOL */
2829 : {
2830 : char **papszDimList =
2831 : CSLTokenizeString2( pszDimList, ",",
2832 0 : CSLT_HONOURSTRINGS );
2833 0 : if( !poDS->ProcessSwathGeolocation( hSW, papszDimList ) )
2834 : {
2835 : CPLDebug( "HDF4Image",
2836 0 : "No geolocation available for this swath." );
2837 : }
2838 0 : CSLDestroy( papszDimList );
2839 : }
2840 :
2841 : /* -------------------------------------------------------------------- */
2842 : /* Cleanup. */
2843 : /* -------------------------------------------------------------------- */
2844 0 : CPLFree( pszDimList );
2845 0 : SWdetach( hSW );
2846 : }
2847 0 : break;
2848 :
2849 : /* -------------------------------------------------------------------- */
2850 : /* HDF-EOS Grid. */
2851 : /* -------------------------------------------------------------------- */
2852 : case H4ST_EOS_GRID:
2853 : {
2854 9 : int32 hGD, iProjCode = 0, iZoneCode = 0, iSphereCode = 0;
2855 : int32 nXSize, nYSize;
2856 : char szDimList[8192];
2857 : double adfUpLeft[2], adfLowRight[2], adfProjParms[15];
2858 :
2859 9 : if( poOpenInfo->eAccess == GA_ReadOnly )
2860 9 : poDS->hHDF4 = GDopen( poDS->pszFilename, DFACC_READ );
2861 : else
2862 0 : poDS->hHDF4 = GDopen( poDS->pszFilename, DFACC_WRITE );
2863 :
2864 9 : if( poDS->hHDF4 <= 0 )
2865 : {
2866 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
2867 0 : delete poDS;
2868 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
2869 0 : return( NULL );
2870 : }
2871 :
2872 9 : hGD = GDattach( poDS->hHDF4, poDS->pszSubdatasetName );
2873 :
2874 : /* -------------------------------------------------------------------- */
2875 : /* Decode the dimension map. */
2876 : /* -------------------------------------------------------------------- */
2877 : GDfieldinfo( hGD, poDS->pszFieldName, &poDS->iRank,
2878 9 : poDS->aiDimSizes, &poDS->iNumType, szDimList );
2879 : #ifdef DEBUG
2880 : CPLDebug( "HDF4Image",
2881 : "List of dimensions in grid %s: %s",
2882 9 : poDS->pszFieldName, szDimList);
2883 : #endif
2884 9 : poDS->GetImageDimensions( szDimList );
2885 :
2886 : int32 tilecode, tilerank;
2887 9 : if( GDtileinfo( hGD , poDS->pszFieldName , &tilecode, &tilerank, NULL ) == 0 )
2888 : {
2889 9 : if ( tilecode == HDFE_TILE )
2890 : {
2891 7 : int32 *tiledims = NULL;
2892 7 : tiledims = (int32 *) CPLCalloc( tilerank , sizeof( int32 ) );
2893 7 : GDtileinfo( hGD , poDS->pszFieldName , &tilecode, &tilerank, tiledims );
2894 14 : if ( ( tilerank == 2 ) && ( poDS->iRank == tilerank ) )
2895 : {
2896 7 : poDS->nBlockPreferredXSize = tiledims[1];
2897 7 : poDS->nBlockPreferredYSize = tiledims[0];
2898 7 : poDS->bReadTile = true;
2899 : #ifdef DEBUG
2900 : CPLDebug( "HDF4_EOS:EOS_GRID:","tilerank in grid %s: %d",
2901 7 : poDS->pszFieldName , (int)tilerank );
2902 : CPLDebug( "HDF4_EOS:EOS_GRID:","tiledimens in grid %s: %d,%d",
2903 7 : poDS->pszFieldName , (int)tiledims[0] , (int)tiledims[1] );
2904 : #endif
2905 : }
2906 : #ifdef DEBUG
2907 : else
2908 : {
2909 : CPLDebug( "HDF4_EOS:EOS_GRID:","tilerank in grid %s: %d not supported",
2910 0 : poDS->pszFieldName , (int)tilerank );
2911 : }
2912 : #endif
2913 7 : CPLFree(tiledims);
2914 : }
2915 : else
2916 : {
2917 : #ifdef DEBUG
2918 : CPLDebug( "HDF4_EOS:EOS_GRID:","tilecode == HDFE_NOTILE in grid %s: %d",
2919 : poDS->pszFieldName ,
2920 2 : (int)poDS->iRank );
2921 : #endif
2922 : }
2923 : }
2924 : #ifdef DEBUG
2925 : else
2926 : {
2927 0 : CPLDebug( "HDF4_EOS:EOS_GRID:","ERROR GDtileinfo %s", poDS->pszFieldName );
2928 : }
2929 : #endif
2930 :
2931 : /* -------------------------------------------------------------------- */
2932 : /* Fetch projection information */
2933 : /* -------------------------------------------------------------------- */
2934 9 : if ( GDprojinfo( hGD, &iProjCode, &iZoneCode,
2935 : &iSphereCode, adfProjParms) >= 0 )
2936 : {
2937 : #ifdef DEBUG
2938 : CPLDebug( "HDF4Image",
2939 : "Grid projection: "
2940 : "projection code: %ld, zone code %ld, "
2941 : "sphere code %ld",
2942 : (long)iProjCode, (long)iZoneCode,
2943 9 : (long)iSphereCode );
2944 : #endif
2945 : poDS->oSRS.importFromUSGS( iProjCode, iZoneCode,
2946 : adfProjParms, iSphereCode,
2947 9 : USGS_ANGLE_RADIANS );
2948 :
2949 9 : if ( poDS->pszProjection )
2950 9 : CPLFree( poDS->pszProjection );
2951 9 : poDS->oSRS.exportToWkt( &poDS->pszProjection );
2952 : }
2953 :
2954 : /* -------------------------------------------------------------------- */
2955 : /* Fetch geotransformation matrix */
2956 : /* -------------------------------------------------------------------- */
2957 9 : if ( GDgridinfo( hGD, &nXSize, &nYSize,
2958 : adfUpLeft, adfLowRight ) >= 0 )
2959 : {
2960 : #ifdef DEBUG
2961 : CPLDebug( "HDF4Image",
2962 : "Grid geolocation: "
2963 : "top left X %f, top left Y %f, "
2964 : "low right X %f, low right Y %f, "
2965 : "cols %ld, rows %ld",
2966 : adfUpLeft[0], adfUpLeft[1],
2967 : adfLowRight[0], adfLowRight[1],
2968 9 : (long)nXSize, (long)nYSize );
2969 : #endif
2970 9 : if ( iProjCode )
2971 : {
2972 : // For projected systems coordinates are in meters
2973 : poDS->adfGeoTransform[1] =
2974 7 : (adfLowRight[0] - adfUpLeft[0]) / nXSize;
2975 : poDS->adfGeoTransform[5] =
2976 7 : (adfLowRight[1] - adfUpLeft[1]) / nYSize;
2977 7 : poDS->adfGeoTransform[0] = adfUpLeft[0];
2978 7 : poDS->adfGeoTransform[3] = adfUpLeft[1];
2979 : }
2980 : else
2981 : {
2982 : // Handle angular geographic coordinates here
2983 : poDS->adfGeoTransform[1] =
2984 : (CPLPackedDMSToDec(adfLowRight[0]) -
2985 2 : CPLPackedDMSToDec(adfUpLeft[0])) / nXSize;
2986 : poDS->adfGeoTransform[5] =
2987 : (CPLPackedDMSToDec(adfLowRight[1]) -
2988 2 : CPLPackedDMSToDec(adfUpLeft[1])) / nYSize;
2989 : poDS->adfGeoTransform[0] =
2990 2 : CPLPackedDMSToDec(adfUpLeft[0]);
2991 : poDS->adfGeoTransform[3] =
2992 2 : CPLPackedDMSToDec(adfUpLeft[1]);
2993 : }
2994 9 : poDS->adfGeoTransform[2] = 0.0;
2995 9 : poDS->adfGeoTransform[4] = 0.0;
2996 9 : poDS->bHasGeoTransform = TRUE;
2997 : }
2998 :
2999 : /* -------------------------------------------------------------------- */
3000 : /* Fetch metadata. */
3001 : /* -------------------------------------------------------------------- */
3002 9 : poDS->GetGridAttrs( hGD );
3003 :
3004 9 : GDdetach( hGD );
3005 :
3006 : /* -------------------------------------------------------------------- */
3007 : /* Fetch NODATA value. */
3008 : /* -------------------------------------------------------------------- */
3009 : pNoDataValue =
3010 9 : CPLMalloc( poDS->GetDataTypeSize(poDS->iNumType) );
3011 9 : if ( GDgetfillvalue( hGD, poDS->pszFieldName,
3012 : pNoDataValue ) != -1 )
3013 : {
3014 : dfNoData = poDS->AnyTypeToDouble( poDS->iNumType,
3015 0 : pNoDataValue );
3016 0 : bNoDataSet = TRUE;
3017 : }
3018 : else
3019 : {
3020 : const char *pszNoData =
3021 : CSLFetchNameValue( poDS->papszLocalMetadata,
3022 9 : "_FillValue" );
3023 9 : if ( pszNoData )
3024 : {
3025 7 : dfNoData = CPLAtof( pszNoData );
3026 7 : bNoDataSet = TRUE;
3027 : }
3028 : }
3029 9 : CPLFree( pNoDataValue );
3030 : }
3031 : break;
3032 :
3033 : default:
3034 : break;
3035 : }
3036 :
3037 : /* -------------------------------------------------------------------- */
3038 : /* Fetch unit type, scale, offset and description */
3039 : /* Should be similar among various HDF-EOS kinds. */
3040 : /* -------------------------------------------------------------------- */
3041 : {
3042 : const char *pszTmp =
3043 : CSLFetchNameValue( poDS->papszLocalMetadata,
3044 9 : "scale_factor" );
3045 9 : if ( pszTmp )
3046 : {
3047 7 : dfScale = CPLAtof( pszTmp );
3048 : // some producers (ie. lndcsm from LEDAPS) emit
3049 : // files with scale_factor=0 which is crazy to carry
3050 : // through.
3051 7 : if( dfScale == 0.0 )
3052 0 : dfScale = 1.0;
3053 :
3054 7 : bHaveScale = (dfScale != 0.0);
3055 : }
3056 :
3057 : pszTmp =
3058 9 : CSLFetchNameValue( poDS->papszLocalMetadata, "add_offset" );
3059 9 : if ( pszTmp )
3060 : {
3061 5 : dfOffset = CPLAtof( pszTmp );
3062 5 : bHaveOffset = TRUE;
3063 : }
3064 :
3065 : pszUnits = CSLFetchNameValue( poDS->papszLocalMetadata,
3066 9 : "units" );
3067 : pszDescription = CSLFetchNameValue( poDS->papszLocalMetadata,
3068 9 : "long_name" );
3069 : }
3070 : }
3071 9 : break;
3072 :
3073 : /* -------------------------------------------------------------------- */
3074 : /* 'Plain' HDF scientific datasets. */
3075 : /* -------------------------------------------------------------------- */
3076 : case HDF4_SDS:
3077 : {
3078 : int32 iSDS;
3079 :
3080 256 : if( poOpenInfo->eAccess == GA_ReadOnly )
3081 256 : poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_READ, 0 );
3082 : else
3083 0 : poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_WRITE, 0 );
3084 :
3085 256 : if( poDS->hHDF4 <= 0 )
3086 : {
3087 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3088 0 : delete poDS;
3089 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3090 0 : return( NULL );
3091 : }
3092 :
3093 256 : poDS->hSD = SDstart( poDS->pszFilename, DFACC_READ );
3094 256 : if ( poDS->hSD == -1 )
3095 : {
3096 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3097 0 : delete poDS;
3098 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3099 0 : return NULL;
3100 : }
3101 :
3102 256 : if ( poDS->ReadGlobalAttributes( poDS->hSD ) != CE_None )
3103 : {
3104 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3105 0 : delete poDS;
3106 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3107 0 : return NULL;
3108 : }
3109 :
3110 : int32 nDatasets, nAttrs;
3111 :
3112 256 : if ( SDfileinfo( poDS->hSD, &nDatasets, &nAttrs ) != 0 )
3113 : {
3114 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3115 0 : delete poDS;
3116 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3117 0 : return NULL;
3118 : }
3119 :
3120 256 : if (poDS->iDataset < 0 || poDS->iDataset >= nDatasets)
3121 : {
3122 : CPLError(CE_Failure, CPLE_AppDefined,
3123 : "Subdataset index should be between 0 and %ld",
3124 0 : (long int)nDatasets - 1);
3125 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3126 0 : delete poDS;
3127 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3128 0 : return NULL;
3129 : }
3130 :
3131 256 : memset( poDS->aiDimSizes, 0, sizeof(int32) * H4_MAX_VAR_DIMS );
3132 256 : iSDS = SDselect( poDS->hSD, poDS->iDataset );
3133 : SDgetinfo( iSDS, poDS->szName, &poDS->iRank, poDS->aiDimSizes,
3134 256 : &poDS->iNumType, &poDS->nAttrs);
3135 :
3136 : // We will duplicate global metadata for every subdataset
3137 : poDS->papszLocalMetadata =
3138 256 : CSLDuplicate( poDS->papszGlobalMetadata );
3139 :
3140 287 : for ( iAttribute = 0; iAttribute < poDS->nAttrs; iAttribute++ )
3141 : {
3142 : char szAttrName[H4_MAX_NC_NAME];
3143 : SDattrinfo( iSDS, iAttribute, szAttrName,
3144 31 : &iAttrNumType, &nValues );
3145 : poDS->papszLocalMetadata =
3146 : poDS->TranslateHDF4Attributes( iSDS, iAttribute,
3147 : szAttrName, iAttrNumType,
3148 : nValues,
3149 31 : poDS->papszLocalMetadata );
3150 : }
3151 256 : poDS->SetMetadata( poDS->papszLocalMetadata, "" );
3152 256 : SDendaccess( iSDS );
3153 :
3154 : #ifdef DEBUG
3155 : CPLDebug( "HDF4Image",
3156 : "aiDimSizes[0]=%ld, aiDimSizes[1]=%ld, "
3157 : "aiDimSizes[2]=%ld, aiDimSizes[3]=%ld",
3158 512 : (long)poDS->aiDimSizes[0], (long)poDS->aiDimSizes[1],
3159 768 : (long)poDS->aiDimSizes[2], (long)poDS->aiDimSizes[3] );
3160 : #endif
3161 :
3162 256 : switch( poDS->iRank )
3163 : {
3164 : case 1:
3165 0 : poDS->nBandCount = 1;
3166 0 : poDS->iXDim = 0;
3167 0 : poDS->iYDim = -1;
3168 0 : break;
3169 : case 2:
3170 110 : poDS->nBandCount = 1;
3171 110 : poDS->iXDim = 1;
3172 110 : poDS->iYDim = 0;
3173 110 : break;
3174 : case 3:
3175 : /* FIXME ? We should probably remove the following test as there are valid datasets */
3176 : /* where the height is lower than the band number : for example
3177 : http://www.iapmw.unibe.ch/research/projects/FriOWL/data/otd/LISOTD_HRAC_V2.2.hdf */
3178 : /* which is a 720x360 x 365 bands */
3179 : /* Use a HACK for now */
3180 146 : if( poDS->aiDimSizes[0] < poDS->aiDimSizes[2] &&
3181 0 : !(poDS->aiDimSizes[0] == 360 &&
3182 0 : poDS->aiDimSizes[1] == 720 &&
3183 0 : poDS->aiDimSizes[2] == 365) )
3184 : {
3185 0 : poDS->iBandDim = 0;
3186 0 : poDS->iXDim = 2;
3187 0 : poDS->iYDim = 1;
3188 : }
3189 : else
3190 : {
3191 292 : if( poDS->aiDimSizes[1] <= poDS->aiDimSizes[0] &&
3192 146 : poDS->aiDimSizes[1] <= poDS->aiDimSizes[2] )
3193 : {
3194 0 : poDS->iBandDim = 1;
3195 0 : poDS->iXDim = 2;
3196 0 : poDS->iYDim = 0;
3197 : }
3198 : else
3199 : {
3200 146 : poDS->iBandDim = 2;
3201 146 : poDS->iXDim = 1;
3202 146 : poDS->iYDim = 0;
3203 :
3204 : }
3205 : }
3206 146 : poDS->nBandCount = poDS->aiDimSizes[poDS->iBandDim];
3207 146 : break;
3208 : case 4: // FIXME
3209 0 : poDS->nBandCount = poDS->aiDimSizes[2] * poDS->aiDimSizes[3];
3210 : break;
3211 : default:
3212 : break;
3213 : }
3214 :
3215 : // We preset this because CaptureNRLGeoTransform needs it.
3216 256 : poDS->nRasterXSize = poDS->aiDimSizes[poDS->iXDim];
3217 256 : if (poDS->iYDim >= 0)
3218 256 : poDS->nRasterYSize = poDS->aiDimSizes[poDS->iYDim];
3219 : else
3220 0 : poDS->nRasterYSize = 1;
3221 :
3222 : // Special case projection info for NRL generated files.
3223 : const char *pszMapProjectionSystem =
3224 : CSLFetchNameValue(poDS->papszGlobalMetadata,
3225 256 : "mapProjectionSystem");
3226 256 : if( pszMapProjectionSystem != NULL
3227 : && EQUAL(pszMapProjectionSystem,"NRL(USGS)") )
3228 : {
3229 0 : poDS->CaptureNRLGeoTransform();
3230 : }
3231 :
3232 : // Special case for coastwatch hdf files.
3233 256 : if( CSLFetchNameValue( poDS->papszGlobalMetadata,
3234 : "gctp_sys" ) != NULL )
3235 0 : poDS->CaptureCoastwatchGCTPInfo();
3236 :
3237 : // Special case for MODIS geolocation
3238 256 : poDS->ProcessModisSDSGeolocation();
3239 :
3240 : // Special case for NASA/CCRS Landsat in HDF
3241 256 : poDS->CaptureL1GMTLInfo();
3242 : }
3243 256 : break;
3244 :
3245 : /* -------------------------------------------------------------------- */
3246 : /* 'Plain' HDF rasters. */
3247 : /* -------------------------------------------------------------------- */
3248 : case HDF4_GR:
3249 0 : if( poOpenInfo->eAccess == GA_ReadOnly )
3250 0 : poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_READ, 0 );
3251 : else
3252 0 : poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_WRITE, 0 );
3253 :
3254 0 : if( poDS->hHDF4 <= 0 )
3255 : {
3256 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3257 0 : delete poDS;
3258 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3259 0 : return( NULL );
3260 : }
3261 :
3262 0 : poDS->hGR = GRstart( poDS->hHDF4 );
3263 0 : if ( poDS->hGR == -1 )
3264 : {
3265 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3266 0 : delete poDS;
3267 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3268 0 : return NULL;
3269 : }
3270 :
3271 0 : poDS->iGR = GRselect( poDS->hGR, poDS->iDataset );
3272 0 : if ( GRgetiminfo( poDS->iGR, poDS->szName,
3273 : &poDS->iRank, &poDS->iNumType,
3274 : &poDS->iInterlaceMode, poDS->aiDimSizes,
3275 : &poDS->nAttrs ) != 0 )
3276 : {
3277 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3278 0 : delete poDS;
3279 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3280 0 : return NULL;
3281 : }
3282 :
3283 : // We will duplicate global metadata for every subdataset
3284 0 : poDS->papszLocalMetadata = CSLDuplicate( poDS->papszGlobalMetadata );
3285 :
3286 0 : for ( iAttribute = 0; iAttribute < poDS->nAttrs; iAttribute++ )
3287 : {
3288 : char szAttrName[H4_MAX_NC_NAME];
3289 : GRattrinfo( poDS->iGR, iAttribute, szAttrName,
3290 0 : &iAttrNumType, &nValues );
3291 : poDS->papszLocalMetadata =
3292 : poDS->TranslateHDF4Attributes( poDS->iGR, iAttribute,
3293 : szAttrName, iAttrNumType,
3294 : nValues,
3295 0 : poDS->papszLocalMetadata );
3296 : }
3297 0 : poDS->SetMetadata( poDS->papszLocalMetadata, "" );
3298 : // Read colour table
3299 : GDALColorEntry oEntry;
3300 :
3301 0 : poDS->iPal = GRgetlutid ( poDS->iGR, poDS->iDataset );
3302 0 : if ( poDS->iPal != -1 )
3303 : {
3304 : GRgetlutinfo( poDS->iPal, &poDS->nComps, &poDS->iPalDataType,
3305 0 : &poDS->iPalInterlaceMode, &poDS->nPalEntries );
3306 0 : GRreadlut( poDS->iPal, poDS->aiPaletteData );
3307 0 : poDS->poColorTable = new GDALColorTable();
3308 0 : for( i = 0; i < N_COLOR_ENTRIES; i++ )
3309 : {
3310 0 : oEntry.c1 = poDS->aiPaletteData[i][0];
3311 0 : oEntry.c2 = poDS->aiPaletteData[i][1];
3312 0 : oEntry.c3 = poDS->aiPaletteData[i][2];
3313 0 : oEntry.c4 = 255;
3314 :
3315 0 : poDS->poColorTable->SetColorEntry( i, &oEntry );
3316 : }
3317 : }
3318 :
3319 0 : poDS->iXDim = 0;
3320 0 : poDS->iYDim = 1;
3321 0 : poDS->nBandCount = poDS->iRank;
3322 0 : break;
3323 : default:
3324 0 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3325 0 : delete poDS;
3326 0 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3327 0 : return NULL;
3328 : }
3329 :
3330 265 : poDS->nRasterXSize = poDS->aiDimSizes[poDS->iXDim];
3331 265 : if (poDS->iYDim >= 0)
3332 265 : poDS->nRasterYSize = poDS->aiDimSizes[poDS->iYDim];
3333 : else
3334 0 : poDS->nRasterYSize = 1;
3335 :
3336 265 : if ( poDS->iSubdatasetType == H4ST_HYPERION_L1 )
3337 : {
3338 : // XXX: Hyperion SDSs has Height x Bands x Width dimensions scheme
3339 0 : if ( poDS->iRank > 2 )
3340 : {
3341 0 : poDS->nBandCount = poDS->aiDimSizes[1];
3342 0 : poDS->nRasterXSize = poDS->aiDimSizes[2];
3343 0 : poDS->nRasterYSize = poDS->aiDimSizes[0];
3344 : }
3345 : else
3346 : {
3347 0 : poDS->nBandCount = poDS->aiDimSizes[0];
3348 0 : poDS->nRasterXSize = poDS->aiDimSizes[1];
3349 0 : poDS->nRasterYSize = 1;
3350 : }
3351 : }
3352 :
3353 : /* -------------------------------------------------------------------- */
3354 : /* Create band information objects. */
3355 : /* -------------------------------------------------------------------- */
3356 574 : for( i = 1; i <= poDS->nBandCount; i++ )
3357 : {
3358 : HDF4ImageRasterBand *poBand =
3359 : new HDF4ImageRasterBand( poDS, i,
3360 309 : poDS->GetDataType( poDS->iNumType ) );
3361 309 : poDS->SetBand( i, poBand );
3362 :
3363 309 : if ( bNoDataSet )
3364 7 : poBand->SetNoDataValue( dfNoData );
3365 309 : if ( bHaveScale )
3366 : {
3367 7 : poBand->bHaveScale = TRUE;
3368 7 : poBand->dfScale = dfScale;
3369 : }
3370 309 : if ( bHaveOffset )
3371 : {
3372 5 : poBand->bHaveOffset = TRUE;
3373 5 : poBand->dfOffset = dfOffset;
3374 : }
3375 309 : if ( pszUnits )
3376 7 : poBand->osUnitType = pszUnits;
3377 309 : if ( pszDescription )
3378 9 : poBand->SetDescription( pszDescription );
3379 : }
3380 :
3381 : /* -------------------------------------------------------------------- */
3382 : /* Now we will handle particular types of HDF products. Every */
3383 : /* HDF product has its own structure. */
3384 : /* -------------------------------------------------------------------- */
3385 :
3386 : // Variables for reading georeferencing
3387 : double dfULX, dfULY, dfLRX, dfLRY;
3388 :
3389 265 : switch ( poDS->iSubdatasetType )
3390 : {
3391 : /* -------------------------------------------------------------------- */
3392 : /* HDF, created by GDAL. */
3393 : /* -------------------------------------------------------------------- */
3394 : case H4ST_GDAL:
3395 : {
3396 : const char *pszValue;
3397 :
3398 : CPLDebug( "HDF4Image",
3399 249 : "Input dataset interpreted as GDAL_HDF4" );
3400 :
3401 249 : if ( (pszValue =
3402 : CSLFetchNameValue(poDS->papszGlobalMetadata,
3403 : "Projection")) )
3404 : {
3405 105 : if ( poDS->pszProjection )
3406 105 : CPLFree( poDS->pszProjection );
3407 105 : poDS->pszProjection = CPLStrdup( pszValue );
3408 : }
3409 249 : if ( (pszValue = CSLFetchNameValue(poDS->papszGlobalMetadata,
3410 : "TransformationMatrix")) )
3411 : {
3412 249 : int i = 0;
3413 249 : char *pszString = (char *) pszValue;
3414 1992 : while ( *pszValue && i < 6 )
3415 : {
3416 1494 : poDS->adfGeoTransform[i++] = strtod(pszString, &pszString);
3417 1494 : pszString++;
3418 : }
3419 249 : poDS->bHasGeoTransform = TRUE;
3420 : }
3421 528 : for( i = 1; i <= poDS->nBands; i++ )
3422 : {
3423 279 : if ( (pszValue =
3424 : CSLFetchNameValue(poDS->papszGlobalMetadata,
3425 : CPLSPrintf("BandDesc%d", i))) )
3426 32 : poDS->GetRasterBand( i )->SetDescription( pszValue );
3427 : }
3428 528 : for( i = 1; i <= poDS->nBands; i++ )
3429 : {
3430 279 : if ( (pszValue =
3431 : CSLFetchNameValue(poDS->papszGlobalMetadata,
3432 : CPLSPrintf("NoDataValue%d", i))) )
3433 32 : poDS->GetRasterBand(i)->SetNoDataValue( CPLAtof(pszValue) );
3434 : }
3435 : }
3436 249 : break;
3437 :
3438 : /* -------------------------------------------------------------------- */
3439 : /* SeaWiFS Level 3 Standard Mapped Image Products. */
3440 : /* Organized similar to MODIS Level 3 products. */
3441 : /* -------------------------------------------------------------------- */
3442 : case H4ST_SEAWIFS_L3:
3443 : {
3444 0 : CPLDebug( "HDF4Image", "Input dataset interpreted as SEAWIFS_L3" );
3445 :
3446 : // Read band description
3447 0 : for ( i = 1; i <= poDS->nBands; i++ )
3448 : {
3449 : poDS->GetRasterBand( i )->SetDescription(
3450 0 : CSLFetchNameValue( poDS->papszGlobalMetadata, "Parameter" ) );
3451 : }
3452 :
3453 : // Read coordinate system and geotransform matrix
3454 0 : poDS->oSRS.SetWellKnownGeogCS( "WGS84" );
3455 :
3456 0 : if ( EQUAL(CSLFetchNameValue(poDS->papszGlobalMetadata,
3457 : "Map Projection"),
3458 : "Equidistant Cylindrical") )
3459 : {
3460 0 : poDS->oSRS.SetEquirectangular( 0.0, 0.0, 0.0, 0.0 );
3461 0 : poDS->oSRS.SetLinearUnits( SRS_UL_METER, 1 );
3462 0 : if ( poDS->pszProjection )
3463 0 : CPLFree( poDS->pszProjection );
3464 0 : poDS->oSRS.exportToWkt( &poDS->pszProjection );
3465 : }
3466 :
3467 : dfULX = CPLAtof( CSLFetchNameValue(poDS->papszGlobalMetadata,
3468 0 : "Westernmost Longitude") );
3469 : dfULY = CPLAtof( CSLFetchNameValue(poDS->papszGlobalMetadata,
3470 0 : "Northernmost Latitude") );
3471 : dfLRX = CPLAtof( CSLFetchNameValue(poDS->papszGlobalMetadata,
3472 0 : "Easternmost Longitude") );
3473 : dfLRY = CPLAtof( CSLFetchNameValue(poDS->papszGlobalMetadata,
3474 0 : "Southernmost Latitude") );
3475 0 : poDS->ToGeoref( &dfULX, &dfULY );
3476 0 : poDS->ToGeoref( &dfLRX, &dfLRY );
3477 0 : poDS->adfGeoTransform[0] = dfULX;
3478 0 : poDS->adfGeoTransform[3] = dfULY;
3479 0 : poDS->adfGeoTransform[1] = (dfLRX - dfULX) / poDS->nRasterXSize;
3480 0 : poDS->adfGeoTransform[5] = (dfULY - dfLRY) / poDS->nRasterYSize;
3481 0 : if ( dfULY > 0) // Northern hemisphere
3482 0 : poDS->adfGeoTransform[5] = - poDS->adfGeoTransform[5];
3483 0 : poDS->adfGeoTransform[2] = 0.0;
3484 0 : poDS->adfGeoTransform[4] = 0.0;
3485 0 : poDS->bHasGeoTransform = TRUE;
3486 : }
3487 0 : break;
3488 :
3489 :
3490 : /* -------------------------------------------------------------------- */
3491 : /* Generic SDS */
3492 : /* -------------------------------------------------------------------- */
3493 : case H4ST_UNKNOWN:
3494 : {
3495 :
3496 : // This is a coastwatch convention.
3497 7 : if( CSLFetchNameValue( poDS->papszLocalMetadata, "missing_value" ) )
3498 : {
3499 : int i;
3500 0 : for( i = 1; i <= poDS->nBands; i++ )
3501 : {
3502 : poDS->GetRasterBand(i)->SetNoDataValue(
3503 : CPLAtof( CSLFetchNameValue(poDS->papszLocalMetadata,
3504 0 : "missing_value") ) );
3505 : }
3506 : }
3507 :
3508 : // Coastwatch offset and scale.
3509 7 : if( CSLFetchNameValue( poDS->papszLocalMetadata, "scale_factor" )
3510 : && CSLFetchNameValue( poDS->papszLocalMetadata, "add_offset" ) )
3511 : {
3512 0 : for( i = 1; i <= poDS->nBands; i++ )
3513 : {
3514 : HDF4ImageRasterBand *poBand =
3515 0 : (HDF4ImageRasterBand *) poDS->GetRasterBand(i);
3516 :
3517 0 : poBand->bHaveScale = poBand->bHaveOffset = TRUE;
3518 : poBand->dfScale =
3519 : CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3520 0 : "scale_factor" ) );
3521 : // See #4891 regarding offset interpretation.
3522 : //poBand->dfOffset = -1 * poBand->dfScale *
3523 : // CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3524 : // "add_offset" ) );
3525 : poBand->dfOffset =
3526 : CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3527 0 : "add_offset" ) );
3528 : }
3529 : }
3530 :
3531 : // this is a modis level3 convention (data from ACT)
3532 : // Eg data/hdf/act/modis/MODAM2004280160000.L3_NOAA_GMX
3533 :
3534 7 : if( CSLFetchNameValue( poDS->papszLocalMetadata,
3535 : "scalingSlope" )
3536 : && CSLFetchNameValue( poDS->papszLocalMetadata,
3537 : "scalingIntercept" ) )
3538 : {
3539 : int i;
3540 0 : CPLString osUnits;
3541 :
3542 0 : if( CSLFetchNameValue( poDS->papszLocalMetadata,
3543 : "productUnits" ) )
3544 : {
3545 : osUnits = CSLFetchNameValue( poDS->papszLocalMetadata,
3546 0 : "productUnits" );
3547 : }
3548 :
3549 0 : for( i = 1; i <= poDS->nBands; i++ )
3550 : {
3551 : HDF4ImageRasterBand *poBand =
3552 0 : (HDF4ImageRasterBand *) poDS->GetRasterBand(i);
3553 :
3554 0 : poBand->bHaveScale = poBand->bHaveOffset = TRUE;
3555 : poBand->dfScale =
3556 : CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3557 0 : "scalingSlope" ) );
3558 : poBand->dfOffset =
3559 : CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3560 0 : "scalingIntercept" ) );
3561 :
3562 0 : poBand->osUnitType = osUnits;
3563 0 : }
3564 : }
3565 : }
3566 7 : break;
3567 :
3568 : /* -------------------------------------------------------------------- */
3569 : /* Hyperion Level 1. */
3570 : /* -------------------------------------------------------------------- */
3571 : case H4ST_HYPERION_L1:
3572 : {
3573 0 : CPLDebug( "HDF4Image", "Input dataset interpreted as HYPERION_L1" );
3574 : }
3575 : break;
3576 :
3577 : default:
3578 : break;
3579 : }
3580 :
3581 : /* -------------------------------------------------------------------- */
3582 : /* Setup PAM info for this subdataset */
3583 : /* -------------------------------------------------------------------- */
3584 265 : poDS->SetPhysicalFilename( poDS->pszFilename );
3585 265 : poDS->SetSubdatasetName( osSubdatasetName );
3586 :
3587 265 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3588 265 : poDS->TryLoadXML();
3589 :
3590 265 : poDS->oOvManager.Initialize( poDS, ":::VIRTUAL:::" );
3591 265 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3592 :
3593 265 : return( poDS );
3594 : }
3595 :
3596 : /************************************************************************/
3597 : /* Create() */
3598 : /************************************************************************/
3599 :
3600 162 : GDALDataset *HDF4ImageDataset::Create( const char * pszFilename,
3601 : int nXSize, int nYSize, int nBands,
3602 : GDALDataType eType,
3603 : char **papszOptions )
3604 :
3605 : {
3606 : /* -------------------------------------------------------------------- */
3607 : /* Create the dataset. */
3608 : /* -------------------------------------------------------------------- */
3609 : HDF4ImageDataset *poDS;
3610 : const char *pszSDSName;
3611 : int iBand;
3612 162 : int32 iSDS = -1;
3613 : int32 aiDimSizes[H4_MAX_VAR_DIMS];
3614 :
3615 162 : if( nBands == 0 )
3616 : {
3617 : CPLError( CE_Failure, CPLE_NotSupported,
3618 1 : "Unable to export files with zero bands." );
3619 1 : return NULL;
3620 : }
3621 :
3622 161 : poDS = new HDF4ImageDataset();
3623 :
3624 161 : CPLMutexHolderD(&hHDF4Mutex);
3625 :
3626 : /* -------------------------------------------------------------------- */
3627 : /* Choose rank for the created dataset. */
3628 : /* -------------------------------------------------------------------- */
3629 161 : poDS->iRank = 3;
3630 322 : if ( CSLFetchNameValue( papszOptions, "RANK" ) != NULL &&
3631 : EQUAL( CSLFetchNameValue( papszOptions, "RANK" ), "2" ) )
3632 48 : poDS->iRank = 2;
3633 :
3634 161 : poDS->hSD = SDstart( pszFilename, DFACC_CREATE );
3635 161 : if ( poDS->hSD == -1 )
3636 : {
3637 : CPLError( CE_Failure, CPLE_OpenFailed,
3638 17 : "Can't create HDF4 file %s", pszFilename );
3639 17 : CPLReleaseMutex(hHDF4Mutex); // Release mutex otherwise we'll deadlock with GDALDataset own mutex
3640 17 : delete poDS;
3641 17 : CPLAcquireMutex(hHDF4Mutex, 1000.0);
3642 17 : return NULL;
3643 : }
3644 144 : poDS->iXDim = 1;
3645 144 : poDS->iYDim = 0;
3646 144 : poDS->iBandDim = 2;
3647 144 : aiDimSizes[poDS->iXDim] = nXSize;
3648 144 : aiDimSizes[poDS->iYDim] = nYSize;
3649 144 : aiDimSizes[poDS->iBandDim] = nBands;
3650 :
3651 144 : if ( poDS->iRank == 2 )
3652 : {
3653 96 : for ( iBand = 0; iBand < nBands; iBand++ )
3654 : {
3655 48 : pszSDSName = CPLSPrintf( "Band%d", iBand );
3656 48 : switch ( eType )
3657 : {
3658 : case GDT_Float64:
3659 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_FLOAT64,
3660 6 : poDS->iRank, aiDimSizes );
3661 6 : break;
3662 : case GDT_Float32:
3663 : iSDS = SDcreate( poDS-> hSD, pszSDSName, DFNT_FLOAT32,
3664 6 : poDS->iRank, aiDimSizes );
3665 6 : break;
3666 : case GDT_UInt32:
3667 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT32,
3668 6 : poDS->iRank, aiDimSizes );
3669 6 : break;
3670 : case GDT_UInt16:
3671 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT16,
3672 6 : poDS->iRank, aiDimSizes );
3673 6 : break;
3674 : case GDT_Int32:
3675 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT32,
3676 6 : poDS->iRank, aiDimSizes );
3677 6 : break;
3678 : case GDT_Int16:
3679 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT16,
3680 6 : poDS->iRank, aiDimSizes );
3681 6 : break;
3682 : case GDT_Byte:
3683 : default:
3684 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT8,
3685 12 : poDS->iRank, aiDimSizes );
3686 : break;
3687 : }
3688 48 : SDendaccess( iSDS );
3689 : }
3690 : }
3691 96 : else if ( poDS->iRank == 3 )
3692 : {
3693 96 : pszSDSName = "3-dimensional Scientific Dataset";
3694 96 : poDS->iDataset = 0;
3695 96 : switch ( eType )
3696 : {
3697 : case GDT_Float64:
3698 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_FLOAT64,
3699 10 : poDS->iRank, aiDimSizes );
3700 10 : break;
3701 : case GDT_Float32:
3702 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_FLOAT32,
3703 10 : poDS->iRank, aiDimSizes );
3704 10 : break;
3705 : case GDT_UInt32:
3706 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT32,
3707 10 : poDS->iRank, aiDimSizes );
3708 10 : break;
3709 : case GDT_UInt16:
3710 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT16,
3711 10 : poDS->iRank, aiDimSizes );
3712 10 : break;
3713 : case GDT_Int32:
3714 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT32,
3715 10 : poDS->iRank, aiDimSizes );
3716 10 : break;
3717 : case GDT_Int16:
3718 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT16,
3719 10 : poDS->iRank, aiDimSizes );
3720 10 : break;
3721 : case GDT_Byte:
3722 : default:
3723 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT8,
3724 36 : poDS->iRank, aiDimSizes );
3725 : break;
3726 : }
3727 : }
3728 : else // Should never happen
3729 0 : return NULL;
3730 :
3731 144 : if ( iSDS < 0 )
3732 : {
3733 : CPLError( CE_Failure, CPLE_AppDefined,
3734 : "Can't create SDS with rank %ld for file %s",
3735 0 : (long)poDS->iRank, pszFilename );
3736 0 : return NULL;
3737 : }
3738 :
3739 144 : poDS->nRasterXSize = nXSize;
3740 144 : poDS->nRasterYSize = nYSize;
3741 144 : poDS->eAccess = GA_Update;
3742 144 : poDS->iDatasetType = HDF4_SDS;
3743 144 : poDS->iSubdatasetType = H4ST_GDAL;
3744 144 : poDS->nBands = nBands;
3745 :
3746 : /* -------------------------------------------------------------------- */
3747 : /* Create band information objects. */
3748 : /* -------------------------------------------------------------------- */
3749 328 : for( iBand = 1; iBand <= nBands; iBand++ )
3750 184 : poDS->SetBand( iBand, new HDF4ImageRasterBand( poDS, iBand, eType ) );
3751 :
3752 : SDsetattr( poDS->hSD, "Signature", DFNT_CHAR8, strlen(pszGDALSignature) + 1,
3753 144 : pszGDALSignature );
3754 :
3755 144 : return (GDALDataset *) poDS;
3756 : }
3757 :
3758 : /************************************************************************/
3759 : /* GDALRegister_HDF4Image() */
3760 : /************************************************************************/
3761 :
3762 610 : void GDALRegister_HDF4Image()
3763 :
3764 : {
3765 : GDALDriver *poDriver;
3766 :
3767 610 : if( GDALGetDriverByName( "HDF4Image" ) == NULL )
3768 : {
3769 588 : poDriver = new GDALDriver();
3770 :
3771 588 : poDriver->SetDescription( "HDF4Image" );
3772 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
3773 588 : "HDF4 Dataset" );
3774 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
3775 588 : "frmt_hdf4.html" );
3776 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
3777 588 : "Byte Int16 UInt16 Int32 UInt32 Float32 Float64" );
3778 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
3779 : "<CreationOptionList>"
3780 : " <Option name='RANK' type='int' description='Rank of output SDS'/>"
3781 588 : "</CreationOptionList>" );
3782 :
3783 588 : poDriver->pfnOpen = HDF4ImageDataset::Open;
3784 588 : poDriver->pfnCreate = HDF4ImageDataset::Create;
3785 :
3786 588 : GetGDALDriverManager()->RegisterDriver( poDriver );
3787 : }
3788 610 : }
|
