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