1 : /******************************************************************************
2 : * $Id: hdf4imagedataset.cpp 19494 2010-04-22 04:19:35Z warmerdam $
3 : *
4 : * Project: Hierarchical Data Format Release 4 (HDF4)
5 : * Purpose: Read subdatasets of HDF4 file.
6 : * This driver initially based on code supplied by Markus Neteler
7 : * Author: Andrey Kiselev, dron@ak4719.spb.edu
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2002, Andrey Kiselev <dron@ak4719.spb.edu>
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ****************************************************************************/
30 :
31 : #include <string.h>
32 : #include <math.h>
33 :
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 19494 2010-04-22 04:19:35Z warmerdam $");
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 : };
76 :
77 : /************************************************************************/
78 : /* ==================================================================== */
79 : /* HDF4ImageDataset */
80 : /* ==================================================================== */
81 : /************************************************************************/
82 :
83 : class HDF4ImageDataset : public HDF4Dataset
84 : {
85 : friend class HDF4ImageRasterBand;
86 :
87 : char *pszFilename;
88 : int32 hHDF4, iGR, iPal, iDataset;
89 : int32 iRank, iNumType, nAttrs,
90 : iInterlaceMode, iPalInterlaceMode, iPalDataType;
91 : int32 nComps, nPalEntries;
92 : int32 aiDimSizes[H4_MAX_VAR_DIMS];
93 : int iXDim, iYDim, iBandDim, i4Dim;
94 : int nBandCount;
95 : char **papszLocalMetadata;
96 : #define N_COLOR_ENTRIES 256
97 : uint8 aiPaletteData[N_COLOR_ENTRIES][3]; // XXX: Static array for now
98 : char szName[HDF4_SDS_MAXNAMELEN];
99 : char *pszSubdatasetName;
100 : char *pszFieldName;
101 :
102 : GDALColorTable *poColorTable;
103 :
104 : OGRSpatialReference oSRS;
105 : int bHasGeoTransform;
106 : double adfGeoTransform[6];
107 : char *pszProjection;
108 : char *pszGCPProjection;
109 : GDAL_GCP *pasGCPList;
110 : int nGCPCount;
111 :
112 : HDF4DatasetType iDatasetType;
113 :
114 : int32 iSDS;
115 :
116 : int nBlockPreferredXSize;
117 : int nBlockPreferredYSize;
118 : bool bReadTile;
119 :
120 : void ToGeoref( double *, double * );
121 : void GetImageDimensions( char * );
122 : void GetSwatAttrs( int32 );
123 : void GetGridAttrs( int32 hGD );
124 : void CaptureNRLGeoTransform(void);
125 : void CaptureL1GMTLInfo(void);
126 : void CaptureCoastwatchGCTPInfo(void);
127 : void ProcessModisSDSGeolocation(void);
128 : int ProcessSwathGeolocation( int32, char ** );
129 :
130 : static long USGSMnemonicToCode( const char* );
131 : static void ReadCoordinates( const char*, double*, double* );
132 :
133 : public:
134 : HDF4ImageDataset();
135 : ~HDF4ImageDataset();
136 :
137 : static GDALDataset *Open( GDALOpenInfo * );
138 : static GDALDataset *Create( const char * pszFilename,
139 : int nXSize, int nYSize, int nBands,
140 : GDALDataType eType, char ** papszParmList );
141 : virtual void FlushCache( void );
142 : CPLErr GetGeoTransform( double * padfTransform );
143 : virtual CPLErr SetGeoTransform( double * );
144 : const char *GetProjectionRef();
145 : virtual CPLErr SetProjection( const char * );
146 : virtual int GetGCPCount();
147 : virtual const char *GetGCPProjection();
148 : virtual const GDAL_GCP *GetGCPs();
149 : };
150 :
151 : /************************************************************************/
152 : /* ==================================================================== */
153 : /* HDF4ImageRasterBand */
154 : /* ==================================================================== */
155 : /************************************************************************/
156 :
157 : class HDF4ImageRasterBand : public GDALPamRasterBand
158 539 : {
159 : friend class HDF4ImageDataset;
160 :
161 : int bNoDataSet;
162 : double dfNoDataValue;
163 :
164 : int bHaveScaleAndOffset;
165 : double dfScale;
166 : double dfOffset;
167 :
168 : CPLString osUnitType;
169 :
170 : public:
171 :
172 : HDF4ImageRasterBand( HDF4ImageDataset *, int, GDALDataType );
173 :
174 : virtual CPLErr IReadBlock( int, int, void * );
175 : virtual CPLErr IWriteBlock( int, int, void * );
176 : virtual GDALColorInterp GetColorInterpretation();
177 : virtual GDALColorTable *GetColorTable();
178 : virtual double GetNoDataValue( int * );
179 : virtual CPLErr SetNoDataValue( double );
180 : virtual double GetOffset( int *pbSuccess );
181 : virtual double GetScale( int *pbSuccess );
182 : virtual const char *GetUnitType();
183 : };
184 :
185 : /************************************************************************/
186 : /* HDF4ImageRasterBand() */
187 : /************************************************************************/
188 :
189 : HDF4ImageRasterBand::HDF4ImageRasterBand( HDF4ImageDataset *poDS, int nBand,
190 539 : GDALDataType eType )
191 :
192 : {
193 539 : this->poDS = poDS;
194 539 : this->nBand = nBand;
195 539 : eDataType = eType;
196 539 : bNoDataSet = FALSE;
197 539 : dfNoDataValue = -9999.0;
198 :
199 539 : bHaveScaleAndOffset = FALSE;
200 539 : dfScale = 1.0;
201 539 : dfOffset = 0.0;
202 :
203 539 : nBlockXSize = poDS->GetRasterXSize();
204 :
205 : // Aim for a block of about 1000000 pixels. Chunking up substantially
206 : // improves performance in some situations. For now we only chunk up for
207 : // SDS and EOS based datasets since other variations haven't been tested. #2208
208 1078 : if( poDS->iDatasetType == HDF4_SDS ||
209 : poDS->iDatasetType == HDF4_EOS)
210 : {
211 : int nChunkSize =
212 539 : atoi( CPLGetConfigOption("HDF4_BLOCK_PIXELS", "1000000") );
213 :
214 539 : nBlockYSize = nChunkSize / poDS->GetRasterXSize();
215 539 : nBlockYSize = MAX(1,MIN(nBlockYSize,poDS->GetRasterYSize()));
216 : }
217 : else
218 : {
219 0 : nBlockYSize = 1;
220 : }
221 :
222 : /* HDF4_EOS:EOS_GRID case. We ensure that the block size matches */
223 : /* the raster width, as the IReadBlock() code can only handle multiple */
224 : /* blocks per row */
225 539 : if ( poDS->nBlockPreferredXSize == nBlockXSize &&
226 : poDS->nBlockPreferredYSize > 0 )
227 : {
228 6 : if (poDS->nBlockPreferredYSize == 1)
229 : {
230 : /* Avoid defaulting to tile reading when the preferred height is 1 */
231 : /* as it leads to very poor performance with : */
232 : /* ftp://e4ftl01u.ecs.nasa.gov/MODIS_Composites/MOLT/MOD13Q1.005/2006.06.10/MOD13Q1.A2006161.h21v13.005.2008234103220.hd */
233 2 : poDS->bReadTile = FALSE;
234 : }
235 : else
236 : {
237 4 : nBlockYSize = poDS->nBlockPreferredYSize;
238 : }
239 : }
240 539 : }
241 :
242 : /************************************************************************/
243 : /* IReadBlock() */
244 : /************************************************************************/
245 :
246 : CPLErr HDF4ImageRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
247 195 : void * pImage )
248 : {
249 195 : HDF4ImageDataset *poGDS = (HDF4ImageDataset *) poDS;
250 : int32 aiStart[H4_MAX_NC_DIMS], aiEdges[H4_MAX_NC_DIMS];
251 195 : CPLErr eErr = CE_None;
252 :
253 195 : if( poGDS->eAccess == GA_Update )
254 : {
255 : memset( pImage, 0,
256 0 : nBlockXSize * nBlockYSize * GDALGetDataTypeSize(eDataType) / 8 );
257 0 : return CE_None;
258 : }
259 :
260 : /* -------------------------------------------------------------------- */
261 : /* Work out some block oriented details. */
262 : /* -------------------------------------------------------------------- */
263 195 : int nYOff = nBlockYOff * nBlockYSize;
264 195 : int nYSize = MIN(nYOff + nBlockYSize, poDS->GetRasterYSize()) - nYOff;
265 195 : CPLAssert( nBlockXOff == 0 );
266 :
267 : /* -------------------------------------------------------------------- */
268 : /* HDF files with external data files, such as some landsat */
269 : /* products (eg. data/hdf/L1G) need to be told what directory */
270 : /* to look in to find the external files. Normally this is the */
271 : /* directory holding the hdf file. */
272 : /* -------------------------------------------------------------------- */
273 195 : HXsetdir(CPLGetPath(poGDS->pszFilename));
274 :
275 : /* -------------------------------------------------------------------- */
276 : /* Handle different configurations. */
277 : /* -------------------------------------------------------------------- */
278 195 : switch ( poGDS->iDatasetType )
279 : {
280 : case HDF4_SDS:
281 : {
282 : /* We avoid doing SDselect() / SDendaccess() for each block access */
283 : /* as this is very slow when zlib compression is used */
284 :
285 94 : if (poGDS->iSDS == FAIL)
286 44 : poGDS->iSDS = SDselect( poGDS->hSD, poGDS->iDataset );
287 :
288 : /* HDF rank:
289 : A rank 2 dataset is an image read in scan-line order (2D).
290 : A rank 3 dataset is a series of images which are read in
291 : an image at a time to form a volume.
292 : A rank 4 dataset may be thought of as a series of volumes.
293 :
294 : The "aiStart" array specifies the multi-dimensional index of the
295 : starting corner of the hyperslab to read. The values are zero
296 : based.
297 :
298 : The "edge" array specifies the number of values to read along
299 : each dimension of the hyperslab.
300 :
301 : The "iStride" array allows for sub-sampling along each
302 : dimension. If a iStride value is specified for a dimension,
303 : that many values will be skipped over when reading along that
304 : dimension. Specifying iStride = NULL in the C interface or
305 : iStride = 1 in either interface specifies contiguous reading
306 : of data. If the iStride values are set to 0, SDreaddata
307 : returns FAIL (or -1). No matter what iStride value is
308 : provided, data is always placed contiguously in buffer.
309 : */
310 94 : switch ( poGDS->iRank )
311 : {
312 : case 4: // 4Dim: volume-time
313 : // FIXME: needs sample file. Does not work currently.
314 0 : aiStart[3] = 0/* range: 0--aiDimSizes[3]-1 */;
315 0 : aiEdges[3] = 1;
316 0 : aiStart[2] = 0/* range: 0--aiDimSizes[2]-1 */;
317 0 : aiEdges[2] = 1;
318 0 : aiStart[1] = nYOff; aiEdges[1] = nYSize;
319 0 : aiStart[0] = nBlockXOff; aiEdges[0] = nBlockXSize;
320 0 : break;
321 : case 3: // 3Dim: volume
322 40 : aiStart[poGDS->iBandDim] = nBand - 1;
323 40 : aiEdges[poGDS->iBandDim] = 1;
324 :
325 40 : aiStart[poGDS->iYDim] = nYOff;
326 40 : aiEdges[poGDS->iYDim] = nYSize;
327 :
328 40 : aiStart[poGDS->iXDim] = nBlockXOff;
329 40 : aiEdges[poGDS->iXDim] = nBlockXSize;
330 40 : break;
331 : case 2: // 2Dim: rows/cols
332 54 : aiStart[poGDS->iYDim] = nYOff;
333 54 : aiEdges[poGDS->iYDim] = nYSize;
334 :
335 54 : aiStart[poGDS->iXDim] = nBlockXOff;
336 54 : aiEdges[poGDS->iXDim] = nBlockXSize;
337 54 : break;
338 : case 1: //1Dim:
339 0 : aiStart[poGDS->iXDim] = nBlockXOff;
340 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
341 : break;
342 : }
343 :
344 : // Read HDF SDS array
345 94 : if( SDreaddata( poGDS->iSDS, aiStart, NULL, aiEdges, pImage ) < 0 )
346 : {
347 : CPLError( CE_Failure, CPLE_AppDefined,
348 0 : "SDreaddata() failed for block." );
349 0 : eErr = CE_Failure;
350 : }
351 :
352 : //SDendaccess( iSDS );
353 : }
354 94 : break;
355 :
356 : case HDF4_GR:
357 : {
358 : int nDataTypeSize =
359 0 : GDALGetDataTypeSize(poGDS->GetDataType(poGDS->iNumType)) / 8;
360 : GByte *pbBuffer = (GByte *)
361 0 : CPLMalloc(nBlockXSize*nBlockYSize*poGDS->iRank*nBlockYSize);
362 : int i, j;
363 :
364 0 : aiStart[poGDS->iYDim] = nYOff;
365 0 : aiEdges[poGDS->iYDim] = nYSize;
366 :
367 0 : aiStart[poGDS->iXDim] = nBlockXOff;
368 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
369 :
370 0 : if( GRreadimage(poGDS->iGR, aiStart, NULL, aiEdges, pbBuffer) < 0 )
371 : {
372 : CPLError( CE_Failure, CPLE_AppDefined,
373 0 : "GRreaddata() failed for block." );
374 0 : eErr = CE_Failure;
375 : }
376 : else
377 : {
378 0 : for ( i = 0, j = (nBand - 1) * nDataTypeSize;
379 : i < nBlockXSize * nDataTypeSize;
380 : i += nDataTypeSize, j += poGDS->nBands * nDataTypeSize )
381 0 : memcpy( (GByte *)pImage + i, pbBuffer + j, nDataTypeSize );
382 : }
383 :
384 0 : CPLFree( pbBuffer );
385 : }
386 0 : break;
387 :
388 : case HDF4_EOS:
389 : {
390 101 : switch ( poGDS->iSubdatasetType )
391 : {
392 : case H4ST_EOS_GRID:
393 : {
394 : int32 hGD;
395 :
396 : hGD = GDattach( poGDS->hHDF4,
397 101 : poGDS->pszSubdatasetName );
398 101 : switch ( poGDS->iRank )
399 : {
400 : case 4: // 4Dim: volume
401 : aiStart[poGDS->i4Dim] =
402 : (nBand - 1)
403 0 : / poGDS->aiDimSizes[poGDS->iBandDim];
404 0 : aiEdges[poGDS->i4Dim] = 1;
405 :
406 : aiStart[poGDS->iBandDim] =
407 : (nBand - 1)
408 0 : % poGDS->aiDimSizes[poGDS->iBandDim];
409 0 : aiEdges[poGDS->iBandDim] = 1;
410 :
411 0 : aiStart[poGDS->iYDim] = nYOff;
412 0 : aiEdges[poGDS->iYDim] = nYSize;
413 :
414 0 : aiStart[poGDS->iXDim] = nBlockXOff;
415 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
416 0 : break;
417 : case 3: // 3Dim: volume
418 0 : aiStart[poGDS->iBandDim] = nBand - 1;
419 0 : aiEdges[poGDS->iBandDim] = 1;
420 :
421 0 : aiStart[poGDS->iYDim] = nYOff;
422 0 : aiEdges[poGDS->iYDim] = nYSize;
423 :
424 0 : aiStart[poGDS->iXDim] = nBlockXOff;
425 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
426 0 : break;
427 : case 2: // 2Dim: rows/cols
428 101 : aiStart[poGDS->iYDim] = nYOff;
429 101 : aiEdges[poGDS->iYDim] = nYSize;
430 :
431 101 : aiStart[poGDS->iXDim] = nBlockXOff;
432 101 : aiEdges[poGDS->iXDim] = nBlockXSize;
433 : break;
434 : }
435 :
436 : /* Ensure that we don't overlap the bottom or right edges */
437 : /* of the dataset in order to use the GDreadtile() API */
438 177 : if ( poGDS->bReadTile &&
439 : (nBlockXOff + 1) * nBlockXSize <= nRasterXSize &&
440 : (nBlockYOff + 1) * nBlockYSize <= nRasterYSize )
441 : {
442 76 : int32 tilecoords[] = { nBlockYOff , nBlockXOff };
443 76 : if( GDreadtile( hGD, poGDS->pszFieldName , tilecoords , pImage ) != 0 )
444 : {
445 0 : CPLError( CE_Failure, CPLE_AppDefined, "GDreadtile() failed for block." );
446 0 : eErr = CE_Failure;
447 : }
448 : }
449 25 : else if( GDreadfield( hGD, poGDS->pszFieldName,
450 : aiStart, NULL, aiEdges, pImage ) < 0 )
451 : {
452 : CPLError( CE_Failure, CPLE_AppDefined,
453 0 : "GDreadfield() failed for block." );
454 0 : eErr = CE_Failure;
455 : }
456 101 : GDdetach( hGD );
457 : }
458 101 : break;
459 :
460 : case H4ST_EOS_SWATH:
461 : case H4ST_EOS_SWATH_GEOL:
462 : {
463 : int32 hSW;
464 :
465 : hSW = SWattach( poGDS->hHDF4,
466 0 : poGDS->pszSubdatasetName );
467 0 : switch ( poGDS->iRank )
468 : {
469 : case 3: // 3Dim: volume
470 0 : aiStart[poGDS->iBandDim] = nBand - 1;
471 0 : aiEdges[poGDS->iBandDim] = 1;
472 :
473 0 : aiStart[poGDS->iYDim] = nYOff;
474 0 : aiEdges[poGDS->iYDim] = nYSize;
475 :
476 0 : aiStart[poGDS->iXDim] = nBlockXOff;
477 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
478 0 : break;
479 : case 2: // 2Dim: rows/cols
480 0 : aiStart[poGDS->iYDim] = nYOff;
481 0 : aiEdges[poGDS->iYDim] = nYSize;
482 :
483 0 : aiStart[poGDS->iXDim] = nBlockXOff;
484 0 : aiEdges[poGDS->iXDim] = nBlockXSize;
485 : break;
486 : }
487 0 : if( SWreadfield( hSW, poGDS->pszFieldName,
488 : aiStart, NULL, aiEdges, pImage ) < 0 )
489 : {
490 : CPLError( CE_Failure, CPLE_AppDefined,
491 0 : "SWreadfield() failed for block." );
492 0 : eErr = CE_Failure;
493 : }
494 0 : SWdetach( hSW );
495 : }
496 : break;
497 : default:
498 : break;
499 : }
500 : }
501 101 : break;
502 :
503 : default:
504 0 : eErr = CE_Failure;
505 : break;
506 : }
507 :
508 195 : return eErr;
509 : }
510 :
511 : /************************************************************************/
512 : /* IWriteBlock() */
513 : /************************************************************************/
514 :
515 : CPLErr HDF4ImageRasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff,
516 65 : void * pImage )
517 : {
518 65 : HDF4ImageDataset *poGDS = (HDF4ImageDataset *)poDS;
519 : int32 aiStart[H4_MAX_NC_DIMS], aiEdges[H4_MAX_NC_DIMS];
520 65 : CPLErr eErr = CE_None;
521 :
522 65 : CPLAssert( poGDS != NULL
523 : && nBlockXOff == 0
524 : && nBlockYOff >= 0
525 : && pImage != NULL );
526 :
527 : /* -------------------------------------------------------------------- */
528 : /* Work out some block oriented details. */
529 : /* -------------------------------------------------------------------- */
530 65 : int nYOff = nBlockYOff * nBlockYSize;
531 65 : int nYSize = MIN(nYOff + nBlockYSize, poDS->GetRasterYSize()) - nYOff;
532 65 : CPLAssert( nBlockXOff == 0 );
533 :
534 : /* -------------------------------------------------------------------- */
535 : /* Process based on rank. */
536 : /* -------------------------------------------------------------------- */
537 65 : switch ( poGDS->iRank )
538 : {
539 : case 3:
540 : {
541 57 : int32 iSDS = SDselect( poGDS->hSD, poGDS->iDataset );
542 :
543 57 : aiStart[poGDS->iBandDim] = nBand - 1;
544 57 : aiEdges[poGDS->iBandDim] = 1;
545 :
546 57 : aiStart[poGDS->iYDim] = nYOff;
547 57 : aiEdges[poGDS->iYDim] = nYSize;
548 :
549 57 : aiStart[poGDS->iXDim] = nBlockXOff;
550 57 : aiEdges[poGDS->iXDim] = nBlockXSize;
551 :
552 57 : if ( (SDwritedata( iSDS, aiStart, NULL,
553 : aiEdges, (VOIDP)pImage )) < 0 )
554 0 : eErr = CE_Failure;
555 :
556 57 : SDendaccess( iSDS );
557 : }
558 57 : break;
559 :
560 : case 2:
561 : {
562 8 : int32 iSDS = SDselect( poGDS->hSD, nBand - 1 );
563 8 : aiStart[poGDS->iYDim] = nYOff;
564 8 : aiEdges[poGDS->iYDim] = nYSize;
565 :
566 8 : aiStart[poGDS->iXDim] = nBlockXOff;
567 8 : aiEdges[poGDS->iXDim] = nBlockXSize;
568 :
569 8 : if ( (SDwritedata( iSDS, aiStart, NULL,
570 : aiEdges, (VOIDP)pImage )) < 0 )
571 0 : eErr = CE_Failure;
572 :
573 8 : SDendaccess( iSDS );
574 : }
575 8 : break;
576 :
577 : default:
578 0 : eErr = CE_Failure;
579 : break;
580 : }
581 :
582 65 : return eErr;
583 : }
584 :
585 : /************************************************************************/
586 : /* GetColorTable() */
587 : /************************************************************************/
588 :
589 0 : GDALColorTable *HDF4ImageRasterBand::GetColorTable()
590 : {
591 0 : HDF4ImageDataset *poGDS = (HDF4ImageDataset *) poDS;
592 :
593 0 : return poGDS->poColorTable;
594 : }
595 :
596 : /************************************************************************/
597 : /* GetColorInterpretation() */
598 : /************************************************************************/
599 :
600 36 : GDALColorInterp HDF4ImageRasterBand::GetColorInterpretation()
601 : {
602 36 : HDF4ImageDataset *poGDS = (HDF4ImageDataset *) poDS;
603 :
604 36 : if ( poGDS->iDatasetType == HDF4_SDS )
605 36 : return GCI_GrayIndex;
606 0 : else if ( poGDS->iDatasetType == HDF4_GR )
607 : {
608 0 : if ( poGDS->poColorTable != NULL )
609 0 : return GCI_PaletteIndex;
610 0 : else if ( poGDS->nBands != 1 )
611 : {
612 0 : if ( nBand == 1 )
613 0 : return GCI_RedBand;
614 0 : else if ( nBand == 2 )
615 0 : return GCI_GreenBand;
616 0 : else if ( nBand == 3 )
617 0 : return GCI_BlueBand;
618 0 : else if ( nBand == 4 )
619 0 : return GCI_AlphaBand;
620 : else
621 0 : return GCI_Undefined;
622 : }
623 : else
624 0 : return GCI_GrayIndex;
625 : }
626 : else
627 0 : return GCI_GrayIndex;
628 : }
629 :
630 : /************************************************************************/
631 : /* GetNoDataValue() */
632 : /************************************************************************/
633 :
634 96 : double HDF4ImageRasterBand::GetNoDataValue( int * pbSuccess )
635 :
636 : {
637 96 : if( pbSuccess )
638 96 : *pbSuccess = bNoDataSet;
639 :
640 96 : return dfNoDataValue;
641 : }
642 :
643 : /************************************************************************/
644 : /* SetNoDataValue() */
645 : /************************************************************************/
646 :
647 54 : CPLErr HDF4ImageRasterBand::SetNoDataValue( double dfNoData )
648 :
649 : {
650 54 : bNoDataSet = TRUE;
651 54 : dfNoDataValue = dfNoData;
652 :
653 54 : return CE_None;
654 : }
655 :
656 : /************************************************************************/
657 : /* GetUnitType() */
658 : /************************************************************************/
659 :
660 0 : const char *HDF4ImageRasterBand::GetUnitType()
661 :
662 : {
663 0 : if( osUnitType.size() > 0 )
664 0 : return osUnitType;
665 : else
666 0 : return GDALRasterBand::GetUnitType();
667 : }
668 :
669 : /************************************************************************/
670 : /* GetOffset() */
671 : /************************************************************************/
672 :
673 0 : double HDF4ImageRasterBand::GetOffset( int *pbSuccess )
674 :
675 : {
676 0 : if( bHaveScaleAndOffset )
677 : {
678 0 : if( pbSuccess != NULL )
679 0 : *pbSuccess = TRUE;
680 0 : return dfOffset;
681 : }
682 : else
683 0 : return GDALRasterBand::GetOffset( pbSuccess );
684 : }
685 :
686 : /************************************************************************/
687 : /* GetScale() */
688 : /************************************************************************/
689 :
690 0 : double HDF4ImageRasterBand::GetScale( int *pbSuccess )
691 :
692 : {
693 0 : if( bHaveScaleAndOffset )
694 : {
695 0 : if( pbSuccess != NULL )
696 0 : *pbSuccess = TRUE;
697 0 : return dfScale;
698 : }
699 : else
700 0 : return GDALRasterBand::GetScale( pbSuccess );
701 : }
702 :
703 : /************************************************************************/
704 : /* ==================================================================== */
705 : /* HDF4ImageDataset */
706 : /* ==================================================================== */
707 : /************************************************************************/
708 :
709 : /************************************************************************/
710 : /* HDF4ImageDataset() */
711 : /************************************************************************/
712 :
713 407 : HDF4ImageDataset::HDF4ImageDataset()
714 : {
715 407 : pszFilename = NULL;
716 407 : hHDF4 = 0;
717 407 : iGR = 0;
718 407 : iPal = 0;
719 407 : iDataset = 0;
720 407 : iRank = 0;
721 407 : iNumType = 0;
722 407 : nAttrs = 0;
723 407 : iInterlaceMode = 0;
724 407 : iPalInterlaceMode = 0;
725 407 : iPalDataType = 0;
726 407 : nComps = 0;
727 407 : nPalEntries = 0;
728 407 : memset(aiDimSizes, 0, sizeof(aiDimSizes));
729 407 : iXDim = 0;
730 407 : iYDim = 0;
731 407 : iBandDim = -1;
732 407 : i4Dim = 0;
733 407 : nBandCount = 0;
734 407 : papszLocalMetadata = NULL;
735 407 : memset(aiPaletteData, 0, sizeof(aiPaletteData));
736 407 : memset(szName, 0, sizeof(szName));
737 407 : pszSubdatasetName = NULL;
738 407 : pszFieldName = NULL;
739 407 : poColorTable = NULL;
740 407 : bHasGeoTransform = FALSE;
741 407 : adfGeoTransform[0] = 0.0;
742 407 : adfGeoTransform[1] = 1.0;
743 407 : adfGeoTransform[2] = 0.0;
744 407 : adfGeoTransform[3] = 0.0;
745 407 : adfGeoTransform[4] = 0.0;
746 407 : adfGeoTransform[5] = 1.0;
747 407 : pszProjection = CPLStrdup( "" );
748 407 : pszGCPProjection = CPLStrdup( "" );
749 407 : pasGCPList = NULL;
750 407 : nGCPCount = 0;
751 :
752 407 : iDatasetType = HDF4_UNKNOWN;
753 407 : iSDS = FAIL;
754 :
755 407 : nBlockPreferredXSize = -1;
756 407 : nBlockPreferredYSize = -1;
757 407 : bReadTile = false;
758 :
759 407 : }
760 :
761 : /************************************************************************/
762 : /* ~HDF4ImageDataset() */
763 : /************************************************************************/
764 :
765 407 : HDF4ImageDataset::~HDF4ImageDataset()
766 : {
767 407 : FlushCache();
768 :
769 407 : if ( pszFilename )
770 263 : CPLFree( pszFilename );
771 407 : if ( iSDS != FAIL )
772 44 : SDendaccess( iSDS );
773 407 : if ( hSD > 0 )
774 407 : SDend( hSD );
775 407 : hSD = 0;
776 407 : if ( iGR > 0 )
777 0 : GRendaccess( iGR );
778 407 : if ( hGR > 0 )
779 0 : GRend( hGR );
780 407 : hGR = 0;
781 407 : if ( pszSubdatasetName )
782 8 : CPLFree( pszSubdatasetName );
783 407 : if ( pszFieldName )
784 8 : CPLFree( pszFieldName );
785 407 : if ( papszLocalMetadata )
786 263 : CSLDestroy( papszLocalMetadata );
787 407 : if ( poColorTable != NULL )
788 0 : delete poColorTable;
789 407 : if ( pszProjection )
790 407 : CPLFree( pszProjection );
791 407 : if ( pszGCPProjection )
792 407 : CPLFree( pszGCPProjection );
793 407 : if( nGCPCount > 0 )
794 : {
795 0 : for( int i = 0; i < nGCPCount; i++ )
796 : {
797 0 : if ( pasGCPList[i].pszId )
798 0 : CPLFree( pasGCPList[i].pszId );
799 0 : if ( pasGCPList[i].pszInfo )
800 0 : CPLFree( pasGCPList[i].pszInfo );
801 : }
802 :
803 0 : CPLFree( pasGCPList );
804 : }
805 407 : if ( hHDF4 > 0 )
806 : {
807 263 : switch ( iDatasetType )
808 : {
809 : case HDF4_EOS:
810 8 : switch ( iSubdatasetType )
811 : {
812 : case H4ST_EOS_SWATH:
813 : case H4ST_EOS_SWATH_GEOL:
814 0 : SWclose( hHDF4 );
815 0 : break;
816 : case H4ST_EOS_GRID:
817 8 : GDclose( hHDF4 );
818 : default:
819 : break;
820 :
821 : }
822 8 : break;
823 : case HDF4_SDS:
824 : case HDF4_GR:
825 255 : hHDF4 = Hclose( hHDF4 );
826 : break;
827 : default:
828 : break;
829 : }
830 : }
831 407 : }
832 :
833 : /************************************************************************/
834 : /* GetGeoTransform() */
835 : /************************************************************************/
836 :
837 17 : CPLErr HDF4ImageDataset::GetGeoTransform( double * padfTransform )
838 : {
839 17 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
840 :
841 17 : if ( !bHasGeoTransform )
842 0 : return CE_Failure;
843 :
844 17 : 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 0 : int HDF4ImageDataset::GetGCPCount()
888 :
889 : {
890 0 : 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 407 : void HDF4ImageDataset::FlushCache()
920 :
921 : {
922 : int iBand;
923 : char *pszName;
924 : const char *pszValue;
925 :
926 407 : GDALDataset::FlushCache();
927 :
928 407 : if( eAccess == GA_ReadOnly )
929 263 : 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 47 : "Cannot write metadata information to output file");
967 : }
968 :
969 47 : CPLFree( pszName );
970 : }
971 : }
972 :
973 : // Write out NoData values
974 344 : for ( iBand = 1; iBand <= nBands; iBand++ )
975 : {
976 : HDF4ImageRasterBand *poBand =
977 200 : (HDF4ImageRasterBand *)GetRasterBand(iBand);
978 :
979 200 : if ( poBand->bNoDataSet )
980 : {
981 16 : pszName = CPLStrdup( CPLSPrintf( "NoDataValue%d", iBand ) );
982 16 : pszValue = CPLSPrintf( "%lf", 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 344 : for ( iBand = 1; iBand <= nBands; iBand++ )
997 : {
998 : HDF4ImageRasterBand *poBand =
999 200 : (HDF4ImageRasterBand *)GetRasterBand(iBand);
1000 :
1001 200 : pszName = CPLStrdup( CPLSPrintf( "BandDesc%d", iBand ) );
1002 200 : pszValue = poBand->GetDescription();
1003 200 : 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 200 : 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 : void HDF4ImageDataset::ReadCoordinates( const char *pszString,
1070 0 : 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 255 : void HDF4ImageDataset::CaptureL1GMTLInfo()
1145 :
1146 : {
1147 : /* -------------------------------------------------------------------- */
1148 : /* Does the physical file look like it matches our expected */
1149 : /* name pattern? */
1150 : /* -------------------------------------------------------------------- */
1151 255 : 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 0 : CPLString osMTLFilename = pszFilename;
1161 0 : osMTLFilename.resize(osMTLFilename.length() - 8);
1162 0 : osMTLFilename += "_MTL.L1G";
1163 :
1164 : /* -------------------------------------------------------------------- */
1165 : /* Ingest the MTL using the NASAKeywordHandler written for the */
1166 : /* PDS driver. */
1167 : /* -------------------------------------------------------------------- */
1168 0 : NASAKeywordHandler oMTL;
1169 :
1170 0 : FILE *fp = VSIFOpenL( osMTLFilename, "r" );
1171 0 : if( fp == NULL )
1172 0 : return;
1173 :
1174 0 : if( !oMTL.Ingest( fp, 0 ) )
1175 : {
1176 0 : VSIFCloseL( fp );
1177 : return;
1178 : }
1179 :
1180 0 : 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 0 : CPLString osPrefix;
1189 :
1190 0 : if( oMTL.GetKeyword( "LPGS_METADATA_FILE.PRODUCT_METADATA"
1191 : ".PRODUCT_UL_CORNER_LON", NULL ) )
1192 0 : osPrefix = "LPGS_METADATA_FILE.PRODUCT_METADATA.PRODUCT_";
1193 0 : else if( oMTL.GetKeyword( "L1_METADATA_FILE.PRODUCT_METADATA"
1194 : ".PRODUCT_UL_CORNER_LON", NULL ) )
1195 0 : osPrefix = "L1_METADATA_FILE.PRODUCT_METADATA.PRODUCT_";
1196 : else
1197 : return;
1198 :
1199 0 : dfULX = atof(oMTL.GetKeyword((osPrefix+"UL_CORNER_LON").c_str(), "0" ));
1200 0 : dfULY = atof(oMTL.GetKeyword((osPrefix+"UL_CORNER_LAT").c_str(), "0" ));
1201 0 : dfLRX = atof(oMTL.GetKeyword((osPrefix+"LR_CORNER_LON").c_str(), "0" ));
1202 0 : dfLRY = atof(oMTL.GetKeyword((osPrefix+"LR_CORNER_LAT").c_str(), "0" ));
1203 0 : dfLLX = atof(oMTL.GetKeyword((osPrefix+"LL_CORNER_LON").c_str(), "0" ));
1204 0 : dfLLY = atof(oMTL.GetKeyword((osPrefix+"LL_CORNER_LAT").c_str(), "0" ));
1205 0 : dfURX = atof(oMTL.GetKeyword((osPrefix+"UR_CORNER_LON").c_str(), "0" ));
1206 0 : dfURY = atof(oMTL.GetKeyword((osPrefix+"UR_CORNER_LAT").c_str(), "0" ));
1207 :
1208 0 : CPLFree( pszGCPProjection );
1209 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\"]]" );
1210 :
1211 0 : nGCPCount = 4;
1212 0 : pasGCPList = (GDAL_GCP *) CPLCalloc( nGCPCount, sizeof( GDAL_GCP ) );
1213 0 : GDALInitGCPs( nGCPCount, pasGCPList );
1214 :
1215 0 : pasGCPList[0].dfGCPX = dfULX;
1216 0 : pasGCPList[0].dfGCPY = dfULY;
1217 0 : pasGCPList[0].dfGCPPixel = 0.0;
1218 0 : pasGCPList[0].dfGCPLine = 0.0;
1219 :
1220 0 : pasGCPList[1].dfGCPX = dfURX;
1221 0 : pasGCPList[1].dfGCPY = dfURY;
1222 0 : pasGCPList[1].dfGCPPixel = GetRasterXSize();
1223 0 : pasGCPList[1].dfGCPLine = 0.0;
1224 :
1225 0 : pasGCPList[2].dfGCPX = dfLLX;
1226 0 : pasGCPList[2].dfGCPY = dfLLY;
1227 0 : pasGCPList[2].dfGCPPixel = 0.0;
1228 0 : pasGCPList[2].dfGCPLine = GetRasterYSize();
1229 :
1230 0 : pasGCPList[3].dfGCPX = dfLRX;
1231 0 : pasGCPList[3].dfGCPY = dfLRY;
1232 0 : pasGCPList[3].dfGCPPixel = GetRasterXSize();
1233 0 : 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 : || 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 : && aiDimSizes[0] >= 29
1379 : && SDreaddata( iSDS, aiStart, NULL, aiEdges, adfGCTP ) == 0
1380 : && oSRS.importFromUSGS( (long) adfGCTP[1], (long) adfGCTP[2],
1381 : adfGCTP+4,
1382 : (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 : && 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 : 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 : 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 : 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 : 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 255 : void HDF4ImageDataset::ProcessModisSDSGeolocation(void)
1877 :
1878 : {
1879 255 : int iDSIndex, iXIndex=-1, iYIndex=-1;
1880 :
1881 : // No point in assigning geolocation to the geolocation SDSes themselves.
1882 255 : if( EQUAL(szName,"longitude") || EQUAL(szName,"latitude") )
1883 0 : return;
1884 :
1885 255 : if (nRasterYSize == 1)
1886 0 : return;
1887 :
1888 : /* -------------------------------------------------------------------- */
1889 : /* Scan for latitude and longitude sections. */
1890 : /* -------------------------------------------------------------------- */
1891 : int32 nDatasets, nAttributes;
1892 :
1893 255 : if ( SDfileinfo( hSD, &nDatasets, &nAttributes ) != 0 )
1894 0 : return;
1895 :
1896 662 : 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 407 : iSDS = SDselect( hSD, iDSIndex );
1903 :
1904 407 : if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType,
1905 : &nAttrs) == 0 )
1906 : {
1907 407 : if( EQUAL(szName,"latitude") )
1908 2 : iYIndex = iDSIndex;
1909 :
1910 407 : if( EQUAL(szName,"longitude") )
1911 2 : iXIndex = iDSIndex;
1912 : }
1913 :
1914 407 : SDendaccess(iSDS);
1915 : }
1916 :
1917 255 : if( iXIndex == -1 || iYIndex == -1 )
1918 253 : 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 0 : char *pszGeoList = NULL;
1960 0 : char szGeoDimList[8192] = "";
1961 : int32 iWrkNumType;
1962 : int32 nDataFields, nDimMaps;
1963 0 : void *pLat = NULL, *pLong = NULL;
1964 0 : void *pLatticeX = NULL, *pLatticeY = NULL;
1965 0 : int32 iLatticeType, iLatticeDataSize = 0, iRank;
1966 0 : int32 nLatCount = 0, nLongCount = 0;
1967 0 : int32 nXPoints=0, nYPoints=0;
1968 : int32 nStrBufSize;
1969 : int32 aiDimSizes[H4_MAX_VAR_DIMS];
1970 0 : int i, j, iDataSize = 0, iPixelDim=-1,iLineDim=-1, iLongDim=-1, iLatDim=-1;
1971 0 : int32 *paiRank = NULL, *paiNumType = NULL,
1972 0 : *paiOffset = NULL, *paiIncrement = NULL;
1973 0 : char **papszGeolocations = NULL;
1974 0 : char *pszDimMaps = NULL;
1975 :
1976 : /* -------------------------------------------------------------------- */
1977 : /* Determine a product name. */
1978 : /* -------------------------------------------------------------------- */
1979 : const char *pszProduct =
1980 0 : CSLFetchNameValue( papszLocalMetadata, "SHORTNAME" );
1981 :
1982 0 : HDF4EOSProduct eProduct = PROD_UNKNOWN;
1983 0 : if ( pszProduct )
1984 : {
1985 0 : if ( EQUALN(pszProduct, "ASTL1A", 6) )
1986 0 : eProduct = PROD_ASTER_L1A;
1987 0 : else if ( EQUALN(pszProduct, "ASTL1B", 6) )
1988 0 : eProduct = PROD_ASTER_L1B;
1989 0 : else if ( EQUALN(pszProduct, "AST_04", 6)
1990 : || EQUALN(pszProduct, "AST_05", 6)
1991 : || EQUALN(pszProduct, "AST_06", 6)
1992 : || EQUALN(pszProduct, "AST_07", 6)
1993 : || EQUALN(pszProduct, "AST_08", 6)
1994 : || EQUALN(pszProduct, "AST_09", 6)
1995 : || EQUALN(pszProduct, "AST13", 5)
1996 : || EQUALN(pszProduct, "AST3", 4) )
1997 0 : eProduct = PROD_ASTER_L2;
1998 0 : else if ( EQUALN(pszProduct, "AST14", 5) )
1999 0 : eProduct = PROD_ASTER_L3;
2000 0 : else if ( EQUALN(pszProduct, "MOD02", 5)
2001 : || EQUALN(pszProduct, "MYD02", 5) )
2002 0 : eProduct = PROD_MODIS_L1B;
2003 : }
2004 :
2005 : /* -------------------------------------------------------------------- */
2006 : /* xx */
2007 : /* -------------------------------------------------------------------- */
2008 0 : nDataFields = SWnentries( hSW, HDFE_NENTGFLD, &nStrBufSize );
2009 0 : pszGeoList = (char *)CPLMalloc( nStrBufSize + 1 );
2010 0 : paiRank = (int32 *)CPLMalloc( nDataFields * sizeof(int32) );
2011 0 : paiNumType = (int32 *)CPLMalloc( nDataFields * sizeof(int32) );
2012 0 : if ( nDataFields !=
2013 : SWinqgeofields(hSW, pszGeoList, paiRank, paiNumType) )
2014 : {
2015 : CPLDebug( "HDF4Image",
2016 : "Can't get the list of geolocation fields in swath %s",
2017 0 : pszSubdatasetName );
2018 : }
2019 :
2020 : #ifdef DEBUG
2021 : else
2022 : {
2023 : char *pszTmp;
2024 : CPLDebug( "HDF4Image",
2025 : "Number of geolocation fields in swath %s: %ld",
2026 0 : pszSubdatasetName, (long)nDataFields );
2027 : CPLDebug( "HDF4Image",
2028 : "List of geolocation fields in swath %s: %s",
2029 0 : pszSubdatasetName, pszGeoList );
2030 : pszTmp = SPrintArray( GDT_UInt32, paiRank,
2031 0 : nDataFields, "," );
2032 : CPLDebug( "HDF4Image",
2033 0 : "Geolocation fields ranks: %s", pszTmp );
2034 0 : CPLFree( pszTmp );
2035 : }
2036 : #endif
2037 :
2038 : papszGeolocations = CSLTokenizeString2( pszGeoList, ",",
2039 0 : CSLT_HONOURSTRINGS );
2040 : // Read geolocation data
2041 0 : nDimMaps = SWnentries( hSW, HDFE_NENTMAP, &nStrBufSize );
2042 0 : if ( nDimMaps <= 0 )
2043 : {
2044 : CPLDebug( "HDF4Image", "No geolocation maps in swath %s",
2045 0 : pszSubdatasetName );
2046 : }
2047 : else
2048 : {
2049 0 : pszDimMaps = (char *)CPLMalloc( nStrBufSize + 1 );
2050 0 : paiOffset = (int32 *)CPLMalloc( nDimMaps * sizeof(int32) );
2051 0 : memset( paiOffset, 0, nDimMaps * sizeof(int32) );
2052 0 : paiIncrement = (int32 *)CPLMalloc( nDimMaps * sizeof(int32) );
2053 0 : memset( paiIncrement, 0, nDimMaps * sizeof(int32) );
2054 :
2055 :
2056 0 : *pszDimMaps = '\0';
2057 0 : if ( nDimMaps != SWinqmaps(hSW, pszDimMaps, paiOffset, paiIncrement) )
2058 : {
2059 : CPLDebug( "HDF4Image",
2060 : "Can't get the list of geolocation maps in swath %s",
2061 0 : pszSubdatasetName );
2062 : }
2063 :
2064 : #ifdef DEBUG
2065 : else
2066 : {
2067 : char *pszTmp;
2068 :
2069 : CPLDebug( "HDF4Image",
2070 : "List of geolocation maps in swath %s: %s",
2071 0 : pszSubdatasetName, pszDimMaps );
2072 : pszTmp = SPrintArray( GDT_Int32, paiOffset,
2073 0 : nDimMaps, "," );
2074 : CPLDebug( "HDF4Image",
2075 0 : "Geolocation map offsets: %s", pszTmp );
2076 0 : CPLFree( pszTmp );
2077 : pszTmp = SPrintArray( GDT_Int32, paiIncrement,
2078 0 : nDimMaps, "," );
2079 : CPLDebug( "HDF4Image",
2080 0 : "Geolocation map increments: %s", pszTmp );
2081 0 : CPLFree( pszTmp );
2082 : }
2083 : #endif
2084 :
2085 : char **papszDimMap;
2086 :
2087 : papszDimMap = CSLTokenizeString2( pszDimMaps, ",",
2088 0 : CSLT_HONOURSTRINGS );
2089 0 : CPLFree( pszDimMaps );
2090 0 : pszDimMaps = NULL;
2091 :
2092 0 : for ( i = 0; i < CSLCount(papszDimMap); i++ )
2093 : {
2094 0 : if ( strstr(papszDimMap[i], papszDimList[iXDim]) )
2095 : {
2096 0 : strncpy( szPixel, papszDimList[iXDim], 8192 );
2097 0 : strncpy( szXGeo, papszDimMap[i], 8192 );
2098 0 : char *pszTemp = strchr( szXGeo, '/' );
2099 0 : if ( pszTemp )
2100 0 : *pszTemp = '\0';
2101 : }
2102 0 : else if ( strstr(papszDimMap[i], papszDimList[iYDim]) )
2103 : {
2104 0 : strncpy( szLine, papszDimList[iYDim], 8192 );
2105 0 : strncpy( szYGeo, papszDimMap[i], 8192 );
2106 0 : char *pszTemp = strchr( szYGeo, '/' );
2107 0 : if ( pszTemp )
2108 0 : *pszTemp = '\0';
2109 : }
2110 : }
2111 0 : CSLDestroy( papszDimMap );
2112 0 : papszDimMap = NULL;
2113 : }
2114 :
2115 0 : for ( i = 0; i < CSLCount(papszGeolocations); i++ )
2116 : {
2117 0 : char **papszGeoDimList = NULL;
2118 :
2119 : SWfieldinfo( hSW, papszGeolocations[i], &iRank,
2120 0 : aiDimSizes, &iWrkNumType, szGeoDimList );
2121 : papszGeoDimList = CSLTokenizeString2( szGeoDimList,
2122 0 : ",", CSLT_HONOURSTRINGS );
2123 :
2124 : #ifdef DEBUG
2125 : CPLDebug( "HDF4Image",
2126 : "List of dimensions in geolocation field %s: %s",
2127 0 : papszGeolocations[i], szGeoDimList );
2128 : #endif
2129 :
2130 0 : if (szXGeo[0] == 0 || szYGeo[0] == 0)
2131 0 : return FALSE;
2132 :
2133 0 : nXPoints = aiDimSizes[CSLFindString( papszGeoDimList, szXGeo )];
2134 0 : nYPoints = aiDimSizes[CSLFindString( papszGeoDimList, szYGeo )];
2135 :
2136 0 : if ( EQUAL(szPixel, papszDimList[iXDim]) )
2137 : {
2138 0 : iPixelDim = 1;
2139 0 : iLineDim = 0;
2140 : }
2141 : else
2142 : {
2143 0 : iPixelDim = 0;
2144 0 : iLineDim = 1;
2145 : }
2146 :
2147 0 : iDataSize = GetDataTypeSize( iWrkNumType );
2148 0 : if ( strstr( papszGeolocations[i], "Latitude" ) )
2149 : {
2150 0 : iLatDim = i;
2151 0 : nLatCount = nXPoints * nYPoints;
2152 0 : pLat = CPLMalloc( nLatCount * iDataSize );
2153 0 : if (SWreadfield( hSW, papszGeolocations[i], NULL,
2154 : NULL, NULL, (VOIDP)pLat ) < 0)
2155 : {
2156 : CPLDebug( "HDF4Image",
2157 : "Can't read geolocation field %s",
2158 0 : papszGeolocations[i]);
2159 0 : CPLFree( pLat );
2160 0 : pLat = NULL;
2161 : }
2162 : }
2163 0 : else if ( strstr( papszGeolocations[i], "Longitude" ) )
2164 : {
2165 0 : iLongDim = i;
2166 0 : nLongCount = nXPoints * nYPoints;
2167 0 : pLong = CPLMalloc( nLongCount * iDataSize );
2168 0 : if (SWreadfield( hSW, papszGeolocations[i], NULL,
2169 : NULL, NULL, (VOIDP)pLong ) < 0)
2170 : {
2171 : CPLDebug( "HDF4Image",
2172 : "Can't read geolocation field %s",
2173 0 : papszGeolocations[i]);
2174 0 : CPLFree( pLong );
2175 0 : pLong = NULL;
2176 : }
2177 : }
2178 :
2179 0 : CSLDestroy( papszGeoDimList );
2180 :
2181 : }
2182 :
2183 : /* -------------------------------------------------------------------- */
2184 : /* Do we have a lattice table? */
2185 : /* -------------------------------------------------------------------- */
2186 0 : if (SWfieldinfo(hSW, "LatticePoint", &iRank, aiDimSizes,
2187 : &iLatticeType, szGeoDimList) == 0
2188 : && iRank == 3
2189 : && nXPoints == aiDimSizes[1]
2190 : && nYPoints == aiDimSizes[0]
2191 : && aiDimSizes[2] == 2 )
2192 : {
2193 : int32 iStart[H4_MAX_NC_DIMS], iEdges[H4_MAX_NC_DIMS];
2194 :
2195 : iLatticeDataSize =
2196 0 : GetDataTypeSize( iLatticeType );
2197 :
2198 0 : iStart[1] = 0;
2199 0 : iEdges[1] = nXPoints;
2200 :
2201 0 : iStart[0] = 0;
2202 0 : iEdges[0] = nYPoints;
2203 :
2204 0 : iStart[2] = 0;
2205 0 : iEdges[2] = 1;
2206 :
2207 0 : pLatticeX = CPLMalloc( nLatCount * iLatticeDataSize );
2208 0 : if (SWreadfield( hSW, "LatticePoint", iStart, NULL,
2209 : iEdges, (VOIDP)pLatticeX ) < 0)
2210 : {
2211 0 : CPLDebug( "HDF4Image", "Can't read lattice field" );
2212 0 : CPLFree( pLatticeX );
2213 0 : pLatticeX = NULL;
2214 : }
2215 :
2216 0 : iStart[2] = 1;
2217 0 : iEdges[2] = 1;
2218 :
2219 0 : pLatticeY = CPLMalloc( nLatCount * iLatticeDataSize );
2220 0 : if (SWreadfield( hSW, "LatticePoint", iStart, NULL,
2221 : iEdges, (VOIDP)pLatticeY ) < 0)
2222 : {
2223 0 : CPLDebug( "HDF4Image", "Can't read lattice field" );
2224 0 : CPLFree( pLatticeY );
2225 0 : pLatticeY = NULL;
2226 : }
2227 :
2228 : }
2229 :
2230 : /* -------------------------------------------------------------------- */
2231 : /* Determine whether to use no, partial or full GCPs. */
2232 : /* -------------------------------------------------------------------- */
2233 : const char *pszGEOL_AS_GCPS = CPLGetConfigOption( "GEOL_AS_GCPS",
2234 0 : "PARTIAL" );
2235 : int iGCPStepX, iGCPStepY;
2236 :
2237 0 : if( EQUAL(pszGEOL_AS_GCPS,"NONE") )
2238 : {
2239 0 : iGCPStepX = iGCPStepY = 0;
2240 : }
2241 0 : else if( EQUAL(pszGEOL_AS_GCPS,"FULL") )
2242 : {
2243 0 : iGCPStepX = iGCPStepY = 1;
2244 : }
2245 : else
2246 : {
2247 : // aim for 10x10 grid or so.
2248 0 : iGCPStepX = MAX(1,((nXPoints-1) / 11));
2249 0 : iGCPStepY = MAX(1,((nYPoints-1) / 11));
2250 : }
2251 :
2252 : /* -------------------------------------------------------------------- */
2253 : /* Fetch projection information for various datasets. */
2254 : /* -------------------------------------------------------------------- */
2255 0 : if ( nLatCount && nLongCount && nLatCount == nLongCount
2256 : && pLat && pLong )
2257 : {
2258 0 : CPLFree( pszGCPProjection );
2259 0 : pszGCPProjection = NULL;
2260 :
2261 : // ASTER Level 1A
2262 0 : if ( eProduct == PROD_ASTER_L1A )
2263 : {
2264 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\"]]" );
2265 : }
2266 :
2267 : // ASTER Level 1B, Level 2
2268 0 : else if ( eProduct == PROD_ASTER_L1B
2269 : || eProduct == PROD_ASTER_L2 )
2270 : {
2271 : // Constuct the metadata keys.
2272 : // A band number is taken from the field name.
2273 0 : const char *pszBand = strpbrk( pszFieldName, "0123456789" );
2274 :
2275 0 : if ( !pszBand )
2276 0 : pszBand = "";
2277 :
2278 : char *pszProjLine =
2279 0 : CPLStrdup(CPLSPrintf("MPMETHOD%s", pszBand));
2280 : char *pszParmsLine =
2281 : CPLStrdup(CPLSPrintf("PROJECTIONPARAMETERS%s",
2282 0 : pszBand));
2283 : char *pszZoneLine =
2284 : CPLStrdup(CPLSPrintf("UTMZONECODE%s",
2285 0 : pszBand));
2286 : char *pszEllipsoidLine =
2287 : CPLStrdup(CPLSPrintf("ELLIPSOIDANDDATUM%s",
2288 0 : pszBand));
2289 :
2290 : // Fetch projection related values from the
2291 : // metadata.
2292 : const char *pszProj =
2293 : CSLFetchNameValue( papszLocalMetadata,
2294 0 : pszProjLine );
2295 : const char *pszParms =
2296 : CSLFetchNameValue( papszLocalMetadata,
2297 0 : pszParmsLine );
2298 : const char *pszZone =
2299 : CSLFetchNameValue( papszLocalMetadata,
2300 0 : pszZoneLine );
2301 : const char* pszEllipsoid =
2302 : CSLFetchNameValue( papszLocalMetadata,
2303 0 : pszEllipsoidLine );
2304 :
2305 : #ifdef DEBUG
2306 : CPLDebug( "HDF4Image",
2307 : "Projection %s=%s, parameters %s=%s, "
2308 : "zone %s=%s",
2309 : pszProjLine, pszProj, pszParmsLine,
2310 0 : pszParms, pszZoneLine, pszZone );
2311 : CPLDebug( "HDF4Image", "Ellipsoid %s=%s",
2312 0 : pszEllipsoidLine, pszEllipsoid );
2313 : #endif
2314 :
2315 : // Transform all mnemonical codes in the values.
2316 : int i, nParms;
2317 : // Projection is UTM by default
2318 : long iProjSys = (pszProj) ?
2319 0 : USGSMnemonicToCode(pszProj) : 1L;
2320 : long iZone =
2321 0 : (pszZone && iProjSys == 1L) ? atoi(pszZone): 0L;
2322 : char **papszEllipsoid = (pszEllipsoid) ?
2323 : CSLTokenizeString2( pszEllipsoid, ",",
2324 0 : CSLT_HONOURSTRINGS ) : NULL;
2325 :
2326 0 : long iEllipsoid = 8L; // WGS84 by default
2327 0 : if ( papszEllipsoid
2328 : && CSLCount(papszEllipsoid) > 0 )
2329 : {
2330 0 : if (EQUAL( papszEllipsoid[0], "WGS84"))
2331 0 : iEllipsoid = 8L;
2332 : }
2333 :
2334 : double adfProjParms[15];
2335 : char **papszParms = (pszParms) ?
2336 : CSLTokenizeString2( pszParms, ",",
2337 0 : CSLT_HONOURSTRINGS ) : NULL;
2338 0 : nParms = CSLCount(papszParms);
2339 0 : if (nParms >= 15)
2340 0 : nParms = 15;
2341 0 : for (i = 0; i < nParms; i++)
2342 0 : adfProjParms[i] = CPLAtof( papszParms[i] );
2343 0 : for (; i < 15; i++)
2344 0 : adfProjParms[i] = 0.0;
2345 :
2346 : // Create projection definition
2347 : oSRS.importFromUSGS( iProjSys, iZone,
2348 0 : adfProjParms, iEllipsoid );
2349 0 : oSRS.SetLinearUnits( SRS_UL_METER, 1.0 );
2350 0 : oSRS.exportToWkt( &pszGCPProjection );
2351 :
2352 0 : CSLDestroy( papszParms );
2353 0 : CPLFree( pszEllipsoidLine );
2354 0 : CPLFree( pszZoneLine );
2355 0 : CPLFree( pszParmsLine );
2356 0 : CPLFree( pszProjLine );
2357 : }
2358 :
2359 : // ASTER Level 3 (DEM)
2360 0 : else if ( eProduct == PROD_ASTER_L3 )
2361 : {
2362 : double dfCenterX, dfCenterY;
2363 : int iZone;
2364 :
2365 : ReadCoordinates( CSLFetchNameValue(
2366 : papszGlobalMetadata, "SCENECENTER" ),
2367 0 : &dfCenterY, &dfCenterX );
2368 :
2369 : // Calculate UTM zone from scene center coordinates
2370 0 : iZone = 30 + (int) ((dfCenterX + 6.0) / 6.0);
2371 :
2372 : // Create projection definition
2373 0 : if( dfCenterY > 0 )
2374 0 : oSRS.SetUTM( iZone, TRUE );
2375 : else
2376 0 : oSRS.SetUTM( - iZone, FALSE );
2377 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
2378 0 : oSRS.SetLinearUnits( SRS_UL_METER, 1.0 );
2379 0 : oSRS.exportToWkt( &pszGCPProjection );
2380 : }
2381 :
2382 : // MODIS L1B
2383 0 : else if ( eProduct == PROD_MODIS_L1B )
2384 : {
2385 0 : pszGCPProjection = CPLStrdup( SRS_WKT_WGS84 );
2386 : }
2387 : }
2388 :
2389 : /* -------------------------------------------------------------------- */
2390 : /* Fill the GCPs list. */
2391 : /* -------------------------------------------------------------------- */
2392 0 : if( iGCPStepX > 0 )
2393 : {
2394 : nGCPCount = (((nXPoints-1) / iGCPStepX) + 1)
2395 0 : * (((nYPoints-1) / iGCPStepY) + 1);
2396 :
2397 0 : pasGCPList = (GDAL_GCP *) CPLCalloc( nGCPCount, sizeof( GDAL_GCP ) );
2398 0 : GDALInitGCPs( nGCPCount, pasGCPList );
2399 :
2400 0 : int iGCP = 0;
2401 0 : for ( i = 0; i < nYPoints; i += iGCPStepY )
2402 : {
2403 0 : for ( j = 0; j < nXPoints; j += iGCPStepX )
2404 : {
2405 0 : int iGeoOff = i * nXPoints + j;
2406 :
2407 : pasGCPList[iGCP].dfGCPX =
2408 : AnyTypeToDouble(iWrkNumType,
2409 0 : (void *)((char*)pLong+ iGeoOff*iDataSize));
2410 : pasGCPList[iGCP].dfGCPY =
2411 : AnyTypeToDouble(iWrkNumType,
2412 0 : (void *)((char*)pLat + iGeoOff*iDataSize));
2413 :
2414 : // GCPs in Level 1A/1B dataset are in geocentric
2415 : // coordinates. Convert them in geodetic (we
2416 : // will convert latitudes only, longitudes
2417 : // do not need to be converted, because
2418 : // they are the same).
2419 : // This calculation valid for WGS84 datum only.
2420 0 : if ( eProduct == PROD_ASTER_L1A
2421 : || eProduct == PROD_ASTER_L1B )
2422 : {
2423 : pasGCPList[iGCP].dfGCPY =
2424 : atan(tan(pasGCPList[iGCP].dfGCPY
2425 0 : *PI/180)/0.99330562)*180/PI;
2426 : }
2427 :
2428 : ToGeoref(&pasGCPList[iGCP].dfGCPX,
2429 0 : &pasGCPList[iGCP].dfGCPY);
2430 :
2431 0 : pasGCPList[iGCP].dfGCPZ = 0.0;
2432 :
2433 0 : if ( pLatticeX && pLatticeY )
2434 : {
2435 : pasGCPList[iGCP].dfGCPPixel =
2436 : AnyTypeToDouble(iLatticeType,
2437 : (void *)((char *)pLatticeX
2438 0 : + iGeoOff*iLatticeDataSize))+0.5;
2439 : pasGCPList[iGCP].dfGCPLine =
2440 : AnyTypeToDouble(iLatticeType,
2441 : (void *)((char *)pLatticeY
2442 0 : + iGeoOff*iLatticeDataSize))+0.5;
2443 : }
2444 0 : else if ( paiOffset && paiIncrement )
2445 : {
2446 : pasGCPList[iGCP].dfGCPPixel =
2447 : paiOffset[iPixelDim] +
2448 0 : j * paiIncrement[iPixelDim] + 0.5;
2449 : pasGCPList[iGCP].dfGCPLine =
2450 : paiOffset[iLineDim] +
2451 0 : i * paiIncrement[iLineDim] + 0.5;
2452 : }
2453 :
2454 0 : iGCP++;
2455 : }
2456 : }
2457 : }
2458 :
2459 : /* -------------------------------------------------------------------- */
2460 : /* Establish geolocation metadata, but only if there is no */
2461 : /* lattice. The lattice destroys the regularity of the grid. */
2462 : /* -------------------------------------------------------------------- */
2463 0 : if( pLatticeX == NULL
2464 : && iLatDim != -1 && iLongDim != -1
2465 : && iPixelDim != -1 && iLineDim != -1 )
2466 : {
2467 0 : CPLString osWrk;
2468 :
2469 0 : SetMetadataItem( "SRS", pszGCPProjection, "GEOLOCATION" );
2470 :
2471 : osWrk.Printf( "HDF4_EOS:EOS_SWATH_GEOL:\"%s\":%s:%s",
2472 : pszFilename, pszSubdatasetName,
2473 0 : papszGeolocations[iLongDim] );
2474 0 : SetMetadataItem( "X_DATASET", osWrk, "GEOLOCATION" );
2475 0 : SetMetadataItem( "X_BAND", "1" , "GEOLOCATION" );
2476 :
2477 : osWrk.Printf( "HDF4_EOS:EOS_SWATH_GEOL:\"%s\":%s:%s",
2478 : pszFilename, pszSubdatasetName,
2479 0 : papszGeolocations[iLatDim] );
2480 0 : SetMetadataItem( "Y_DATASET", osWrk, "GEOLOCATION" );
2481 0 : SetMetadataItem( "Y_BAND", "1" , "GEOLOCATION" );
2482 :
2483 0 : if ( paiOffset && paiIncrement )
2484 : {
2485 0 : osWrk.Printf( "%ld", (long)paiOffset[iPixelDim] );
2486 0 : SetMetadataItem( "PIXEL_OFFSET", osWrk, "GEOLOCATION" );
2487 0 : osWrk.Printf( "%ld", (long)paiIncrement[iPixelDim] );
2488 0 : SetMetadataItem( "PIXEL_STEP", osWrk, "GEOLOCATION" );
2489 :
2490 0 : osWrk.Printf( "%ld", (long)paiOffset[iLineDim] );
2491 0 : SetMetadataItem( "LINE_OFFSET", osWrk, "GEOLOCATION" );
2492 0 : osWrk.Printf( "%ld", (long)paiIncrement[iLineDim] );
2493 0 : SetMetadataItem( "LINE_STEP", osWrk, "GEOLOCATION" );
2494 0 : }
2495 : }
2496 :
2497 : /* -------------------------------------------------------------------- */
2498 : /* Cleanup */
2499 : /* -------------------------------------------------------------------- */
2500 0 : CPLFree( pLatticeX );
2501 0 : CPLFree( pLatticeY );
2502 0 : CPLFree( pLat );
2503 0 : CPLFree( pLong );
2504 :
2505 0 : CPLFree( paiOffset );
2506 0 : CPLFree( paiIncrement );
2507 0 : CPLFree( paiNumType );
2508 0 : CPLFree( paiRank );
2509 :
2510 0 : CSLDestroy( papszGeolocations );
2511 0 : CPLFree( pszGeoList );
2512 :
2513 0 : if( iGCPStepX == 0 )
2514 : {
2515 0 : CPLFree( pszGCPProjection );
2516 0 : pszGCPProjection = NULL;
2517 : }
2518 :
2519 0 : return TRUE;
2520 : }
2521 :
2522 : /************************************************************************/
2523 : /* Open() */
2524 : /************************************************************************/
2525 :
2526 10562 : GDALDataset *HDF4ImageDataset::Open( GDALOpenInfo * poOpenInfo )
2527 : {
2528 : int i;
2529 :
2530 10562 : if( !EQUALN( poOpenInfo->pszFilename, "HDF4_SDS:", 9 ) &&
2531 : !EQUALN( poOpenInfo->pszFilename, "HDF4_GR:", 8 ) &&
2532 : !EQUALN( poOpenInfo->pszFilename, "HDF4_GD:", 8 ) &&
2533 : !EQUALN( poOpenInfo->pszFilename, "HDF4_EOS:", 9 ) )
2534 10299 : return NULL;
2535 :
2536 : /* -------------------------------------------------------------------- */
2537 : /* Create a corresponding GDALDataset. */
2538 : /* -------------------------------------------------------------------- */
2539 : char **papszSubdatasetName;
2540 : HDF4ImageDataset *poDS;
2541 :
2542 263 : poDS = new HDF4ImageDataset( );
2543 263 : poDS->fp = poOpenInfo->fp;
2544 263 : poOpenInfo->fp = NULL;
2545 :
2546 : papszSubdatasetName = CSLTokenizeString2( poOpenInfo->pszFilename,
2547 263 : ":", CSLT_HONOURSTRINGS | CSLT_PRESERVEESCAPES);
2548 271 : if ( CSLCount( papszSubdatasetName ) != 4
2549 : && CSLCount( papszSubdatasetName ) != 5
2550 : && CSLCount( papszSubdatasetName ) != 6 )
2551 : {
2552 0 : CSLDestroy( papszSubdatasetName );
2553 0 : delete poDS;
2554 0 : return NULL;
2555 : }
2556 :
2557 : /* -------------------------------------------------------------------- */
2558 : /* Check for drive name in windows HDF4:"D:\... */
2559 : /* -------------------------------------------------------------------- */
2560 263 : if (strlen(papszSubdatasetName[2]) == 1)
2561 : {
2562 0 : char* pszFilename = (char*) CPLMalloc( 2 + strlen(papszSubdatasetName[3]) + 1);
2563 0 : sprintf(pszFilename, "%s:%s", papszSubdatasetName[2], papszSubdatasetName[3]);
2564 0 : CPLFree(papszSubdatasetName[2]);
2565 0 : CPLFree(papszSubdatasetName[3]);
2566 0 : papszSubdatasetName[2] = pszFilename;
2567 :
2568 : /* Move following arguments one rank upper */
2569 0 : papszSubdatasetName[3] = papszSubdatasetName[4];
2570 0 : if (papszSubdatasetName[4] != NULL)
2571 : {
2572 0 : papszSubdatasetName[4] = papszSubdatasetName[5];
2573 0 : papszSubdatasetName[5] = NULL;
2574 : }
2575 : }
2576 :
2577 263 : poDS->pszFilename = CPLStrdup( papszSubdatasetName[2] );
2578 :
2579 263 : if( EQUAL( papszSubdatasetName[0], "HDF4_SDS" ) )
2580 255 : poDS->iDatasetType = HDF4_SDS;
2581 8 : else if ( EQUAL( papszSubdatasetName[0], "HDF4_GR" ) )
2582 0 : poDS->iDatasetType = HDF4_GR;
2583 8 : else if ( EQUAL( papszSubdatasetName[0], "HDF4_EOS" ) )
2584 8 : poDS->iDatasetType = HDF4_EOS;
2585 : else
2586 0 : poDS->iDatasetType = HDF4_UNKNOWN;
2587 :
2588 263 : if( EQUAL( papszSubdatasetName[1], "GDAL_HDF4" ) )
2589 249 : poDS->iSubdatasetType = H4ST_GDAL;
2590 14 : else if( EQUAL( papszSubdatasetName[1], "EOS_GRID" ) )
2591 8 : poDS->iSubdatasetType = H4ST_EOS_GRID;
2592 6 : else if( EQUAL( papszSubdatasetName[1], "EOS_SWATH" ) )
2593 0 : poDS->iSubdatasetType = H4ST_EOS_SWATH;
2594 6 : else if( EQUAL( papszSubdatasetName[1], "EOS_SWATH_GEOL" ) )
2595 0 : poDS->iSubdatasetType = H4ST_EOS_SWATH_GEOL;
2596 6 : else if( EQUAL( papszSubdatasetName[1], "SEAWIFS_L3" ) )
2597 0 : poDS->iSubdatasetType= H4ST_SEAWIFS_L3;
2598 6 : else if( EQUAL( papszSubdatasetName[1], "HYPERION_L1" ) )
2599 0 : poDS->iSubdatasetType= H4ST_HYPERION_L1;
2600 : else
2601 6 : poDS->iSubdatasetType = H4ST_UNKNOWN;
2602 :
2603 : // Is our file still here?
2604 263 : if ( !Hishdf( poDS->pszFilename ) )
2605 : {
2606 0 : CSLDestroy( papszSubdatasetName );
2607 0 : delete poDS;
2608 0 : return NULL;
2609 : }
2610 :
2611 : /* -------------------------------------------------------------------- */
2612 : /* Collect the remain (post filename) components to treat as */
2613 : /* the subdataset name. */
2614 : /* -------------------------------------------------------------------- */
2615 263 : CPLString osSubdatasetName;
2616 :
2617 263 : osSubdatasetName = papszSubdatasetName[3];
2618 263 : if( papszSubdatasetName[4] != NULL )
2619 : {
2620 8 : osSubdatasetName += ":";
2621 8 : osSubdatasetName += papszSubdatasetName[4];
2622 : }
2623 :
2624 : /* -------------------------------------------------------------------- */
2625 : /* Try opening the dataset. */
2626 : /* -------------------------------------------------------------------- */
2627 : int32 iAttribute, nValues, iAttrNumType;
2628 263 : double dfNoData = 0.0;
2629 263 : int bNoDataSet = FALSE;
2630 :
2631 : /* -------------------------------------------------------------------- */
2632 : /* Select SDS or GR to read from. */
2633 : /* -------------------------------------------------------------------- */
2634 263 : if ( poDS->iDatasetType == HDF4_EOS )
2635 : {
2636 8 : poDS->pszSubdatasetName = CPLStrdup( papszSubdatasetName[3] );
2637 8 : if (papszSubdatasetName[4] == NULL)
2638 : {
2639 0 : delete poDS;
2640 263 : return NULL;
2641 : }
2642 8 : poDS->pszFieldName = CPLStrdup( papszSubdatasetName[4] );
2643 : }
2644 : else
2645 255 : poDS->iDataset = atoi( papszSubdatasetName[3] );
2646 263 : CSLDestroy( papszSubdatasetName );
2647 :
2648 263 : switch ( poDS->iDatasetType )
2649 : {
2650 : case HDF4_EOS:
2651 : {
2652 8 : void *pNoDataValue = NULL;
2653 :
2654 8 : switch ( poDS->iSubdatasetType )
2655 : {
2656 :
2657 : /* -------------------------------------------------------------------- */
2658 : /* HDF-EOS Swath. */
2659 : /* -------------------------------------------------------------------- */
2660 : case H4ST_EOS_SWATH:
2661 : case H4ST_EOS_SWATH_GEOL:
2662 : {
2663 : int32 hSW, nStrBufSize;
2664 0 : char *pszDimList = NULL;
2665 :
2666 0 : if( poOpenInfo->eAccess == GA_ReadOnly )
2667 0 : poDS->hHDF4 = SWopen( poDS->pszFilename, DFACC_READ );
2668 : else
2669 0 : poDS->hHDF4 = SWopen( poDS->pszFilename, DFACC_WRITE );
2670 :
2671 0 : if( poDS->hHDF4 <= 0 )
2672 : {
2673 : CPLDebug( "HDF4Image", "Can't open HDF4 file %s",
2674 0 : poDS->pszFilename );
2675 0 : delete poDS;
2676 0 : return( NULL );
2677 : }
2678 :
2679 0 : hSW = SWattach( poDS->hHDF4, poDS->pszSubdatasetName );
2680 0 : if( hSW < 0 )
2681 : {
2682 : CPLDebug( "HDF4Image", "Can't attach to subdataset %s",
2683 0 : poDS->pszSubdatasetName );
2684 0 : delete poDS;
2685 0 : return( NULL );
2686 : }
2687 :
2688 : /* -------------------------------------------------------------------- */
2689 : /* Decode the dimension map. */
2690 : /* -------------------------------------------------------------------- */
2691 0 : SWnentries( hSW, HDFE_NENTDIM, &nStrBufSize );
2692 0 : pszDimList = (char *)CPLMalloc( nStrBufSize + 1 );
2693 : SWfieldinfo( hSW, poDS->pszFieldName, &poDS->iRank,
2694 0 : poDS->aiDimSizes, &poDS->iNumType, pszDimList );
2695 : #ifdef DEBUG
2696 : CPLDebug( "HDF4Image",
2697 : "List of dimensions in swath %s: %s",
2698 0 : poDS->pszFieldName, pszDimList );
2699 : #endif
2700 0 : poDS->GetImageDimensions( pszDimList );
2701 :
2702 : #ifdef DEBUG
2703 : CPLDebug( "HDF4Image",
2704 : "X dimension is %d, Y dimension is %d",
2705 0 : poDS->iXDim, poDS->iYDim );
2706 : #endif
2707 :
2708 : /* -------------------------------------------------------------------- */
2709 : /* Fetch metadata. */
2710 : /* -------------------------------------------------------------------- */
2711 0 : if( poDS->iSubdatasetType == H4ST_EOS_SWATH ) /* Not H4ST_EOS_SWATH_GEOL */
2712 0 : poDS->GetSwatAttrs( hSW );
2713 :
2714 : /* -------------------------------------------------------------------- */
2715 : /* Fetch NODATA value. */
2716 : /* -------------------------------------------------------------------- */
2717 : pNoDataValue =
2718 0 : CPLMalloc( poDS->GetDataTypeSize(poDS->iNumType) );
2719 0 : if ( SWgetfillvalue( hSW, poDS->pszFieldName,
2720 : pNoDataValue ) != -1 )
2721 : {
2722 : dfNoData = poDS->AnyTypeToDouble( poDS->iNumType,
2723 0 : pNoDataValue );
2724 0 : bNoDataSet = TRUE;
2725 : }
2726 : else
2727 : {
2728 : const char *pszNoData =
2729 : CSLFetchNameValue( poDS->papszLocalMetadata,
2730 0 : "_FillValue" );
2731 0 : if ( pszNoData )
2732 : {
2733 0 : dfNoData = CPLAtof( pszNoData );
2734 0 : bNoDataSet = TRUE;
2735 : }
2736 : }
2737 0 : CPLFree( pNoDataValue );
2738 :
2739 : /* -------------------------------------------------------------------- */
2740 : /* Handle Geolocation processing. */
2741 : /* -------------------------------------------------------------------- */
2742 0 : if( poDS->iSubdatasetType == H4ST_EOS_SWATH ) /* Not H4ST_SWATH_GEOL */
2743 : {
2744 : char **papszDimList =
2745 : CSLTokenizeString2( pszDimList, ",",
2746 0 : CSLT_HONOURSTRINGS );
2747 0 : if( !poDS->ProcessSwathGeolocation( hSW, papszDimList ) )
2748 : {
2749 : CPLDebug( "HDF4Image",
2750 0 : "No geolocation available for this swath." );
2751 : }
2752 0 : CSLDestroy( papszDimList );
2753 : }
2754 :
2755 : /* -------------------------------------------------------------------- */
2756 : /* Cleanup. */
2757 : /* -------------------------------------------------------------------- */
2758 0 : CPLFree( pszDimList );
2759 0 : SWdetach( hSW );
2760 : }
2761 0 : break;
2762 :
2763 : /* -------------------------------------------------------------------- */
2764 : /* HDF-EOS Grid. */
2765 : /* -------------------------------------------------------------------- */
2766 : case H4ST_EOS_GRID:
2767 : {
2768 8 : int32 hGD, iProjCode = 0, iZoneCode = 0, iSphereCode = 0;
2769 : int32 nXSize, nYSize;
2770 : char szDimList[8192];
2771 : double adfUpLeft[2], adfLowRight[2], adfProjParms[15];
2772 :
2773 8 : if( poOpenInfo->eAccess == GA_ReadOnly )
2774 8 : poDS->hHDF4 = GDopen( poDS->pszFilename, DFACC_READ );
2775 : else
2776 0 : poDS->hHDF4 = GDopen( poDS->pszFilename, DFACC_WRITE );
2777 :
2778 8 : if( poDS->hHDF4 <= 0 )
2779 : {
2780 0 : delete poDS;
2781 0 : return( NULL );
2782 : }
2783 :
2784 8 : hGD = GDattach( poDS->hHDF4, poDS->pszSubdatasetName );
2785 :
2786 : /* -------------------------------------------------------------------- */
2787 : /* Decode the dimension map. */
2788 : /* -------------------------------------------------------------------- */
2789 : GDfieldinfo( hGD, poDS->pszFieldName, &poDS->iRank,
2790 8 : poDS->aiDimSizes, &poDS->iNumType, szDimList );
2791 : #ifdef DEBUG
2792 : CPLDebug( "HDF4Image",
2793 : "List of dimensions in grid %s: %s",
2794 8 : poDS->pszFieldName, szDimList);
2795 : #endif
2796 8 : poDS->GetImageDimensions( szDimList );
2797 :
2798 : int32 tilecode, tilerank;
2799 8 : if( GDtileinfo( hGD , poDS->pszFieldName , &tilecode, &tilerank, NULL ) == 0 )
2800 : {
2801 8 : if ( tilecode == HDFE_TILE )
2802 : {
2803 6 : int32 *tiledims = NULL;
2804 6 : tiledims = (int32 *) CPLCalloc( tilerank , sizeof( int32 ) );
2805 6 : GDtileinfo( hGD , poDS->pszFieldName , &tilecode, &tilerank, tiledims );
2806 12 : if ( ( tilerank == 2 ) && ( poDS->iRank == tilerank ) )
2807 : {
2808 6 : poDS->nBlockPreferredXSize = tiledims[1];
2809 6 : poDS->nBlockPreferredYSize = tiledims[0];
2810 6 : poDS->bReadTile = true;
2811 : #ifdef DEBUG
2812 : CPLDebug( "HDF4_EOS:EOS_GRID:","tilerank in grid %s: %d",
2813 6 : poDS->pszFieldName , (int)tilerank );
2814 : CPLDebug( "HDF4_EOS:EOS_GRID:","tiledimens in grid %s: %d,%d",
2815 6 : poDS->pszFieldName , (int)tiledims[0] , (int)tiledims[1] );
2816 : #endif
2817 : }
2818 : #ifdef DEBUG
2819 : else
2820 : {
2821 : CPLDebug( "HDF4_EOS:EOS_GRID:","tilerank in grid %s: %d not supported",
2822 0 : poDS->pszFieldName , (int)tilerank );
2823 : }
2824 : #endif
2825 6 : CPLFree(tiledims);
2826 : }
2827 : else
2828 : {
2829 : #ifdef DEBUG
2830 : CPLDebug( "HDF4_EOS:EOS_GRID:","tilecode == HDFE_NOTILE in grid %s: %d",
2831 : poDS->pszFieldName ,
2832 2 : (int)poDS->iRank );
2833 : #endif
2834 : }
2835 : }
2836 : #ifdef DEBUG
2837 : else
2838 : {
2839 0 : CPLDebug( "HDF4_EOS:EOS_GRID:","ERROR GDtileinfo %s", poDS->pszFieldName );
2840 : }
2841 : #endif
2842 :
2843 : /* -------------------------------------------------------------------- */
2844 : /* Fetch projection information */
2845 : /* -------------------------------------------------------------------- */
2846 8 : if ( GDprojinfo( hGD, &iProjCode, &iZoneCode,
2847 : &iSphereCode, adfProjParms) >= 0 )
2848 : {
2849 : #ifdef DEBUG
2850 : CPLDebug( "HDF4Image",
2851 : "Grid projection: "
2852 : "projection code: %ld, zone code %ld, "
2853 : "sphere code %ld",
2854 : (long)iProjCode, (long)iZoneCode,
2855 8 : (long)iSphereCode );
2856 : #endif
2857 : poDS->oSRS.importFromUSGS( iProjCode, iZoneCode,
2858 8 : adfProjParms, iSphereCode );
2859 :
2860 8 : if ( poDS->pszProjection )
2861 8 : CPLFree( poDS->pszProjection );
2862 8 : poDS->oSRS.exportToWkt( &poDS->pszProjection );
2863 : }
2864 :
2865 : /* -------------------------------------------------------------------- */
2866 : /* Fetch geotransformation matrix */
2867 : /* -------------------------------------------------------------------- */
2868 8 : if ( GDgridinfo( hGD, &nXSize, &nYSize,
2869 : adfUpLeft, adfLowRight ) >= 0 )
2870 : {
2871 : #ifdef DEBUG
2872 : CPLDebug( "HDF4Image",
2873 : "Grid geolocation: "
2874 : "top left X %f, top left Y %f, "
2875 : "low right X %f, low right Y %f, "
2876 : "cols %ld, rows %ld",
2877 : adfUpLeft[0], adfUpLeft[1],
2878 : adfLowRight[0], adfLowRight[1],
2879 8 : (long)nXSize, (long)nYSize );
2880 : #endif
2881 8 : if ( iProjCode )
2882 : {
2883 : // For projected systems coordinates are in meters
2884 : poDS->adfGeoTransform[1] =
2885 6 : (adfLowRight[0] - adfUpLeft[0]) / nXSize;
2886 : poDS->adfGeoTransform[5] =
2887 6 : (adfLowRight[1] - adfUpLeft[1]) / nYSize;
2888 6 : poDS->adfGeoTransform[0] = adfUpLeft[0];
2889 6 : poDS->adfGeoTransform[3] = adfUpLeft[1];
2890 : }
2891 : else
2892 : {
2893 : // Handle angular geographic coordinates here
2894 : poDS->adfGeoTransform[1] =
2895 : (CPLPackedDMSToDec(adfLowRight[0]) -
2896 2 : CPLPackedDMSToDec(adfUpLeft[0])) / nXSize;
2897 : poDS->adfGeoTransform[5] =
2898 : (CPLPackedDMSToDec(adfLowRight[1]) -
2899 2 : CPLPackedDMSToDec(adfUpLeft[1])) / nYSize;
2900 : poDS->adfGeoTransform[0] =
2901 2 : CPLPackedDMSToDec(adfUpLeft[0]);
2902 : poDS->adfGeoTransform[3] =
2903 2 : CPLPackedDMSToDec(adfUpLeft[1]);
2904 : }
2905 8 : poDS->adfGeoTransform[2] = 0.0;
2906 8 : poDS->adfGeoTransform[4] = 0.0;
2907 8 : poDS->bHasGeoTransform = TRUE;
2908 : }
2909 :
2910 : /* -------------------------------------------------------------------- */
2911 : /* Fetch metadata. */
2912 : /* -------------------------------------------------------------------- */
2913 8 : poDS->GetGridAttrs( hGD );
2914 :
2915 8 : GDdetach( hGD );
2916 :
2917 : /* -------------------------------------------------------------------- */
2918 : /* Fetch NODATA value. */
2919 : /* -------------------------------------------------------------------- */
2920 : pNoDataValue =
2921 8 : CPLMalloc( poDS->GetDataTypeSize(poDS->iNumType) );
2922 8 : if ( GDgetfillvalue( hGD, poDS->pszFieldName,
2923 : pNoDataValue ) != -1 )
2924 : {
2925 : dfNoData = poDS->AnyTypeToDouble( poDS->iNumType,
2926 0 : pNoDataValue );
2927 0 : bNoDataSet = TRUE;
2928 : }
2929 : else
2930 : {
2931 : const char *pszNoData =
2932 : CSLFetchNameValue( poDS->papszLocalMetadata,
2933 8 : "_FillValue" );
2934 8 : if ( pszNoData )
2935 : {
2936 6 : dfNoData = CPLAtof( pszNoData );
2937 6 : bNoDataSet = TRUE;
2938 : }
2939 : }
2940 8 : CPLFree( pNoDataValue );
2941 : }
2942 : break;
2943 :
2944 : default:
2945 : break;
2946 : }
2947 : }
2948 8 : break;
2949 :
2950 : /* -------------------------------------------------------------------- */
2951 : /* 'Plain' HDF scientific datasets. */
2952 : /* -------------------------------------------------------------------- */
2953 : case HDF4_SDS:
2954 : {
2955 : int32 iSDS;
2956 :
2957 255 : if( poOpenInfo->eAccess == GA_ReadOnly )
2958 255 : poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_READ, 0 );
2959 : else
2960 0 : poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_WRITE, 0 );
2961 :
2962 255 : if( poDS->hHDF4 <= 0 )
2963 : {
2964 0 : delete poDS;
2965 0 : return( NULL );
2966 : }
2967 :
2968 255 : poDS->hSD = SDstart( poDS->pszFilename, DFACC_READ );
2969 255 : if ( poDS->hSD == -1 )
2970 : {
2971 0 : delete poDS;
2972 0 : return NULL;
2973 : }
2974 :
2975 255 : if ( poDS->ReadGlobalAttributes( poDS->hSD ) != CE_None )
2976 : {
2977 0 : delete poDS;
2978 0 : return NULL;
2979 : }
2980 :
2981 : int32 nDatasets, nAttrs;
2982 :
2983 255 : if ( SDfileinfo( poDS->hSD, &nDatasets, &nAttrs ) != 0 )
2984 : {
2985 0 : delete poDS;
2986 0 : return NULL;
2987 : }
2988 :
2989 255 : if (poDS->iDataset < 0 || poDS->iDataset >= nDatasets)
2990 : {
2991 : CPLError(CE_Failure, CPLE_AppDefined,
2992 0 : "Subdataset index should be between 0 and %d", nDatasets - 1);
2993 0 : delete poDS;
2994 0 : return NULL;
2995 : }
2996 :
2997 255 : memset( poDS->aiDimSizes, 0, sizeof(int32) * H4_MAX_VAR_DIMS );
2998 255 : iSDS = SDselect( poDS->hSD, poDS->iDataset );
2999 : SDgetinfo( iSDS, poDS->szName, &poDS->iRank, poDS->aiDimSizes,
3000 255 : &poDS->iNumType, &poDS->nAttrs);
3001 :
3002 : // We will duplicate global metadata for every subdataset
3003 : poDS->papszLocalMetadata =
3004 255 : CSLDuplicate( poDS->papszGlobalMetadata );
3005 :
3006 286 : for ( iAttribute = 0; iAttribute < poDS->nAttrs; iAttribute++ )
3007 : {
3008 : char szAttrName[H4_MAX_NC_NAME];
3009 : SDattrinfo( iSDS, iAttribute, szAttrName,
3010 31 : &iAttrNumType, &nValues );
3011 : poDS->papszLocalMetadata =
3012 : poDS->TranslateHDF4Attributes( iSDS, iAttribute,
3013 : szAttrName, iAttrNumType,
3014 : nValues,
3015 31 : poDS->papszLocalMetadata );
3016 : }
3017 255 : poDS->SetMetadata( poDS->papszLocalMetadata, "" );
3018 255 : SDendaccess( iSDS );
3019 :
3020 : #ifdef DEBUG
3021 : CPLDebug( "HDF4Image",
3022 : "aiDimSizes[0]=%ld, aiDimSizes[1]=%ld, "
3023 : "aiDimSizes[2]=%ld, aiDimSizes[3]=%ld",
3024 : (long)poDS->aiDimSizes[0], (long)poDS->aiDimSizes[1],
3025 255 : (long)poDS->aiDimSizes[2], (long)poDS->aiDimSizes[3] );
3026 : #endif
3027 :
3028 255 : switch( poDS->iRank )
3029 : {
3030 : case 1:
3031 0 : poDS->nBandCount = 1;
3032 0 : poDS->iXDim = 0;
3033 0 : poDS->iYDim = -1;
3034 0 : break;
3035 : case 2:
3036 109 : poDS->nBandCount = 1;
3037 109 : poDS->iXDim = 1;
3038 109 : poDS->iYDim = 0;
3039 109 : break;
3040 : case 3:
3041 : /* FIXME ? We should probably remove the following test as there are valid datasets */
3042 : /* where the height is lower than the band number : for example
3043 : http://www.iapmw.unibe.ch/research/projects/FriOWL/data/otd/LISOTD_HRAC_V2.2.hdf */
3044 : /* which is a 720x360 x 365 bands */
3045 : /* Use a HACK for now */
3046 146 : if( poDS->aiDimSizes[0] < poDS->aiDimSizes[2] &&
3047 : !(poDS->aiDimSizes[0] == 360 &&
3048 : poDS->aiDimSizes[1] == 720 &&
3049 : poDS->aiDimSizes[2] == 365) )
3050 : {
3051 0 : poDS->iBandDim = 0;
3052 0 : poDS->iXDim = 2;
3053 0 : poDS->iYDim = 1;
3054 : }
3055 : else
3056 : {
3057 146 : if( poDS->aiDimSizes[1] <= poDS->aiDimSizes[0] &&
3058 : poDS->aiDimSizes[1] <= poDS->aiDimSizes[2] )
3059 : {
3060 0 : poDS->iBandDim = 1;
3061 0 : poDS->iXDim = 2;
3062 0 : poDS->iYDim = 0;
3063 : }
3064 : else
3065 : {
3066 146 : poDS->iBandDim = 2;
3067 146 : poDS->iXDim = 1;
3068 146 : poDS->iYDim = 0;
3069 :
3070 : }
3071 : }
3072 146 : poDS->nBandCount = poDS->aiDimSizes[poDS->iBandDim];
3073 146 : break;
3074 : case 4: // FIXME
3075 0 : poDS->nBandCount = poDS->aiDimSizes[2] * poDS->aiDimSizes[3];
3076 : break;
3077 : default:
3078 : break;
3079 : }
3080 :
3081 : // We preset this because CaptureNRLGeoTransform needs it.
3082 255 : poDS->nRasterXSize = poDS->aiDimSizes[poDS->iXDim];
3083 255 : if (poDS->iYDim >= 0)
3084 255 : poDS->nRasterYSize = poDS->aiDimSizes[poDS->iYDim];
3085 : else
3086 0 : poDS->nRasterYSize = 1;
3087 :
3088 : // Special case projection info for NRL generated files.
3089 : const char *pszMapProjectionSystem =
3090 : CSLFetchNameValue(poDS->papszGlobalMetadata,
3091 255 : "mapProjectionSystem");
3092 255 : if( pszMapProjectionSystem != NULL
3093 : && EQUAL(pszMapProjectionSystem,"NRL(USGS)") )
3094 : {
3095 0 : poDS->CaptureNRLGeoTransform();
3096 : }
3097 :
3098 : // Special case for coastwatch hdf files.
3099 255 : if( CSLFetchNameValue( poDS->papszGlobalMetadata,
3100 : "gctp_sys" ) != NULL )
3101 0 : poDS->CaptureCoastwatchGCTPInfo();
3102 :
3103 : // Special case for MODIS geolocation
3104 255 : poDS->ProcessModisSDSGeolocation();
3105 :
3106 : // Special case for NASA/CCRS Landsat in HDF
3107 255 : poDS->CaptureL1GMTLInfo();
3108 : }
3109 255 : break;
3110 :
3111 : /* -------------------------------------------------------------------- */
3112 : /* 'Plain' HDF rasters. */
3113 : /* -------------------------------------------------------------------- */
3114 : case HDF4_GR:
3115 0 : if( poOpenInfo->eAccess == GA_ReadOnly )
3116 0 : poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_READ, 0 );
3117 : else
3118 0 : poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_WRITE, 0 );
3119 :
3120 0 : if( poDS->hHDF4 <= 0 )
3121 : {
3122 0 : delete poDS;
3123 0 : return( NULL );
3124 : }
3125 :
3126 0 : poDS->hGR = GRstart( poDS->hHDF4 );
3127 0 : if ( poDS->hGR == -1 )
3128 : {
3129 0 : delete poDS;
3130 0 : return NULL;
3131 : }
3132 :
3133 0 : poDS->iGR = GRselect( poDS->hGR, poDS->iDataset );
3134 0 : if ( GRgetiminfo( poDS->iGR, poDS->szName,
3135 : &poDS->iRank, &poDS->iNumType,
3136 : &poDS->iInterlaceMode, poDS->aiDimSizes,
3137 : &poDS->nAttrs ) != 0 )
3138 : {
3139 0 : delete poDS;
3140 0 : return NULL;
3141 : }
3142 :
3143 : // We will duplicate global metadata for every subdataset
3144 0 : poDS->papszLocalMetadata = CSLDuplicate( poDS->papszGlobalMetadata );
3145 :
3146 0 : for ( iAttribute = 0; iAttribute < poDS->nAttrs; iAttribute++ )
3147 : {
3148 : char szAttrName[H4_MAX_NC_NAME];
3149 : GRattrinfo( poDS->iGR, iAttribute, szAttrName,
3150 0 : &iAttrNumType, &nValues );
3151 : poDS->papszLocalMetadata =
3152 : poDS->TranslateHDF4Attributes( poDS->iGR, iAttribute,
3153 : szAttrName, iAttrNumType,
3154 : nValues,
3155 0 : poDS->papszLocalMetadata );
3156 : }
3157 0 : poDS->SetMetadata( poDS->papszLocalMetadata, "" );
3158 : // Read colour table
3159 : GDALColorEntry oEntry;
3160 :
3161 0 : poDS->iPal = GRgetlutid ( poDS->iGR, poDS->iDataset );
3162 0 : if ( poDS->iPal != -1 )
3163 : {
3164 : GRgetlutinfo( poDS->iPal, &poDS->nComps, &poDS->iPalDataType,
3165 0 : &poDS->iPalInterlaceMode, &poDS->nPalEntries );
3166 0 : GRreadlut( poDS->iPal, poDS->aiPaletteData );
3167 0 : poDS->poColorTable = new GDALColorTable();
3168 0 : for( i = 0; i < N_COLOR_ENTRIES; i++ )
3169 : {
3170 0 : oEntry.c1 = poDS->aiPaletteData[i][0];
3171 0 : oEntry.c2 = poDS->aiPaletteData[i][1];
3172 0 : oEntry.c3 = poDS->aiPaletteData[i][2];
3173 0 : oEntry.c4 = 255;
3174 :
3175 0 : poDS->poColorTable->SetColorEntry( i, &oEntry );
3176 : }
3177 : }
3178 :
3179 0 : poDS->iXDim = 0;
3180 0 : poDS->iYDim = 1;
3181 0 : poDS->nBandCount = poDS->iRank;
3182 0 : break;
3183 : default:
3184 0 : delete poDS;
3185 0 : return NULL;
3186 : }
3187 :
3188 263 : poDS->nRasterXSize = poDS->aiDimSizes[poDS->iXDim];
3189 263 : if (poDS->iYDim >= 0)
3190 263 : poDS->nRasterYSize = poDS->aiDimSizes[poDS->iYDim];
3191 : else
3192 0 : poDS->nRasterYSize = 1;
3193 :
3194 263 : if ( poDS->iSubdatasetType == H4ST_HYPERION_L1 )
3195 : {
3196 : // XXX: Hyperion SDSs has Height x Bands x Width dimensions scheme
3197 0 : if ( poDS->iRank > 2 )
3198 : {
3199 0 : poDS->nBandCount = poDS->aiDimSizes[1];
3200 0 : poDS->nRasterXSize = poDS->aiDimSizes[2];
3201 0 : poDS->nRasterYSize = poDS->aiDimSizes[0];
3202 : }
3203 : else
3204 : {
3205 0 : poDS->nBandCount = poDS->aiDimSizes[0];
3206 0 : poDS->nRasterXSize = poDS->aiDimSizes[1];
3207 0 : poDS->nRasterYSize = 1;
3208 : }
3209 : }
3210 :
3211 : /* -------------------------------------------------------------------- */
3212 : /* Create band information objects. */
3213 : /* -------------------------------------------------------------------- */
3214 602 : for( i = 1; i <= poDS->nBandCount; i++ )
3215 : {
3216 : HDF4ImageRasterBand *poBand =
3217 : new HDF4ImageRasterBand( poDS, i,
3218 339 : poDS->GetDataType( poDS->iNumType ) );
3219 339 : poDS->SetBand( i, poBand );
3220 :
3221 339 : if ( bNoDataSet )
3222 6 : poBand->SetNoDataValue( dfNoData );
3223 : }
3224 :
3225 : /* -------------------------------------------------------------------- */
3226 : /* Now we will handle particular types of HDF products. Every */
3227 : /* HDF product has its own structure. */
3228 : /* -------------------------------------------------------------------- */
3229 :
3230 : // Variables for reading georeferencing
3231 : double dfULX, dfULY, dfLRX, dfLRY;
3232 :
3233 263 : switch ( poDS->iSubdatasetType )
3234 : {
3235 : /* -------------------------------------------------------------------- */
3236 : /* HDF, created by GDAL. */
3237 : /* -------------------------------------------------------------------- */
3238 : case H4ST_GDAL:
3239 : {
3240 : const char *pszValue;
3241 :
3242 : CPLDebug( "HDF4Image",
3243 249 : "Input dataset interpreted as GDAL_HDF4" );
3244 :
3245 249 : if ( (pszValue =
3246 : CSLFetchNameValue(poDS->papszGlobalMetadata,
3247 : "Projection")) )
3248 : {
3249 105 : if ( poDS->pszProjection )
3250 105 : CPLFree( poDS->pszProjection );
3251 105 : poDS->pszProjection = CPLStrdup( pszValue );
3252 : }
3253 249 : if ( (pszValue = CSLFetchNameValue(poDS->papszGlobalMetadata,
3254 : "TransformationMatrix")) )
3255 : {
3256 249 : int i = 0;
3257 249 : char *pszString = (char *) pszValue;
3258 1992 : while ( *pszValue && i < 6 )
3259 : {
3260 1494 : poDS->adfGeoTransform[i++] = strtod(pszString, &pszString);
3261 1494 : pszString++;
3262 : }
3263 249 : poDS->bHasGeoTransform = TRUE;
3264 : }
3265 560 : for( i = 1; i <= poDS->nBands; i++ )
3266 : {
3267 311 : if ( (pszValue =
3268 : CSLFetchNameValue(poDS->papszGlobalMetadata,
3269 : CPLSPrintf("BandDesc%d", i))) )
3270 32 : poDS->GetRasterBand( i )->SetDescription( pszValue );
3271 : }
3272 560 : for( i = 1; i <= poDS->nBands; i++ )
3273 : {
3274 311 : if ( (pszValue =
3275 : CSLFetchNameValue(poDS->papszGlobalMetadata,
3276 : CPLSPrintf("NoDataValue%d", i))) )
3277 32 : poDS->GetRasterBand(i)->SetNoDataValue( CPLAtof(pszValue) );
3278 : }
3279 : }
3280 249 : break;
3281 :
3282 : /* -------------------------------------------------------------------- */
3283 : /* SeaWiFS Level 3 Standard Mapped Image Products. */
3284 : /* Organized similar to MODIS Level 3 products. */
3285 : /* -------------------------------------------------------------------- */
3286 : case H4ST_SEAWIFS_L3:
3287 : {
3288 0 : CPLDebug( "HDF4Image", "Input dataset interpreted as SEAWIFS_L3" );
3289 :
3290 : // Read band description
3291 0 : for ( i = 1; i <= poDS->nBands; i++ )
3292 : {
3293 : poDS->GetRasterBand( i )->SetDescription(
3294 0 : CSLFetchNameValue( poDS->papszGlobalMetadata, "Parameter" ) );
3295 : }
3296 :
3297 : // Read coordinate system and geotransform matrix
3298 0 : poDS->oSRS.SetWellKnownGeogCS( "WGS84" );
3299 :
3300 0 : if ( EQUAL(CSLFetchNameValue(poDS->papszGlobalMetadata,
3301 : "Map Projection"),
3302 : "Equidistant Cylindrical") )
3303 : {
3304 0 : poDS->oSRS.SetEquirectangular( 0.0, 0.0, 0.0, 0.0 );
3305 0 : poDS->oSRS.SetLinearUnits( SRS_UL_METER, 1 );
3306 0 : if ( poDS->pszProjection )
3307 0 : CPLFree( poDS->pszProjection );
3308 0 : poDS->oSRS.exportToWkt( &poDS->pszProjection );
3309 : }
3310 :
3311 : dfULX = CPLAtof( CSLFetchNameValue(poDS->papszGlobalMetadata,
3312 0 : "Westernmost Longitude") );
3313 : dfULY = CPLAtof( CSLFetchNameValue(poDS->papszGlobalMetadata,
3314 0 : "Northernmost Latitude") );
3315 : dfLRX = CPLAtof( CSLFetchNameValue(poDS->papszGlobalMetadata,
3316 0 : "Easternmost Longitude") );
3317 : dfLRY = CPLAtof( CSLFetchNameValue(poDS->papszGlobalMetadata,
3318 0 : "Southernmost Latitude") );
3319 0 : poDS->ToGeoref( &dfULX, &dfULY );
3320 0 : poDS->ToGeoref( &dfLRX, &dfLRY );
3321 0 : poDS->adfGeoTransform[0] = dfULX;
3322 0 : poDS->adfGeoTransform[3] = dfULY;
3323 0 : poDS->adfGeoTransform[1] = (dfLRX - dfULX) / poDS->nRasterXSize;
3324 0 : poDS->adfGeoTransform[5] = (dfULY - dfLRY) / poDS->nRasterYSize;
3325 0 : if ( dfULY > 0) // Northern hemisphere
3326 0 : poDS->adfGeoTransform[5] = - poDS->adfGeoTransform[5];
3327 0 : poDS->adfGeoTransform[2] = 0.0;
3328 0 : poDS->adfGeoTransform[4] = 0.0;
3329 0 : poDS->bHasGeoTransform = TRUE;
3330 : }
3331 0 : break;
3332 :
3333 :
3334 : /* -------------------------------------------------------------------- */
3335 : /* Generic SDS */
3336 : /* -------------------------------------------------------------------- */
3337 : case H4ST_UNKNOWN:
3338 : {
3339 :
3340 : // This is a coastwatch convention.
3341 6 : if( CSLFetchNameValue( poDS->papszLocalMetadata, "missing_value" ) )
3342 : {
3343 : int i;
3344 0 : for( i = 1; i <= poDS->nBands; i++ )
3345 : {
3346 : poDS->GetRasterBand(i)->SetNoDataValue(
3347 : CPLAtof( CSLFetchNameValue(poDS->papszLocalMetadata,
3348 0 : "missing_value") ) );
3349 : }
3350 : }
3351 :
3352 : // Coastwatch offset and scale.
3353 6 : if( CSLFetchNameValue( poDS->papszLocalMetadata, "scale_factor" )
3354 : && CSLFetchNameValue( poDS->papszLocalMetadata, "add_offset" ) )
3355 : {
3356 0 : for( i = 1; i <= poDS->nBands; i++ )
3357 : {
3358 : HDF4ImageRasterBand *poBand =
3359 0 : (HDF4ImageRasterBand *) poDS->GetRasterBand(i);
3360 :
3361 0 : poBand->bHaveScaleAndOffset = TRUE;
3362 : poBand->dfScale =
3363 : CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3364 0 : "scale_factor" ) );
3365 : poBand->dfOffset = -1 * poBand->dfScale *
3366 : CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3367 0 : "add_offset" ) );
3368 : }
3369 : }
3370 :
3371 : // this is a modis level3 convention (data from ACT)
3372 : // Eg data/hdf/act/modis/MODAM2004280160000.L3_NOAA_GMX
3373 :
3374 6 : if( CSLFetchNameValue( poDS->papszLocalMetadata,
3375 : "scalingSlope" )
3376 : && CSLFetchNameValue( poDS->papszLocalMetadata,
3377 : "scalingIntercept" ) )
3378 : {
3379 : int i;
3380 0 : CPLString osUnits;
3381 :
3382 0 : if( CSLFetchNameValue( poDS->papszLocalMetadata,
3383 : "productUnits" ) )
3384 : {
3385 : osUnits = CSLFetchNameValue( poDS->papszLocalMetadata,
3386 0 : "productUnits" );
3387 : }
3388 :
3389 0 : for( i = 1; i <= poDS->nBands; i++ )
3390 : {
3391 : HDF4ImageRasterBand *poBand =
3392 0 : (HDF4ImageRasterBand *) poDS->GetRasterBand(i);
3393 :
3394 0 : poBand->bHaveScaleAndOffset = TRUE;
3395 : poBand->dfScale =
3396 : CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3397 0 : "scalingSlope" ) );
3398 : poBand->dfOffset =
3399 : CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3400 0 : "scalingIntercept" ) );
3401 :
3402 0 : poBand->osUnitType = osUnits;
3403 0 : }
3404 : }
3405 : }
3406 6 : break;
3407 :
3408 : /* -------------------------------------------------------------------- */
3409 : /* Hyperion Level 1. */
3410 : /* -------------------------------------------------------------------- */
3411 : case H4ST_HYPERION_L1:
3412 : {
3413 0 : CPLDebug( "HDF4Image", "Input dataset interpreted as HYPERION_L1" );
3414 : }
3415 : break;
3416 :
3417 : default:
3418 : break;
3419 : }
3420 :
3421 : /* -------------------------------------------------------------------- */
3422 : /* Setup PAM info for this subdataset */
3423 : /* -------------------------------------------------------------------- */
3424 263 : poDS->SetPhysicalFilename( poDS->pszFilename );
3425 263 : poDS->SetSubdatasetName( osSubdatasetName );
3426 :
3427 263 : poDS->TryLoadXML();
3428 :
3429 263 : poDS->oOvManager.Initialize( poDS, ":::VIRTUAL:::" );
3430 :
3431 263 : return( poDS );
3432 : }
3433 :
3434 : /************************************************************************/
3435 : /* Create() */
3436 : /************************************************************************/
3437 :
3438 : GDALDataset *HDF4ImageDataset::Create( const char * pszFilename,
3439 : int nXSize, int nYSize, int nBands,
3440 : GDALDataType eType,
3441 145 : char **papszOptions )
3442 :
3443 : {
3444 : /* -------------------------------------------------------------------- */
3445 : /* Create the dataset. */
3446 : /* -------------------------------------------------------------------- */
3447 : HDF4ImageDataset *poDS;
3448 : const char *pszSDSName;
3449 : int iBand;
3450 145 : int32 iSDS = -1;
3451 : int32 aiDimSizes[H4_MAX_VAR_DIMS];
3452 :
3453 145 : if( nBands == 0 )
3454 : {
3455 : CPLError( CE_Failure, CPLE_NotSupported,
3456 1 : "Unable to export files with zero bands." );
3457 1 : return NULL;
3458 : }
3459 :
3460 144 : poDS = new HDF4ImageDataset();
3461 :
3462 : /* -------------------------------------------------------------------- */
3463 : /* Choose rank for the created dataset. */
3464 : /* -------------------------------------------------------------------- */
3465 144 : poDS->iRank = 3;
3466 248 : if ( CSLFetchNameValue( papszOptions, "RANK" ) != NULL &&
3467 : EQUAL( CSLFetchNameValue( papszOptions, "RANK" ), "2" ) )
3468 48 : poDS->iRank = 2;
3469 :
3470 144 : poDS->hSD = SDstart( pszFilename, DFACC_CREATE );
3471 144 : if ( poDS->hSD == -1 )
3472 : {
3473 : CPLError( CE_Failure, CPLE_OpenFailed,
3474 0 : "Can't create HDF4 file %s", pszFilename );
3475 0 : return NULL;
3476 : }
3477 144 : poDS->iXDim = 1;
3478 144 : poDS->iYDim = 0;
3479 144 : poDS->iBandDim = 2;
3480 144 : aiDimSizes[poDS->iXDim] = nXSize;
3481 144 : aiDimSizes[poDS->iYDim] = nYSize;
3482 144 : aiDimSizes[poDS->iBandDim] = nBands;
3483 :
3484 144 : if ( poDS->iRank == 2 )
3485 : {
3486 96 : for ( iBand = 0; iBand < nBands; iBand++ )
3487 : {
3488 48 : pszSDSName = CPLSPrintf( "Band%d", iBand );
3489 48 : switch ( eType )
3490 : {
3491 : case GDT_Float64:
3492 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_FLOAT64,
3493 6 : poDS->iRank, aiDimSizes );
3494 6 : break;
3495 : case GDT_Float32:
3496 : iSDS = SDcreate( poDS-> hSD, pszSDSName, DFNT_FLOAT32,
3497 6 : poDS->iRank, aiDimSizes );
3498 6 : break;
3499 : case GDT_UInt32:
3500 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT32,
3501 6 : poDS->iRank, aiDimSizes );
3502 6 : break;
3503 : case GDT_UInt16:
3504 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT16,
3505 6 : poDS->iRank, aiDimSizes );
3506 6 : break;
3507 : case GDT_Int32:
3508 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT32,
3509 6 : poDS->iRank, aiDimSizes );
3510 6 : break;
3511 : case GDT_Int16:
3512 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT16,
3513 6 : poDS->iRank, aiDimSizes );
3514 6 : break;
3515 : case GDT_Byte:
3516 : default:
3517 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT8,
3518 12 : poDS->iRank, aiDimSizes );
3519 : break;
3520 : }
3521 48 : SDendaccess( iSDS );
3522 : }
3523 : }
3524 96 : else if ( poDS->iRank == 3 )
3525 : {
3526 96 : pszSDSName = "3-dimensional Scientific Dataset";
3527 96 : poDS->iDataset = 0;
3528 96 : switch ( eType )
3529 : {
3530 : case GDT_Float64:
3531 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_FLOAT64,
3532 10 : poDS->iRank, aiDimSizes );
3533 10 : break;
3534 : case GDT_Float32:
3535 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_FLOAT32,
3536 10 : poDS->iRank, aiDimSizes );
3537 10 : break;
3538 : case GDT_UInt32:
3539 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT32,
3540 10 : poDS->iRank, aiDimSizes );
3541 10 : break;
3542 : case GDT_UInt16:
3543 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT16,
3544 10 : poDS->iRank, aiDimSizes );
3545 10 : break;
3546 : case GDT_Int32:
3547 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT32,
3548 10 : poDS->iRank, aiDimSizes );
3549 10 : break;
3550 : case GDT_Int16:
3551 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_INT16,
3552 10 : poDS->iRank, aiDimSizes );
3553 10 : break;
3554 : case GDT_Byte:
3555 : default:
3556 : iSDS = SDcreate( poDS->hSD, pszSDSName, DFNT_UINT8,
3557 36 : poDS->iRank, aiDimSizes );
3558 : break;
3559 : }
3560 : }
3561 : else // Should never happen
3562 0 : return NULL;
3563 :
3564 144 : if ( iSDS < 0 )
3565 : {
3566 : CPLError( CE_Failure, CPLE_AppDefined,
3567 : "Can't create SDS with rank %ld for file %s",
3568 0 : (long)poDS->iRank, pszFilename );
3569 0 : return NULL;
3570 : }
3571 :
3572 144 : poDS->nRasterXSize = nXSize;
3573 144 : poDS->nRasterYSize = nYSize;
3574 144 : poDS->eAccess = GA_Update;
3575 144 : poDS->iDatasetType = HDF4_SDS;
3576 144 : poDS->iSubdatasetType = H4ST_GDAL;
3577 144 : poDS->nBands = nBands;
3578 :
3579 : /* -------------------------------------------------------------------- */
3580 : /* Create band information objects. */
3581 : /* -------------------------------------------------------------------- */
3582 688 : for( iBand = 1; iBand <= nBands; iBand++ )
3583 200 : poDS->SetBand( iBand, new HDF4ImageRasterBand( poDS, iBand, eType ) );
3584 :
3585 : SDsetattr( poDS->hSD, "Signature", DFNT_CHAR8, strlen(pszGDALSignature) + 1,
3586 144 : pszGDALSignature );
3587 :
3588 144 : return (GDALDataset *) poDS;
3589 : }
3590 :
3591 : /************************************************************************/
3592 : /* GDALRegister_HDF4Image() */
3593 : /************************************************************************/
3594 :
3595 409 : void GDALRegister_HDF4Image()
3596 :
3597 : {
3598 : GDALDriver *poDriver;
3599 :
3600 409 : if( GDALGetDriverByName( "HDF4Image" ) == NULL )
3601 : {
3602 392 : poDriver = new GDALDriver();
3603 :
3604 392 : poDriver->SetDescription( "HDF4Image" );
3605 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
3606 392 : "HDF4 Dataset" );
3607 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
3608 392 : "frmt_hdf4.html" );
3609 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
3610 392 : "Byte Int16 UInt16 Int32 UInt32 Float32 Float64" );
3611 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
3612 : "<CreationOptionList>"
3613 : " <Option name='RANK' type='int' description='Rank of output SDS'/>"
3614 392 : "</CreationOptionList>" );
3615 :
3616 392 : poDriver->pfnOpen = HDF4ImageDataset::Open;
3617 392 : poDriver->pfnCreate = HDF4ImageDataset::Create;
3618 :
3619 392 : GetGDALDriverManager()->RegisterDriver( poDriver );
3620 : }
3621 409 : }
3622 :
|
