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