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