← Back to branch summary

~ubuntu-branches/debian/sid/gdal/sid

« back to all changes in this revision

Viewing changes to frmts/hdf4/hdf4imagedataset.cpp

  • Committer: Package Import Robot
  • Author(s): Francesco Paolo Lovergine
  • Date: 2012-05-07 15:04:42 UTC
  • mfrom: (5.5.16 experimental)
  • Revision ID: package-import@ubuntu.com-20120507150442-2eks97loeh6rq005
Tags: 1.9.0-1
* Ready for sid, starting transition.
* All symfiles updated to latest builds.
* Added dh_numpy call in debian/rules to depend on numpy ABI.
* Policy bumped to 3.9.3, no changes required.

Show diffs side-by-side

added added

removed removed

Lines of Context:
frmts/hdf4/hdf4imagedataset.cpp
1
1
/******************************************************************************
2
 
 * $Id: hdf4imagedataset.cpp 18946 2010-02-27 20:39:51Z rouault $
 
2
 * $Id: hdf4imagedataset.cpp 22765 2011-07-23 09:56:53Z rouault $
3
3
 *
4
4
 * Project:  Hierarchical Data Format Release 4 (HDF4)
5
5
 * Purpose:  Read subdatasets of HDF4 file.
43
43
#include "hdf4compat.h"
44
44
#include "hdf4dataset.h"
45
45
 
46
 
CPL_CVSID("$Id: hdf4imagedataset.cpp 18946 2010-02-27 20:39:51Z rouault $");
 
46
#include "nasakeywordhandler.h"
 
47
 
 
48
CPL_CVSID("$Id: hdf4imagedataset.cpp 22765 2011-07-23 09:56:53Z rouault $");
47
49
 
48
50
CPL_C_START
49
51
void    GDALRegister_HDF4(void);
69
71
    PROD_ASTER_L2,
70
72
    PROD_ASTER_L3,
71
73
    PROD_AST14DEM,
72
 
    PROD_MODIS_L1B
 
74
    PROD_MODIS_L1B,
 
75
    PROD_MODIS_L2
73
76
};
74
77
 
75
78
/************************************************************************/
109
112
 
110
113
    HDF4DatasetType iDatasetType;
111
114
 
 
115
    int32       iSDS;
 
116
 
 
117
    int         nBlockPreferredXSize;
 
118
    int         nBlockPreferredYSize;
 
119
    bool        bReadTile;
 
120
 
112
121
    void                ToGeoref( double *, double * );
113
122
    void                GetImageDimensions( char * );
114
123
    void                GetSwatAttrs( int32 );
115
124
    void                GetGridAttrs( int32 hGD );
116
125
    void                CaptureNRLGeoTransform(void);
 
126
    void                CaptureL1GMTLInfo(void);
117
127
    void                CaptureCoastwatchGCTPInfo(void);
118
128
    void                ProcessModisSDSGeolocation(void);
119
129
    int                 ProcessSwathGeolocation( int32, char ** );
124
134
  public:
125
135
                HDF4ImageDataset();
126
136
                ~HDF4ImageDataset();
127
 
    
 
137
 
128
138
    static GDALDataset  *Open( GDALOpenInfo * );
129
139
    static GDALDataset  *Create( const char * pszFilename,
130
140
                                 int nXSize, int nYSize, int nBands,
152
162
    int         bNoDataSet;
153
163
    double      dfNoDataValue;
154
164
 
155
 
    int         bHaveScaleAndOffset;
 
165
    int         bHaveScale, bHaveOffset;
156
166
    double      dfScale;
157
167
    double      dfOffset;
158
168
 
161
171
  public:
162
172
 
163
173
                HDF4ImageRasterBand( HDF4ImageDataset *, int, GDALDataType );
164
 
    
 
174
 
165
175
    virtual CPLErr          IReadBlock( int, int, void * );
166
176
    virtual CPLErr          IWriteBlock( int, int, void * );
167
177
    virtual GDALColorInterp GetColorInterpretation();
187
197
    bNoDataSet = FALSE;
188
198
    dfNoDataValue = -9999.0;
189
199
 
190
 
    bHaveScaleAndOffset = FALSE;
 
200
    bHaveScale = bHaveOffset = FALSE;
191
201
    dfScale = 1.0;
192
202
    dfOffset = 0.0;
193
203
 
194
204
    nBlockXSize = poDS->GetRasterXSize();
195
205
 
196
 
    // Aim for a block of about 100000 pixels.  Chunking up substantially
 
206
    // Aim for a block of about 1000000 pixels.  Chunking up substantially
197
207
    // improves performance in some situations.  For now we only chunk up for
198
 
    // SDS based datasets since other variations haven't been tested. #2208  
199
 
    if( poDS->iDatasetType == HDF4_SDS )
 
208
    // SDS and EOS based datasets since other variations haven't been tested. #2208
 
209
    if( poDS->iDatasetType == HDF4_SDS ||
 
210
        poDS->iDatasetType == HDF4_EOS)
200
211
    {
201
 
        int nChunkSize = 
202
 
            atoi( CPLGetConfigOption("HDF4_BLOCK_PIXELS", "100000") );
 
212
        int nChunkSize =
 
213
            atoi( CPLGetConfigOption("HDF4_BLOCK_PIXELS", "1000000") );
203
214
 
204
215
        nBlockYSize = nChunkSize / poDS->GetRasterXSize();
205
216
        nBlockYSize = MAX(1,MIN(nBlockYSize,poDS->GetRasterYSize()));
206
217
    }
207
218
    else
 
219
    {
208
220
        nBlockYSize = 1;
 
221
    }
 
222
 
 
223
    /* HDF4_EOS:EOS_GRID case. We ensure that the block size matches */
 
224
    /* the raster width, as the IReadBlock() code can only handle multiple */
 
225
    /* blocks per row */
 
226
    if ( poDS->nBlockPreferredXSize == nBlockXSize &&
 
227
         poDS->nBlockPreferredYSize > 0 )
 
228
    {
 
229
        if (poDS->nBlockPreferredYSize == 1)
 
230
        {
 
231
            /* Avoid defaulting to tile reading when the preferred height is 1 */
 
232
            /* as it leads to very poor performance with : */
 
233
            /* ftp://e4ftl01u.ecs.nasa.gov/MODIS_Composites/MOLT/MOD13Q1.005/2006.06.10/MOD13Q1.A2006161.h21v13.005.2008234103220.hd */
 
234
            poDS->bReadTile = FALSE;
 
235
        }
 
236
        else
 
237
        {
 
238
            nBlockYSize = poDS->nBlockPreferredYSize;
 
239
        }
 
240
    }
209
241
}
210
242
 
211
243
/************************************************************************/
232
264
    int nYOff = nBlockYOff * nBlockYSize;
233
265
    int nYSize = MIN(nYOff + nBlockYSize, poDS->GetRasterYSize()) - nYOff;
234
266
    CPLAssert( nBlockXOff == 0 );
235
 
    
 
267
 
236
268
/* -------------------------------------------------------------------- */
237
269
/*      HDF files with external data files, such as some landsat        */
238
270
/*      products (eg. data/hdf/L1G) need to be told what directory      */
248
280
    {
249
281
      case HDF4_SDS:
250
282
      {
251
 
          int32   iSDS = SDselect( poGDS->hSD, poGDS->iDataset );
 
283
          /* We avoid doing SDselect() / SDendaccess() for each block access */
 
284
          /* as this is very slow when zlib compression is used */
 
285
 
 
286
          if (poGDS->iSDS == FAIL)
 
287
              poGDS->iSDS = SDselect( poGDS->hSD, poGDS->iDataset );
 
288
 
252
289
          /* HDF rank:
253
 
             A rank 2 dataset is an image read in scan-line order (2D). 
 
290
             A rank 2 dataset is an image read in scan-line order (2D).
254
291
             A rank 3 dataset is a series of images which are read in
255
292
             an image at a time to form a volume.
256
293
             A rank 4 dataset may be thought of as a series of volumes.
285
322
            case 3: // 3Dim: volume
286
323
              aiStart[poGDS->iBandDim] = nBand - 1;
287
324
              aiEdges[poGDS->iBandDim] = 1;
288
 
                    
 
325
 
289
326
              aiStart[poGDS->iYDim] = nYOff;
290
327
              aiEdges[poGDS->iYDim] = nYSize;
291
 
                    
 
328
 
292
329
              aiStart[poGDS->iXDim] = nBlockXOff;
293
330
              aiEdges[poGDS->iXDim] = nBlockXSize;
294
331
              break;
295
332
            case 2: // 2Dim: rows/cols
296
333
              aiStart[poGDS->iYDim] = nYOff;
297
334
              aiEdges[poGDS->iYDim] = nYSize;
298
 
                    
 
335
 
299
336
              aiStart[poGDS->iXDim] = nBlockXOff;
300
337
              aiEdges[poGDS->iXDim] = nBlockXSize;
301
338
              break;
304
341
              aiEdges[poGDS->iXDim] = nBlockXSize;
305
342
              break;
306
343
          }
307
 
                
 
344
 
308
345
          // Read HDF SDS array
309
 
          if( SDreaddata( iSDS, aiStart, NULL, aiEdges, pImage ) < 0 )
 
346
          if( SDreaddata( poGDS->iSDS, aiStart, NULL, aiEdges, pImage ) < 0 )
310
347
          {
311
 
              CPLError( CE_Failure, CPLE_AppDefined, 
 
348
              CPLError( CE_Failure, CPLE_AppDefined,
312
349
                        "SDreaddata() failed for block." );
313
350
              eErr = CE_Failure;
314
351
          }
315
 
                
316
 
          SDendaccess( iSDS );
 
352
 
 
353
          //SDendaccess( iSDS );
317
354
      }
318
355
      break;
319
356
 
324
361
          GByte    *pbBuffer = (GByte *)
325
362
              CPLMalloc(nBlockXSize*nBlockYSize*poGDS->iRank*nBlockYSize);
326
363
          int     i, j;
327
 
            
 
364
 
328
365
          aiStart[poGDS->iYDim] = nYOff;
329
366
          aiEdges[poGDS->iYDim] = nYSize;
330
 
            
 
367
 
331
368
          aiStart[poGDS->iXDim] = nBlockXOff;
332
369
          aiEdges[poGDS->iXDim] = nBlockXSize;
333
370
 
334
371
          if( GRreadimage(poGDS->iGR, aiStart, NULL, aiEdges, pbBuffer) < 0 )
335
372
          {
336
 
              CPLError( CE_Failure, CPLE_AppDefined, 
 
373
              CPLError( CE_Failure, CPLE_AppDefined,
337
374
                        "GRreaddata() failed for block." );
338
375
              eErr = CE_Failure;
339
376
          }
371
408
                        (nBand - 1)
372
409
                        % poGDS->aiDimSizes[poGDS->iBandDim];
373
410
                    aiEdges[poGDS->iBandDim] = 1;
374
 
                                
 
411
 
375
412
                    aiStart[poGDS->iYDim] = nYOff;
376
413
                    aiEdges[poGDS->iYDim] = nYSize;
377
 
                                
 
414
 
378
415
                    aiStart[poGDS->iXDim] = nBlockXOff;
379
416
                    aiEdges[poGDS->iXDim] = nBlockXSize;
380
417
                    break;
381
418
                  case 3: // 3Dim: volume
382
419
                    aiStart[poGDS->iBandDim] = nBand - 1;
383
420
                    aiEdges[poGDS->iBandDim] = 1;
384
 
                                
 
421
 
385
422
                    aiStart[poGDS->iYDim] = nYOff;
386
423
                    aiEdges[poGDS->iYDim] = nYSize;
387
 
                                
 
424
 
388
425
                    aiStart[poGDS->iXDim] = nBlockXOff;
389
426
                    aiEdges[poGDS->iXDim] = nBlockXSize;
390
427
                    break;
391
428
                  case 2: // 2Dim: rows/cols
392
429
                    aiStart[poGDS->iYDim] = nYOff;
393
430
                    aiEdges[poGDS->iYDim] = nYSize;
394
 
                                
 
431
 
395
432
                    aiStart[poGDS->iXDim] = nBlockXOff;
396
433
                    aiEdges[poGDS->iXDim] = nBlockXSize;
397
434
                    break;
398
435
                }
399
 
                if( GDreadfield( hGD, poGDS->pszFieldName,
400
 
                                 aiStart, NULL, aiEdges, pImage ) < 0 )
401
 
                {
402
 
                    CPLError( CE_Failure, CPLE_AppDefined, 
 
436
 
 
437
                /* Ensure that we don't overlap the bottom or right edges */
 
438
                /* of the dataset in order to use the GDreadtile() API */
 
439
                if ( poGDS->bReadTile &&
 
440
                     (nBlockXOff + 1) * nBlockXSize <= nRasterXSize &&
 
441
                     (nBlockYOff + 1) * nBlockYSize <= nRasterYSize )
 
442
                {
 
443
                    int32 tilecoords[] = { nBlockYOff , nBlockXOff };
 
444
                    if( GDreadtile( hGD, poGDS->pszFieldName , tilecoords , pImage ) != 0 )
 
445
                    {
 
446
                        CPLError( CE_Failure, CPLE_AppDefined, "GDreadtile() failed for block." );
 
447
                        eErr = CE_Failure;
 
448
                    }
 
449
                }
 
450
                else if( GDreadfield( hGD, poGDS->pszFieldName,
 
451
                                aiStart, NULL, aiEdges, pImage ) < 0 )
 
452
                {
 
453
                    CPLError( CE_Failure, CPLE_AppDefined,
403
454
                              "GDreadfield() failed for block." );
404
455
                    eErr = CE_Failure;
405
456
                }
406
457
                GDdetach( hGD );
407
458
            }
408
459
            break;
409
 
                    
 
460
 
410
461
            case H4ST_EOS_SWATH:
411
462
            case H4ST_EOS_SWATH_GEOL:
412
463
            {
419
470
                  case 3: // 3Dim: volume
420
471
                    aiStart[poGDS->iBandDim] = nBand - 1;
421
472
                    aiEdges[poGDS->iBandDim] = 1;
422
 
                                
 
473
 
423
474
                    aiStart[poGDS->iYDim] = nYOff;
424
475
                    aiEdges[poGDS->iYDim] = nYSize;
425
 
                                
 
476
 
426
477
                    aiStart[poGDS->iXDim] = nBlockXOff;
427
478
                    aiEdges[poGDS->iXDim] = nBlockXSize;
428
479
                    break;
429
480
                  case 2: // 2Dim: rows/cols
430
481
                    aiStart[poGDS->iYDim] = nYOff;
431
482
                    aiEdges[poGDS->iYDim] = nYSize;
432
 
                                
 
483
 
433
484
                    aiStart[poGDS->iXDim] = nBlockXOff;
434
485
                    aiEdges[poGDS->iXDim] = nBlockXSize;
435
486
                    break;
437
488
                if( SWreadfield( hSW, poGDS->pszFieldName,
438
489
                                 aiStart, NULL, aiEdges, pImage ) < 0 )
439
490
                {
440
 
                    CPLError( CE_Failure, CPLE_AppDefined, 
 
491
                    CPLError( CE_Failure, CPLE_AppDefined,
441
492
                              "SWreadfield() failed for block." );
442
493
                    eErr = CE_Failure;
443
494
                }
480
531
    int nYOff = nBlockYOff * nBlockYSize;
481
532
    int nYSize = MIN(nYOff + nBlockYSize, poDS->GetRasterYSize()) - nYOff;
482
533
    CPLAssert( nBlockXOff == 0 );
483
 
    
 
534
 
484
535
/* -------------------------------------------------------------------- */
485
536
/*      Process based on rank.                                          */
486
537
/* -------------------------------------------------------------------- */
492
543
 
493
544
                aiStart[poGDS->iBandDim] = nBand - 1;
494
545
                aiEdges[poGDS->iBandDim] = 1;
495
 
            
 
546
 
496
547
                aiStart[poGDS->iYDim] = nYOff;
497
548
                aiEdges[poGDS->iYDim] = nYSize;
498
 
            
 
549
 
499
550
                aiStart[poGDS->iXDim] = nBlockXOff;
500
551
                aiEdges[poGDS->iXDim] = nBlockXSize;
501
552
 
512
563
                int32 iSDS = SDselect( poGDS->hSD, nBand - 1 );
513
564
                aiStart[poGDS->iYDim] = nYOff;
514
565
                aiEdges[poGDS->iYDim] = nYSize;
515
 
            
 
566
 
516
567
                aiStart[poGDS->iXDim] = nBlockXOff;
517
568
                aiEdges[poGDS->iXDim] = nBlockXSize;
518
569
 
623
674
double HDF4ImageRasterBand::GetOffset( int *pbSuccess )
624
675
 
625
676
{
626
 
    if( bHaveScaleAndOffset )
 
677
    if( bHaveOffset )
627
678
    {
628
679
        if( pbSuccess != NULL )
629
680
            *pbSuccess = TRUE;
640
691
double HDF4ImageRasterBand::GetScale( int *pbSuccess )
641
692
 
642
693
{
643
 
    if( bHaveScaleAndOffset )
 
694
    if( bHaveScale )
644
695
    {
645
696
        if( pbSuccess != NULL )
646
697
            *pbSuccess = TRUE;
700
751
    nGCPCount = 0;
701
752
 
702
753
    iDatasetType = HDF4_UNKNOWN;
 
754
    iSDS = FAIL;
 
755
 
 
756
    nBlockPreferredXSize = -1;
 
757
    nBlockPreferredYSize = -1;
 
758
    bReadTile = false;
 
759
 
703
760
}
704
761
 
705
762
/************************************************************************/
709
766
HDF4ImageDataset::~HDF4ImageDataset()
710
767
{
711
768
    FlushCache();
712
 
    
 
769
 
713
770
    if ( pszFilename )
714
771
        CPLFree( pszFilename );
 
772
    if ( iSDS != FAIL )
 
773
        SDendaccess( iSDS );
715
774
    if ( hSD > 0 )
716
775
        SDend( hSD );
717
776
    hSD = 0;
759
818
                        GDclose( hHDF4 );
760
819
                    default:
761
820
                        break;
762
 
                        
763
821
                }
764
822
                break;
765
823
            case HDF4_SDS:
864
922
    int         iBand;
865
923
    char        *pszName;
866
924
    const char  *pszValue;
867
 
    
 
925
 
868
926
    GDALDataset::FlushCache();
869
927
 
870
928
    if( eAccess == GA_ReadOnly )
902
960
        {
903
961
            pszValue = CPLParseNameValue( *papszMeta++, &pszName );
904
962
            if ( (SDsetattr( hSD, pszName, DFNT_CHAR8,
905
 
                             strlen(pszValue) + 1, pszValue )) < 0 );
 
963
                             strlen(pszValue) + 1, pszValue )) < 0 )
906
964
            {
907
965
                CPLDebug( "HDF4Image",
908
966
                          "Cannot write metadata information to output file");
921
979
        if ( poBand->bNoDataSet )
922
980
        {
923
981
            pszName = CPLStrdup( CPLSPrintf( "NoDataValue%d", iBand ) );
924
 
            pszValue = CPLSPrintf( "%lf", poBand->dfNoDataValue );
 
982
            pszValue = CPLSPrintf( "%f", poBand->dfNoDataValue );
925
983
            if ( (SDsetattr( hSD, pszName, DFNT_CHAR8,
926
984
                             strlen(pszValue) + 1, pszValue )) < 0 )
927
985
                {
993
1051
    OGRSpatialReference *poLatLong = NULL;
994
1052
    poLatLong = oSRS.CloneGeogCS();
995
1053
    poTransform = OGRCreateCoordinateTransformation( poLatLong, &oSRS );
996
 
    
 
1054
 
997
1055
    if( poTransform != NULL )
998
1056
        poTransform->Transform( 1, pdfGeoX, pdfGeoY, NULL );
999
 
        
 
1057
 
1000
1058
    if( poTransform != NULL )
1001
1059
        delete poTransform;
1002
1060
 
1019
1077
}
1020
1078
 
1021
1079
/************************************************************************/
 
1080
/*                         CaptureL1GMTLInfo()                          */
 
1081
/************************************************************************/
 
1082
 
 
1083
/*  FILE L71002025_02520010722_M_MTL.L1G
 
1084
 
 
1085
GROUP = L1_METADATA_FILE
 
1086
  ...
 
1087
  GROUP = PRODUCT_METADATA
 
1088
    PRODUCT_TYPE = "L1G"
 
1089
    PROCESSING_SOFTWARE = "IAS_5.1"
 
1090
    EPHEMERIS_TYPE = "DEFINITIVE"
 
1091
    SPACECRAFT_ID = "Landsat7"
 
1092
    SENSOR_ID = "ETM+"
 
1093
    ACQUISITION_DATE = 2001-07-22
 
1094
    WRS_PATH = 002
 
1095
    STARTING_ROW = 025
 
1096
    ENDING_ROW = 025
 
1097
    BAND_COMBINATION = "12345--7-"
 
1098
    PRODUCT_UL_CORNER_LAT = 51.2704805
 
1099
    PRODUCT_UL_CORNER_LON = -53.8914311
 
1100
    PRODUCT_UR_CORNER_LAT = 50.8458100
 
1101
    PRODUCT_UR_CORNER_LON = -50.9869091
 
1102
    PRODUCT_LL_CORNER_LAT = 49.6960897
 
1103
    PRODUCT_LL_CORNER_LON = -54.4047933
 
1104
    PRODUCT_LR_CORNER_LAT = 49.2841436
 
1105
    PRODUCT_LR_CORNER_LON = -51.5900428
 
1106
    PRODUCT_UL_CORNER_MAPX = 298309.894
 
1107
    PRODUCT_UL_CORNER_MAPY = 5683875.631
 
1108
    PRODUCT_UR_CORNER_MAPX = 500921.624
 
1109
    PRODUCT_UR_CORNER_MAPY = 5632678.683
 
1110
    PRODUCT_LL_CORNER_MAPX = 254477.193
 
1111
    PRODUCT_LL_CORNER_MAPY = 5510407.880
 
1112
    PRODUCT_LR_CORNER_MAPX = 457088.923
 
1113
    PRODUCT_LR_CORNER_MAPY = 5459210.932
 
1114
    PRODUCT_SAMPLES_REF = 6967
 
1115
    PRODUCT_LINES_REF = 5965
 
1116
    BAND1_FILE_NAME = "L71002025_02520010722_B10.L1G"
 
1117
    BAND2_FILE_NAME = "L71002025_02520010722_B20.L1G"
 
1118
    BAND3_FILE_NAME = "L71002025_02520010722_B30.L1G"
 
1119
    BAND4_FILE_NAME = "L71002025_02520010722_B40.L1G"
 
1120
    BAND5_FILE_NAME = "L71002025_02520010722_B50.L1G"
 
1121
    BAND7_FILE_NAME = "L72002025_02520010722_B70.L1G"
 
1122
    METADATA_L1_FILE_NAME = "L71002025_02520010722_MTL.L1G"
 
1123
    CPF_FILE_NAME = "L7CPF20010701_20010930_06"
 
1124
    HDF_DIR_FILE_NAME = "L71002025_02520010722_HDF.L1G"
 
1125
  END_GROUP = PRODUCT_METADATA
 
1126
  ...
 
1127
  GROUP = PROJECTION_PARAMETERS
 
1128
    REFERENCE_DATUM = "NAD83"
 
1129
    REFERENCE_ELLIPSOID = "GRS80"
 
1130
    GRID_CELL_SIZE_PAN = 15.000000
 
1131
    GRID_CELL_SIZE_THM = 60.000000
 
1132
    GRID_CELL_SIZE_REF = 30.000000
 
1133
    ORIENTATION = "NOM"
 
1134
    RESAMPLING_OPTION = "CC"
 
1135
    MAP_PROJECTION = "UTM"
 
1136
  END_GROUP = PROJECTION_PARAMETERS
 
1137
  GROUP = UTM_PARAMETERS
 
1138
    ZONE_NUMBER = 22
 
1139
  END_GROUP = UTM_PARAMETERS
 
1140
END_GROUP = L1_METADATA_FILE
 
1141
END
 
1142
*/
 
1143
 
 
1144
void HDF4ImageDataset::CaptureL1GMTLInfo()
 
1145
 
 
1146
{
 
1147
/* -------------------------------------------------------------------- */
 
1148
/*      Does the physical file look like it matches our expected        */
 
1149
/*      name pattern?                                                   */
 
1150
/* -------------------------------------------------------------------- */
 
1151
    if( strlen(pszFilename) < 8
 
1152
        || !EQUAL(pszFilename+strlen(pszFilename)-8,"_HDF.L1G") )
 
1153
        return;
 
1154
 
 
1155
/* -------------------------------------------------------------------- */
 
1156
/*      Construct the name of the corresponding MTL file.  We should    */
 
1157
/*      likely be able to extract that from the HDF itself but I'm      */
 
1158
/*      not sure where to find it.                                      */
 
1159
/* -------------------------------------------------------------------- */
 
1160
    CPLString osMTLFilename = pszFilename;
 
1161
    osMTLFilename.resize(osMTLFilename.length() - 8);
 
1162
    osMTLFilename += "_MTL.L1G";
 
1163
 
 
1164
/* -------------------------------------------------------------------- */
 
1165
/*      Ingest the MTL using the NASAKeywordHandler written for the     */
 
1166
/*      PDS driver.                                                     */
 
1167
/* -------------------------------------------------------------------- */
 
1168
    NASAKeywordHandler oMTL;
 
1169
 
 
1170
    VSILFILE *fp = VSIFOpenL( osMTLFilename, "r" );
 
1171
    if( fp == NULL )
 
1172
        return;
 
1173
 
 
1174
    if( !oMTL.Ingest( fp, 0 ) )
 
1175
    {
 
1176
        VSIFCloseL( fp );
 
1177
        return;
 
1178
    }
 
1179
 
 
1180
    VSIFCloseL( fp );
 
1181
 
 
1182
/* -------------------------------------------------------------------- */
 
1183
/*  Note: Different variation of MTL files use different group names.   */
 
1184
/*            Check for LPGS_METADATA_FILE and L1_METADATA_FILE.        */
 
1185
/* -------------------------------------------------------------------- */
 
1186
    double dfULX, dfULY, dfLRX, dfLRY;
 
1187
    double dfLLX, dfLLY, dfURX, dfURY;
 
1188
    CPLString osPrefix;
 
1189
 
 
1190
    if( oMTL.GetKeyword( "LPGS_METADATA_FILE.PRODUCT_METADATA"
 
1191
                         ".PRODUCT_UL_CORNER_LON", NULL ) )
 
1192
        osPrefix = "LPGS_METADATA_FILE.PRODUCT_METADATA.PRODUCT_";
 
1193
    else if( oMTL.GetKeyword( "L1_METADATA_FILE.PRODUCT_METADATA"
 
1194
                              ".PRODUCT_UL_CORNER_LON", NULL ) )
 
1195
        osPrefix = "L1_METADATA_FILE.PRODUCT_METADATA.PRODUCT_";
 
1196
    else
 
1197
        return;
 
1198
 
 
1199
    dfULX = atof(oMTL.GetKeyword((osPrefix+"UL_CORNER_LON").c_str(), "0" ));
 
1200
    dfULY = atof(oMTL.GetKeyword((osPrefix+"UL_CORNER_LAT").c_str(), "0" ));
 
1201
    dfLRX = atof(oMTL.GetKeyword((osPrefix+"LR_CORNER_LON").c_str(), "0" ));
 
1202
    dfLRY = atof(oMTL.GetKeyword((osPrefix+"LR_CORNER_LAT").c_str(), "0" ));
 
1203
    dfLLX = atof(oMTL.GetKeyword((osPrefix+"LL_CORNER_LON").c_str(), "0" ));
 
1204
    dfLLY = atof(oMTL.GetKeyword((osPrefix+"LL_CORNER_LAT").c_str(), "0" ));
 
1205
    dfURX = atof(oMTL.GetKeyword((osPrefix+"UR_CORNER_LON").c_str(), "0" ));
 
1206
    dfURY = atof(oMTL.GetKeyword((osPrefix+"UR_CORNER_LAT").c_str(), "0" ));
 
1207
 
 
1208
    CPLFree( pszGCPProjection );
 
1209
    pszGCPProjection = CPLStrdup( "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9108\"]],AXIS[\"Lat\",NORTH],AXIS[\"Long\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]" );
 
1210
 
 
1211
    nGCPCount = 4;
 
1212
    pasGCPList = (GDAL_GCP *) CPLCalloc( nGCPCount, sizeof( GDAL_GCP ) );
 
1213
    GDALInitGCPs( nGCPCount, pasGCPList );
 
1214
 
 
1215
    pasGCPList[0].dfGCPX = dfULX;
 
1216
    pasGCPList[0].dfGCPY = dfULY;
 
1217
    pasGCPList[0].dfGCPPixel = 0.0;
 
1218
    pasGCPList[0].dfGCPLine = 0.0;
 
1219
 
 
1220
    pasGCPList[1].dfGCPX = dfURX;
 
1221
    pasGCPList[1].dfGCPY = dfURY;
 
1222
    pasGCPList[1].dfGCPPixel = GetRasterXSize();
 
1223
    pasGCPList[1].dfGCPLine = 0.0;
 
1224
 
 
1225
    pasGCPList[2].dfGCPX = dfLLX;
 
1226
    pasGCPList[2].dfGCPY = dfLLY;
 
1227
    pasGCPList[2].dfGCPPixel = 0.0;
 
1228
    pasGCPList[2].dfGCPLine = GetRasterYSize();
 
1229
 
 
1230
    pasGCPList[3].dfGCPX = dfLRX;
 
1231
    pasGCPList[3].dfGCPY = dfLRY;
 
1232
    pasGCPList[3].dfGCPPixel = GetRasterXSize();
 
1233
    pasGCPList[3].dfGCPLine = GetRasterYSize();
 
1234
}
 
1235
 
 
1236
/************************************************************************/
1022
1237
/*                       CaptureNRLGeoTransform()                       */
1023
1238
/*                                                                      */
1024
1239
/*      Capture geotransform and coordinate system from NRL (Naval      */
1090
1305
 
1091
1306
    for( iCorner = 0; iCorner < 4; iCorner++ )
1092
1307
    {
1093
 
        const char *pszCornerLoc = 
 
1308
        const char *pszCornerLoc =
1094
1309
            CSLFetchNameValue( papszGlobalMetadata, apszItems[iCorner] );
1095
1310
 
1096
1311
        if( pszCornerLoc == NULL )
1104
1319
        adfXY[iCorner*2+0] = CPLAtof( papszTokens[1] );
1105
1320
        adfXY[iCorner*2+1] = CPLAtof( papszTokens[0] );
1106
1321
 
1107
 
        if( adfXY[iCorner*2+0] < -360 || adfXY[iCorner*2+0] > 360 
 
1322
        if( adfXY[iCorner*2+0] < -360 || adfXY[iCorner*2+0] > 360
1108
1323
            || adfXY[iCorner*2+1] < -90 || adfXY[iCorner*2+1] > 90 )
1109
1324
            bLLPossible = FALSE;
1110
1325
 
1114
1329
/* -------------------------------------------------------------------- */
1115
1330
/*      Does this look like nice clean "northup" lat/long data?         */
1116
1331
/* -------------------------------------------------------------------- */
1117
 
    if( adfXY[0*2+0] == adfXY[2*2+0] && adfXY[0*2+1] == adfXY[1*2+1] 
 
1332
    if( adfXY[0*2+0] == adfXY[2*2+0] && adfXY[0*2+1] == adfXY[1*2+1]
1118
1333
        && bLLPossible )
1119
1334
    {
1120
1335
        bHasGeoTransform = TRUE;
1124
1339
        adfGeoTransform[3] = adfXY[0*2+1];
1125
1340
        adfGeoTransform[4] = 0.0;
1126
1341
        adfGeoTransform[5] = (adfXY[2*2+1] - adfXY[0*2+1]) / nRasterYSize;
1127
 
        
 
1342
 
1128
1343
        oSRS.SetWellKnownGeogCS( "WGS84" );
1129
1344
        CPLFree( pszProjection );
1130
1345
        oSRS.exportToWkt( &pszProjection );
1135
1350
/* -------------------------------------------------------------------- */
1136
1351
    int  bGotGCTPProjection = FALSE;
1137
1352
    int  iSDSIndex = FAIL, iSDS = FAIL;
1138
 
    const char *mapProjection = CSLFetchNameValue( papszGlobalMetadata, 
 
1353
    const char *mapProjection = CSLFetchNameValue( papszGlobalMetadata,
1139
1354
                                                   "mapProjection" );
1140
 
    
 
1355
 
1141
1356
    if( mapProjection )
1142
1357
        iSDSIndex = SDnametoindex( hSD, mapProjection );
1143
1358
 
1144
1359
    if( iSDSIndex != FAIL )
1145
1360
        iSDS = SDselect( hSD, iSDSIndex );
1146
 
       
 
1361
 
1147
1362
    if( iSDS != FAIL )
1148
1363
    {
1149
1364
        char        szName[HDF4_SDS_MAXNAMELEN];
1156
1371
        aiStart[0] = 0;
1157
1372
        aiEdges[0] = 29;
1158
1373
 
1159
 
        if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType, 
1160
 
                       &nAttrs) == 0 
1161
 
            && iNumType == DFNT_FLOAT64 
 
1374
        if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType,
 
1375
                       &nAttrs) == 0
 
1376
            && iNumType == DFNT_FLOAT64
1162
1377
            && iRank == 1
1163
 
            && aiDimSizes[0] >= 29 
1164
 
            && SDreaddata( iSDS, aiStart, NULL, aiEdges, adfGCTP ) == 0 
1165
 
            && oSRS.importFromUSGS( (long) adfGCTP[1], (long) adfGCTP[2], 
1166
 
                                    adfGCTP+4, 
 
1378
            && aiDimSizes[0] >= 29
 
1379
            && SDreaddata( iSDS, aiStart, NULL, aiEdges, adfGCTP ) == 0
 
1380
            && oSRS.importFromUSGS( (long) adfGCTP[1], (long) adfGCTP[2],
 
1381
                                    adfGCTP+4,
1167
1382
                                    (long) adfGCTP[3] ) == OGRERR_NONE )
1168
1383
        {
1169
1384
            CPLDebug( "HDF4Image", "GCTP Parms = %g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g",
1170
 
                      adfGCTP[0], 
1171
 
                      adfGCTP[1], 
1172
 
                      adfGCTP[2], 
1173
 
                      adfGCTP[3], 
1174
 
                      adfGCTP[4], 
1175
 
                      adfGCTP[5], 
1176
 
                      adfGCTP[6], 
1177
 
                      adfGCTP[7], 
1178
 
                      adfGCTP[8], 
1179
 
                      adfGCTP[9], 
1180
 
                      adfGCTP[10], 
1181
 
                      adfGCTP[11], 
1182
 
                      adfGCTP[12], 
1183
 
                      adfGCTP[13], 
1184
 
                      adfGCTP[14], 
1185
 
                      adfGCTP[15], 
1186
 
                      adfGCTP[16], 
1187
 
                      adfGCTP[17], 
1188
 
                      adfGCTP[18], 
1189
 
                      adfGCTP[19], 
1190
 
                      adfGCTP[20], 
1191
 
                      adfGCTP[21], 
1192
 
                      adfGCTP[22], 
1193
 
                      adfGCTP[23], 
1194
 
                      adfGCTP[24], 
1195
 
                      adfGCTP[25], 
1196
 
                      adfGCTP[26], 
1197
 
                      adfGCTP[27], 
 
1385
                      adfGCTP[0],
 
1386
                      adfGCTP[1],
 
1387
                      adfGCTP[2],
 
1388
                      adfGCTP[3],
 
1389
                      adfGCTP[4],
 
1390
                      adfGCTP[5],
 
1391
                      adfGCTP[6],
 
1392
                      adfGCTP[7],
 
1393
                      adfGCTP[8],
 
1394
                      adfGCTP[9],
 
1395
                      adfGCTP[10],
 
1396
                      adfGCTP[11],
 
1397
                      adfGCTP[12],
 
1398
                      adfGCTP[13],
 
1399
                      adfGCTP[14],
 
1400
                      adfGCTP[15],
 
1401
                      adfGCTP[16],
 
1402
                      adfGCTP[17],
 
1403
                      adfGCTP[18],
 
1404
                      adfGCTP[19],
 
1405
                      adfGCTP[20],
 
1406
                      adfGCTP[21],
 
1407
                      adfGCTP[22],
 
1408
                      adfGCTP[23],
 
1409
                      adfGCTP[24],
 
1410
                      adfGCTP[25],
 
1411
                      adfGCTP[26],
 
1412
                      adfGCTP[27],
1198
1413
                      adfGCTP[28] );
1199
 
            
 
1414
 
1200
1415
            CPLFree( pszProjection );
1201
1416
            oSRS.exportToWkt( &pszProjection );
1202
1417
            bGotGCTPProjection = TRUE;
1203
1418
        }
 
1419
 
 
1420
        SDendaccess(iSDS);
1204
1421
    }
1205
1422
 
1206
1423
/* -------------------------------------------------------------------- */
1215
1432
 
1216
1433
        oWGS84.SetWellKnownGeogCS( "WGS84" );
1217
1434
 
1218
 
        OGRCoordinateTransformation *poCT = 
 
1435
        OGRCoordinateTransformation *poCT =
1219
1436
            OGRCreateCoordinateTransformation( &oWGS84, &oSRS );
1220
1437
 
1221
1438
        dfULX = adfXY[0*2+0];
1222
1439
        dfULY = adfXY[0*2+1];
1223
 
        
 
1440
 
1224
1441
        dfLRX = adfXY[3*2+0];
1225
1442
        dfLRY = adfXY[3*2+1];
1226
 
        
1227
 
        if( poCT->Transform( 1, &dfULX, &dfULY ) 
 
1443
 
 
1444
        if( poCT->Transform( 1, &dfULX, &dfULY )
1228
1445
            && poCT->Transform( 1, &dfLRX, &dfLRY ) )
1229
1446
        {
1230
1447
            bHasGeoTransform = TRUE;
1286
1503
void HDF4ImageDataset::CaptureCoastwatchGCTPInfo()
1287
1504
 
1288
1505
{
1289
 
    if( CSLFetchNameValue( papszGlobalMetadata, "gctp_sys" ) == NULL 
1290
 
        || CSLFetchNameValue( papszGlobalMetadata, "gctp_zone" ) == NULL 
1291
 
        || CSLFetchNameValue( papszGlobalMetadata, "gctp_parm" ) == NULL 
1292
 
        || CSLFetchNameValue( papszGlobalMetadata, "gctp_datum" ) == NULL 
 
1506
    if( CSLFetchNameValue( papszGlobalMetadata, "gctp_sys" ) == NULL
 
1507
        || CSLFetchNameValue( papszGlobalMetadata, "gctp_zone" ) == NULL
 
1508
        || CSLFetchNameValue( papszGlobalMetadata, "gctp_parm" ) == NULL
 
1509
        || CSLFetchNameValue( papszGlobalMetadata, "gctp_datum" ) == NULL
1293
1510
        || CSLFetchNameValue( papszGlobalMetadata, "et_affine" ) == NULL )
1294
1511
        return;
1295
1512
 
1304
1521
    nZone = atoi( CSLFetchNameValue( papszGlobalMetadata, "gctp_zone" ) );
1305
1522
    nDatum = atoi( CSLFetchNameValue( papszGlobalMetadata, "gctp_datum" ) );
1306
1523
 
1307
 
    papszTokens = CSLTokenizeStringComplex( 
1308
 
        CSLFetchNameValue( papszGlobalMetadata, "gctp_parm" ), ",", 
 
1524
    papszTokens = CSLTokenizeStringComplex(
 
1525
        CSLFetchNameValue( papszGlobalMetadata, "gctp_parm" ), ",",
1309
1526
        FALSE, FALSE );
1310
1527
    if( CSLCount(papszTokens) < 15 )
1311
1528
        return;
1327
1544
/* -------------------------------------------------------------------- */
1328
1545
/*      Capture the affine transform info.                              */
1329
1546
/* -------------------------------------------------------------------- */
1330
 
    
1331
 
    papszTokens = CSLTokenizeStringComplex( 
1332
 
        CSLFetchNameValue( papszGlobalMetadata, "et_affine" ), ",", 
 
1547
 
 
1548
    papszTokens = CSLTokenizeStringComplex(
 
1549
        CSLFetchNameValue( papszGlobalMetadata, "et_affine" ), ",",
1333
1550
        FALSE, FALSE );
1334
1551
    if( CSLCount(papszTokens) != 6 )
1335
1552
        return;
1336
1553
 
1337
 
    // We don't seem to have proper ef_affine docs so I don't 
1338
 
    // know which of these two coefficients goes where. 
 
1554
    // We don't seem to have proper ef_affine docs so I don't
 
1555
    // know which of these two coefficients goes where.
1339
1556
    if( CPLAtof(papszTokens[0]) != 0.0 || CPLAtof(papszTokens[3]) != 0.0 )
1340
1557
        return;
1341
 
        
 
1558
 
1342
1559
    bHasGeoTransform = TRUE;
1343
1560
    adfGeoTransform[0] = CPLAtof( papszTokens[4] );
1344
1561
    adfGeoTransform[1] = CPLAtof( papszTokens[2] );
1347
1564
    adfGeoTransform[4] = 0.0;
1348
1565
    adfGeoTransform[5] = CPLAtof( papszTokens[1] );
1349
1566
 
1350
 
    // Middle of pixel adjustment. 
 
1567
    // Middle of pixel adjustment.
1351
1568
    adfGeoTransform[0] -= adfGeoTransform[1] * 0.5;
1352
1569
    adfGeoTransform[3] -= adfGeoTransform[5] * 0.5;
1353
1570
}
1403
1620
    }
1404
1621
 
1405
1622
    // If didn't get a band dimension yet, but have an extra
1406
 
    // dimension, use it as the band dimension. 
 
1623
    // dimension, use it as the band dimension.
1407
1624
    if ( iRank > 2 && iBandDim == -1 )
1408
1625
    {
1409
1626
        if( iXDim != 0 && iYDim != 0 )
1449
1666
        SWinqattrs( hSW, pszAttrList, &nStrBufSize );
1450
1667
 
1451
1668
#ifdef DEBUG
1452
 
        CPLDebug( "HDF4Image", "List of attributes in swath %s: %s",
 
1669
        CPLDebug( "HDF4Image", "List of attributes in swath \"%s\": %s",
1453
1670
                  pszFieldName, pszAttrList );
1454
1671
#endif
1455
1672
 
1475
1692
            {
1476
1693
                ((char *)pData)[nValues] = '\0';
1477
1694
                papszLocalMetadata = CSLAddNameValue( papszLocalMetadata,
1478
 
                                                      papszAttributes[i], 
 
1695
                                                      papszAttributes[i],
1479
1696
                                                      (const char *) pData );
1480
1697
            }
1481
1698
            else
1488
1705
                if ( pszTemp )
1489
1706
                    CPLFree( pszTemp );
1490
1707
            }
1491
 
            
 
1708
 
1492
1709
            if ( pData )
1493
1710
                CPLFree( pData );
1494
1711
 
1508
1725
        int32       iRank, iNumType, iAttribute, nAttrs, nValues;
1509
1726
        char        szName[HDF4_SDS_MAXNAMELEN];
1510
1727
        int32       aiDimSizes[H4_MAX_VAR_DIMS];
1511
 
        
1512
 
        if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType, 
 
1728
 
 
1729
        if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType,
1513
1730
                       &nAttrs) == 0 )
1514
1731
        {
1515
1732
            for ( iAttribute = 0; iAttribute < nAttrs; iAttribute++ )
1598
1815
                if ( pszTemp )
1599
1816
                    CPLFree( pszTemp );
1600
1817
            }
1601
 
            
 
1818
 
1602
1819
            if ( pData )
1603
1820
                CPLFree( pData );
1604
1821
 
1618
1835
        int32       iRank, iNumType, iAttribute, nAttrs, nValues;
1619
1836
        char        szName[HDF4_SDS_MAXNAMELEN];
1620
1837
        int32       aiDimSizes[H4_MAX_VAR_DIMS];
1621
 
        
1622
 
        if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType, 
 
1838
 
 
1839
        if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType,
1623
1840
                       &nAttrs) == 0 )
1624
1841
        {
1625
1842
            for ( iAttribute = 0; iAttribute < nAttrs; iAttribute++ )
1681
1898
        int32       iRank, iNumType, nAttrs, iSDS;
1682
1899
        char        szName[HDF4_SDS_MAXNAMELEN];
1683
1900
        int32       aiDimSizes[H4_MAX_VAR_DIMS];
1684
 
        
 
1901
 
1685
1902
        iSDS = SDselect( hSD, iDSIndex );
1686
1903
 
1687
 
        if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType, 
1688
 
                       &nAttrs) != 0 )
1689
 
            continue;
1690
 
 
1691
 
        if( EQUAL(szName,"latitude") )
1692
 
            iYIndex = iDSIndex;
1693
 
 
1694
 
        if( EQUAL(szName,"longitude") )
1695
 
            iXIndex = iDSIndex;
1696
 
 
1697
 
        
 
1904
        if( SDgetinfo( iSDS, szName, &iRank, aiDimSizes, &iNumType,
 
1905
                       &nAttrs) == 0 )
 
1906
        {
 
1907
            if( EQUAL(szName,"latitude") )
 
1908
                iYIndex = iDSIndex;
 
1909
 
 
1910
            if( EQUAL(szName,"longitude") )
 
1911
                iXIndex = iDSIndex;
 
1912
        }
 
1913
 
 
1914
        SDendaccess(iSDS);
1698
1915
    }
1699
1916
 
1700
1917
    if( iXIndex == -1 || iYIndex == -1 )
1704
1921
/*      We found geolocation information.  Record it as metadata.       */
1705
1922
/* -------------------------------------------------------------------- */
1706
1923
    CPLString  osWrk;
1707
 
    
 
1924
 
1708
1925
    SetMetadataItem( "SRS", SRS_WKT_WGS84, "GEOLOCATION" );
1709
 
    
1710
 
    osWrk.Printf( "HDF4_SDS:UNKNOWN:\"%s\":%d", 
 
1926
 
 
1927
    osWrk.Printf( "HDF4_SDS:UNKNOWN:\"%s\":%d",
1711
1928
                  pszFilename, iXIndex );
1712
1929
    SetMetadataItem( "X_DATASET", osWrk, "GEOLOCATION" );
1713
1930
    SetMetadataItem( "X_BAND", "1" , "GEOLOCATION" );
1714
1931
 
1715
 
    osWrk.Printf( "HDF4_SDS:UNKNOWN:\"%s\":%d", 
 
1932
    osWrk.Printf( "HDF4_SDS:UNKNOWN:\"%s\":%d",
1716
1933
                  pszFilename, iYIndex );
1717
1934
    SetMetadataItem( "Y_DATASET", osWrk, "GEOLOCATION" );
1718
1935
    SetMetadataItem( "Y_BAND", "1" , "GEOLOCATION" );
1739
1956
    char    szYGeo[8192] = "";
1740
1957
    char    szPixel[8192]= "";
1741
1958
    char    szLine[8192] = "";
1742
 
    char    *pszGeoList = NULL;
1743
 
    char    szGeoDimList[8192] = "";
1744
1959
    int32   iWrkNumType;
1745
 
    int32   nDataFields, nDimMaps;
1746
1960
    void    *pLat = NULL, *pLong = NULL;
1747
1961
    void    *pLatticeX = NULL, *pLatticeY = NULL;
1748
1962
    int32   iLatticeType, iLatticeDataSize = 0, iRank;
1749
1963
    int32   nLatCount = 0, nLongCount = 0;
1750
1964
    int32   nXPoints=0, nYPoints=0;
1751
1965
    int32   nStrBufSize;
1752
 
    int32   aiDimSizes[H4_MAX_VAR_DIMS];
1753
1966
    int     i, j, iDataSize = 0, iPixelDim=-1,iLineDim=-1, iLongDim=-1, iLatDim=-1;
1754
1967
    int32   *paiRank = NULL, *paiNumType = NULL,
1755
1968
        *paiOffset = NULL, *paiIncrement = NULL;
1756
 
    char    **papszGeolocations = NULL;
1757
 
    char    *pszDimMaps = NULL;
1758
1969
 
1759
1970
/* -------------------------------------------------------------------- */
1760
1971
/*  Determine a product name.                                           */
1783
1994
        else if ( EQUALN(pszProduct, "MOD02", 5)
1784
1995
                  || EQUALN(pszProduct, "MYD02", 5) )
1785
1996
            eProduct = PROD_MODIS_L1B;
 
1997
        else if ( EQUALN(pszProduct, "MOD07_L2", 8) )
 
1998
            eProduct = PROD_MODIS_L2;
1786
1999
    }
1787
 
    
1788
 
/* -------------------------------------------------------------------- */
1789
 
/*      xx                                                              */
1790
 
/* -------------------------------------------------------------------- */
1791
 
    nDataFields = SWnentries( hSW, HDFE_NENTGFLD, &nStrBufSize );
1792
 
    pszGeoList = (char *)CPLMalloc( nStrBufSize + 1 );
 
2000
 
 
2001
/* -------------------------------------------------------------------- */
 
2002
/*      Read names of geolocation fields and corresponding              */
 
2003
/*      geolocation maps.                                               */
 
2004
/* -------------------------------------------------------------------- */
 
2005
    int32   nDataFields = SWnentries( hSW, HDFE_NENTGFLD, &nStrBufSize );
 
2006
    char    *pszGeoList = (char *)CPLMalloc( nStrBufSize + 1 );
1793
2007
    paiRank = (int32 *)CPLMalloc( nDataFields * sizeof(int32) );
1794
2008
    paiNumType = (int32 *)CPLMalloc( nDataFields * sizeof(int32) );
 
2009
 
1795
2010
    if ( nDataFields !=
1796
2011
         SWinqgeofields(hSW, pszGeoList, paiRank, paiNumType) )
1797
2012
    {
1798
2013
        CPLDebug( "HDF4Image",
1799
 
                  "Can't get the list of geolocation fields in swath %s",
 
2014
                  "Can't get the list of geolocation fields in swath \"%s\"",
1800
2015
                  pszSubdatasetName );
1801
2016
    }
1802
2017
 
1805
2020
    {
1806
2021
        char    *pszTmp;
1807
2022
        CPLDebug( "HDF4Image",
1808
 
                  "Number of geolocation fields in swath %s: %ld",
 
2023
                  "Number of geolocation fields in swath \"%s\": %ld",
1809
2024
                  pszSubdatasetName, (long)nDataFields );
1810
2025
        CPLDebug( "HDF4Image",
1811
 
                  "List of geolocation fields in swath %s: %s",
 
2026
                  "List of geolocation fields in swath \"%s\": %s",
1812
2027
                  pszSubdatasetName, pszGeoList );
1813
2028
        pszTmp = SPrintArray( GDT_UInt32, paiRank,
1814
2029
                              nDataFields, "," );
1818
2033
    }
1819
2034
#endif
1820
2035
 
1821
 
    papszGeolocations = CSLTokenizeString2( pszGeoList, ",",
1822
 
                                            CSLT_HONOURSTRINGS );
1823
2036
    // Read geolocation data
1824
 
    nDimMaps = SWnentries( hSW, HDFE_NENTMAP, &nStrBufSize );
 
2037
    int32   nDimMaps = SWnentries( hSW, HDFE_NENTMAP, &nStrBufSize );
1825
2038
    if ( nDimMaps <= 0 )
1826
2039
    {
1827
 
        CPLDebug( "HDF4Image", "No geolocation maps in swath %s",
 
2040
 
 
2041
#ifdef DEBUG
 
2042
        CPLDebug( "HDF4Image", "No geolocation maps in swath \"%s\"",
1828
2043
                  pszSubdatasetName );
 
2044
        CPLDebug( "HDF4Image",
 
2045
                  "Suppose one-to-one mapping. X field is \"%s\", Y field is \"%s\"",
 
2046
                  papszDimList[iXDim], papszDimList[iYDim] );
 
2047
#endif
 
2048
 
 
2049
        strncpy( szPixel, papszDimList[iXDim], 8192 );
 
2050
        strncpy( szLine, papszDimList[iYDim], 8192 );
 
2051
        strncpy( szXGeo, papszDimList[iXDim], 8192 );
 
2052
        strncpy( szYGeo, papszDimList[iYDim], 8192 );
 
2053
        paiOffset = (int32 *)CPLCalloc( 2, sizeof(int32) );
 
2054
        paiIncrement = (int32 *)CPLCalloc( 2, sizeof(int32) );
 
2055
        paiOffset[0] = paiOffset[1] = 0;
 
2056
        paiIncrement[0] = paiIncrement[1] = 1;
1829
2057
    }
1830
2058
    else
1831
2059
    {
1832
 
        pszDimMaps = (char *)CPLMalloc( nStrBufSize + 1 );
1833
 
        paiOffset = (int32 *)CPLMalloc( nDimMaps * sizeof(int32) );
1834
 
        memset( paiOffset, 0, nDimMaps * sizeof(int32) );
1835
 
        paiIncrement = (int32 *)CPLMalloc( nDimMaps * sizeof(int32) );
1836
 
        memset( paiIncrement, 0, nDimMaps * sizeof(int32) );
1837
 
 
 
2060
        char *pszDimMaps = (char *)CPLMalloc( nStrBufSize + 1 );
 
2061
        paiOffset = (int32 *)CPLCalloc( nDimMaps, sizeof(int32) );
 
2062
        paiIncrement = (int32 *)CPLCalloc( nDimMaps, sizeof(int32) );
1838
2063
 
1839
2064
        *pszDimMaps = '\0';
1840
2065
        if ( nDimMaps != SWinqmaps(hSW, pszDimMaps, paiOffset, paiIncrement) )
1841
2066
        {
1842
2067
            CPLDebug( "HDF4Image",
1843
 
                      "Can't get the list of geolocation maps in swath %s",
 
2068
                      "Can't get the list of geolocation maps in swath \"%s\"",
1844
2069
                      pszSubdatasetName );
1845
2070
        }
1846
2071
 
1848
2073
        else
1849
2074
        {
1850
2075
            char    *pszTmp;
1851
 
                
 
2076
 
1852
2077
            CPLDebug( "HDF4Image",
1853
 
                      "List of geolocation maps in swath %s: %s",
 
2078
                      "List of geolocation maps in swath \"%s\": %s",
1854
2079
                      pszSubdatasetName, pszDimMaps );
1855
2080
            pszTmp = SPrintArray( GDT_Int32, paiOffset,
1856
2081
                                  nDimMaps, "," );
1865
2090
        }
1866
2091
#endif
1867
2092
 
1868
 
        char **papszDimMap;
1869
 
 
1870
 
        papszDimMap = CSLTokenizeString2( pszDimMaps, ",",
1871
 
                                          CSLT_HONOURSTRINGS );
1872
 
        CPLFree( pszDimMaps );
1873
 
        pszDimMaps = NULL;
1874
 
 
1875
 
        for ( i = 0; i < CSLCount(papszDimMap); i++ )
 
2093
        char    **papszDimMap = CSLTokenizeString2( pszDimMaps, ",",
 
2094
                                                    CSLT_HONOURSTRINGS );
 
2095
        int     nDimMapCount = CSLCount(papszDimMap);
 
2096
 
 
2097
        for ( i = 0; i < nDimMapCount; i++ )
1876
2098
        {
1877
2099
            if ( strstr(papszDimMap[i], papszDimList[iXDim]) )
1878
2100
            {
1891
2113
                    *pszTemp = '\0';
1892
2114
            }
1893
2115
        }
 
2116
 
1894
2117
        CSLDestroy( papszDimMap );
1895
 
        papszDimMap = NULL;
 
2118
        CPLFree( pszDimMaps );
1896
2119
    }
1897
2120
 
1898
 
    for ( i = 0; i < CSLCount(papszGeolocations); i++ )
 
2121
    if ( *szXGeo == 0 || *szYGeo == 0 )
 
2122
        return FALSE;
 
2123
 
 
2124
/* -------------------------------------------------------------------- */
 
2125
/*      Read geolocation fields.                                        */
 
2126
/* -------------------------------------------------------------------- */
 
2127
    char    szGeoDimList[8192] = "";
 
2128
    char    **papszGeolocations = CSLTokenizeString2( pszGeoList, ",",
 
2129
                                                      CSLT_HONOURSTRINGS );
 
2130
    int     nGeolocationsCount = CSLCount( papszGeolocations );
 
2131
    int32   aiDimSizes[H4_MAX_VAR_DIMS];
 
2132
 
 
2133
    for ( i = 0; i < nGeolocationsCount; i++ )
1899
2134
    {
1900
2135
        char    **papszGeoDimList = NULL;
1901
2136
 
1902
 
        SWfieldinfo( hSW, papszGeolocations[i], &iRank,
1903
 
                     aiDimSizes, &iWrkNumType, szGeoDimList );
 
2137
        if ( SWfieldinfo( hSW, papszGeolocations[i], &iRank,
 
2138
                          aiDimSizes, &iWrkNumType, szGeoDimList ) < 0 )
 
2139
        {
 
2140
 
 
2141
#ifdef DEBUG
 
2142
        CPLDebug( "HDF4Image",
 
2143
                  "Can't read attributes of geolocation field \"%s\"",
 
2144
                  papszGeolocations[i] );
 
2145
#endif
 
2146
 
 
2147
            return FALSE;
 
2148
        }
 
2149
 
1904
2150
        papszGeoDimList = CSLTokenizeString2( szGeoDimList,
1905
2151
                                              ",", CSLT_HONOURSTRINGS );
1906
2152
 
 
2153
        if ( CSLCount(papszGeoDimList) > H4_MAX_VAR_DIMS )
 
2154
        {
 
2155
            CSLDestroy( papszGeoDimList );
 
2156
            return FALSE;
 
2157
        }
 
2158
 
1907
2159
#ifdef DEBUG
1908
2160
        CPLDebug( "HDF4Image",
1909
 
                  "List of dimensions in geolocation field %s: %s",
 
2161
                  "List of dimensions in geolocation field \"%s\": %s",
1910
2162
                  papszGeolocations[i], szGeoDimList );
1911
2163
#endif
1912
2164
 
1913
 
        if (szXGeo[0] == 0 || szYGeo[0] == 0)
1914
 
            return FALSE;
1915
 
 
1916
2165
        nXPoints = aiDimSizes[CSLFindString( papszGeoDimList, szXGeo )];
1917
2166
        nYPoints = aiDimSizes[CSLFindString( papszGeoDimList, szYGeo )];
1918
2167
 
1960
2209
        }
1961
2210
 
1962
2211
        CSLDestroy( papszGeoDimList );
1963
 
 
1964
2212
    }
1965
2213
 
1966
2214
/* -------------------------------------------------------------------- */
1986
2234
 
1987
2235
        iStart[2] = 0;
1988
2236
        iEdges[2] = 1;
1989
 
                        
 
2237
 
1990
2238
        pLatticeX = CPLMalloc( nLatCount * iLatticeDataSize );
1991
2239
        if (SWreadfield( hSW, "LatticePoint", iStart, NULL,
1992
2240
                         iEdges, (VOIDP)pLatticeX ) < 0)
2013
2261
/* -------------------------------------------------------------------- */
2014
2262
/*      Determine whether to use no, partial or full GCPs.              */
2015
2263
/* -------------------------------------------------------------------- */
2016
 
    const char *pszGEOL_AS_GCPS = CPLGetConfigOption( "GEOL_AS_GCPS", 
 
2264
    const char *pszGEOL_AS_GCPS = CPLGetConfigOption( "GEOL_AS_GCPS",
2017
2265
                                                      "PARTIAL" );
2018
2266
    int iGCPStepX, iGCPStepY;
2019
2267
 
2060
2308
 
2061
2309
            char *pszProjLine =
2062
2310
                CPLStrdup(CPLSPrintf("MPMETHOD%s", pszBand));
2063
 
            char *pszParmsLine = 
 
2311
            char *pszParmsLine =
2064
2312
                CPLStrdup(CPLSPrintf("PROJECTIONPARAMETERS%s",
2065
2313
                                     pszBand));
2066
 
            char *pszZoneLine = 
 
2314
            char *pszZoneLine =
2067
2315
                CPLStrdup(CPLSPrintf("UTMZONECODE%s",
2068
2316
                                     pszBand));
2069
2317
            char *pszEllipsoidLine =
2105
2353
            char **papszEllipsoid = (pszEllipsoid) ?
2106
2354
                CSLTokenizeString2( pszEllipsoid, ",",
2107
2355
                                    CSLT_HONOURSTRINGS ) : NULL;
2108
 
                            
 
2356
 
2109
2357
            long iEllipsoid = 8L; // WGS84 by default
2110
2358
            if ( papszEllipsoid
2111
2359
                 && CSLCount(papszEllipsoid) > 0 )
2113
2361
                if (EQUAL( papszEllipsoid[0], "WGS84"))
2114
2362
                    iEllipsoid = 8L;
2115
2363
            }
2116
 
                            
 
2364
 
2117
2365
            double adfProjParms[15];
2118
2366
            char **papszParms = (pszParms) ?
2119
2367
                CSLTokenizeString2( pszParms, ",",
2144
2392
        {
2145
2393
            double  dfCenterX, dfCenterY;
2146
2394
            int     iZone;
2147
 
                            
2148
 
            ReadCoordinates( CSLFetchNameValue( 
 
2395
 
 
2396
            ReadCoordinates( CSLFetchNameValue(
2149
2397
                                 papszGlobalMetadata, "SCENECENTER" ),
2150
2398
                             &dfCenterY, &dfCenterX );
2151
 
                            
 
2399
 
2152
2400
            // Calculate UTM zone from scene center coordinates
2153
2401
            iZone = 30 + (int) ((dfCenterX + 6.0) / 6.0);
2154
 
           
 
2402
 
2155
2403
            // Create projection definition
2156
2404
            if( dfCenterY > 0 )
2157
2405
                oSRS.SetUTM( iZone, TRUE );
2163
2411
        }
2164
2412
 
2165
2413
        // MODIS L1B
2166
 
        else if ( eProduct == PROD_MODIS_L1B )
 
2414
        else if ( eProduct == PROD_MODIS_L1B
 
2415
                  || eProduct == PROD_MODIS_L2 )
2167
2416
        {
2168
2417
            pszGCPProjection = CPLStrdup( SRS_WKT_WGS84 );
2169
2418
        }
2180
2429
        pasGCPList = (GDAL_GCP *) CPLCalloc( nGCPCount, sizeof( GDAL_GCP ) );
2181
2430
        GDALInitGCPs( nGCPCount, pasGCPList );
2182
2431
 
2183
 
        int iGCP = 0; 
 
2432
        int iGCP = 0;
2184
2433
        for ( i = 0; i < nYPoints; i += iGCPStepY )
2185
2434
        {
2186
2435
            for ( j = 0; j < nXPoints; j += iGCPStepX )
2187
2436
            {
2188
2437
                int iGeoOff =  i * nXPoints + j;
2189
 
 
 
2438
 
2190
2439
                pasGCPList[iGCP].dfGCPX =
2191
2440
                    AnyTypeToDouble(iWrkNumType,
2192
2441
                                    (void *)((char*)pLong+ iGeoOff*iDataSize));
2193
2442
                pasGCPList[iGCP].dfGCPY =
2194
2443
                    AnyTypeToDouble(iWrkNumType,
2195
2444
                                    (void *)((char*)pLat + iGeoOff*iDataSize));
2196
 
                                
 
2445
 
2197
2446
                // GCPs in Level 1A/1B dataset are in geocentric
2198
2447
                // coordinates. Convert them in geodetic (we
2199
2448
                // will convert latitudes only, longitudes
2203
2452
                if ( eProduct == PROD_ASTER_L1A
2204
2453
                     || eProduct == PROD_ASTER_L1B )
2205
2454
                {
2206
 
                    pasGCPList[iGCP].dfGCPY = 
 
2455
                    pasGCPList[iGCP].dfGCPY =
2207
2456
                        atan(tan(pasGCPList[iGCP].dfGCPY
2208
2457
                                 *PI/180)/0.99330562)*180/PI;
2209
2458
                }
2243
2492
/*      Establish geolocation metadata, but only if there is no         */
2244
2493
/*      lattice.  The lattice destroys the regularity of the grid.      */
2245
2494
/* -------------------------------------------------------------------- */
2246
 
    if( pLatticeX == NULL 
2247
 
        && iLatDim != -1 && iLongDim != -1 
 
2495
    if( pLatticeX == NULL
 
2496
        && iLatDim != -1 && iLongDim != -1
2248
2497
        && iPixelDim != -1 && iLineDim != -1 )
2249
2498
    {
2250
2499
        CPLString  osWrk;
2251
2500
 
2252
2501
        SetMetadataItem( "SRS", pszGCPProjection, "GEOLOCATION" );
2253
 
        
2254
 
        osWrk.Printf( "HDF4_EOS:EOS_SWATH_GEOL:\"%s\":%s:%s", 
 
2502
 
 
2503
        osWrk.Printf( "HDF4_EOS:EOS_SWATH_GEOL:\"%s\":%s:%s",
2255
2504
                      pszFilename, pszSubdatasetName,
2256
2505
                      papszGeolocations[iLongDim] );
2257
2506
        SetMetadataItem( "X_DATASET", osWrk, "GEOLOCATION" );
2258
2507
        SetMetadataItem( "X_BAND", "1" , "GEOLOCATION" );
2259
2508
 
2260
 
        osWrk.Printf( "HDF4_EOS:EOS_SWATH_GEOL:\"%s\":%s:%s", 
 
2509
        osWrk.Printf( "HDF4_EOS:EOS_SWATH_GEOL:\"%s\":%s:%s",
2261
2510
                      pszFilename, pszSubdatasetName,
2262
2511
                      papszGeolocations[iLatDim] );
2263
2512
        SetMetadataItem( "Y_DATASET", osWrk, "GEOLOCATION" );
2276
2525
            SetMetadataItem( "LINE_STEP", osWrk, "GEOLOCATION" );
2277
2526
        }
2278
2527
    }
2279
 
                    
 
2528
 
2280
2529
/* -------------------------------------------------------------------- */
2281
2530
/*      Cleanup                                                         */
2282
2531
/* -------------------------------------------------------------------- */
2289
2538
    CPLFree( paiIncrement );
2290
2539
    CPLFree( paiNumType );
2291
2540
    CPLFree( paiRank );
2292
 
    
 
2541
 
2293
2542
    CSLDestroy( papszGeolocations );
2294
2543
    CPLFree( pszGeoList );
2295
2544
 
2298
2547
        CPLFree( pszGCPProjection );
2299
2548
        pszGCPProjection = NULL;
2300
2549
    }
2301
 
    
 
2550
 
2302
2551
    return TRUE;
2303
2552
}
2304
2553
 
2309
2558
GDALDataset *HDF4ImageDataset::Open( GDALOpenInfo * poOpenInfo )
2310
2559
{
2311
2560
    int     i;
2312
 
    
 
2561
 
2313
2562
    if( !EQUALN( poOpenInfo->pszFilename, "HDF4_SDS:", 9 ) &&
2314
2563
        !EQUALN( poOpenInfo->pszFilename, "HDF4_GR:", 8 ) &&
2315
2564
        !EQUALN( poOpenInfo->pszFilename, "HDF4_GD:", 8 ) &&
2323
2572
    HDF4ImageDataset    *poDS;
2324
2573
 
2325
2574
    poDS = new HDF4ImageDataset( );
2326
 
 
2327
2575
    poDS->fp = poOpenInfo->fp;
2328
2576
    poOpenInfo->fp = NULL;
2329
2577
 
2368
2616
        poDS->iDatasetType = HDF4_EOS;
2369
2617
    else
2370
2618
        poDS->iDatasetType = HDF4_UNKNOWN;
2371
 
    
 
2619
 
2372
2620
    if( EQUAL( papszSubdatasetName[1], "GDAL_HDF4" ) )
2373
2621
        poDS->iSubdatasetType = H4ST_GDAL;
2374
2622
    else if( EQUAL( papszSubdatasetName[1], "EOS_GRID" ) )
2404
2652
        osSubdatasetName += ":";
2405
2653
        osSubdatasetName += papszSubdatasetName[4];
2406
2654
    }
2407
 
    
 
2655
 
2408
2656
/* -------------------------------------------------------------------- */
2409
2657
/*      Try opening the dataset.                                        */
2410
2658
/* -------------------------------------------------------------------- */
2411
2659
    int32       iAttribute, nValues, iAttrNumType;
2412
 
    double      dfNoData = 0.0;
2413
 
    int         bNoDataSet = FALSE;
2414
 
    
 
2660
    double      dfNoData, dfScale, dfOffset;
 
2661
    int         bNoDataSet = FALSE, bHaveScale = FALSE, bHaveOffset = FALSE;
 
2662
    const char  *pszUnits = NULL, *pszDescription = NULL;
 
2663
 
2415
2664
/* -------------------------------------------------------------------- */
2416
2665
/*      Select SDS or GR to read from.                                  */
2417
2666
/* -------------------------------------------------------------------- */
2444
2693
            case H4ST_EOS_SWATH:
2445
2694
            case H4ST_EOS_SWATH_GEOL:
2446
2695
            {
2447
 
                int32   hSW, nStrBufSize;
 
2696
                int32   hSW;
2448
2697
                char    *pszDimList = NULL;
2449
 
                    
 
2698
 
2450
2699
                if( poOpenInfo->eAccess == GA_ReadOnly )
2451
2700
                    poDS->hHDF4 = SWopen( poDS->pszFilename, DFACC_READ );
2452
2701
                else
2453
2702
                    poDS->hHDF4 = SWopen( poDS->pszFilename, DFACC_WRITE );
2454
 
                    
 
2703
 
2455
2704
                if( poDS->hHDF4 <= 0 )
2456
2705
                {
2457
2706
                    CPLDebug( "HDF4Image", "Can't open HDF4 file %s",
2472
2721
/* -------------------------------------------------------------------- */
2473
2722
/*      Decode the dimension map.                                       */
2474
2723
/* -------------------------------------------------------------------- */
2475
 
                SWnentries( hSW, HDFE_NENTDIM, &nStrBufSize );
 
2724
                int32   nStrBufSize = 0;
 
2725
 
 
2726
                if ( SWnentries( hSW, HDFE_NENTDIM, &nStrBufSize ) < 0
 
2727
                     || nStrBufSize <= 0 )
 
2728
                {
 
2729
                    CPLDebug( "HDF4Image",
 
2730
                              "Can't read a number of dimension maps." );
 
2731
                    delete poDS;
 
2732
                    return NULL;
 
2733
                }
 
2734
 
2476
2735
                pszDimList = (char *)CPLMalloc( nStrBufSize + 1 );
2477
 
                SWfieldinfo( hSW, poDS->pszFieldName, &poDS->iRank,
2478
 
                             poDS->aiDimSizes, &poDS->iNumType, pszDimList );
 
2736
                if ( SWfieldinfo( hSW, poDS->pszFieldName, &poDS->iRank,
 
2737
                                  poDS->aiDimSizes, &poDS->iNumType,
 
2738
                                  pszDimList ) < 0 )
 
2739
                {
 
2740
                    CPLDebug( "HDF4Image", "Can't read dimension maps." );
 
2741
                    CPLFree( pszDimList );
 
2742
                    delete poDS;
 
2743
                    return NULL;
 
2744
                }
 
2745
                pszDimList[nStrBufSize] = '\0';
2479
2746
#ifdef DEBUG
2480
2747
                CPLDebug( "HDF4Image",
2481
 
                          "List of dimensions in swath %s: %s",
 
2748
                          "List of dimensions in swath \"%s\": %s",
2482
2749
                          poDS->pszFieldName, pszDimList );
2483
2750
#endif
2484
2751
                poDS->GetImageDimensions( pszDimList );
2530
2797
                                            CSLT_HONOURSTRINGS );
2531
2798
                    if( !poDS->ProcessSwathGeolocation( hSW, papszDimList ) )
2532
2799
                    {
2533
 
                        CPLDebug( "HDF4Image", 
 
2800
                        CPLDebug( "HDF4Image",
2534
2801
                                  "No geolocation available for this swath." );
2535
2802
                    }
2536
2803
                    CSLDestroy( papszDimList );
2553
2820
                int32   nXSize, nYSize;
2554
2821
                char    szDimList[8192];
2555
2822
                double  adfUpLeft[2], adfLowRight[2], adfProjParms[15];
2556
 
                    
 
2823
 
2557
2824
                if( poOpenInfo->eAccess == GA_ReadOnly )
2558
2825
                    poDS->hHDF4 = GDopen( poDS->pszFilename, DFACC_READ );
2559
2826
                else
2560
2827
                    poDS->hHDF4 = GDopen( poDS->pszFilename, DFACC_WRITE );
2561
 
                    
 
2828
 
2562
2829
                if( poDS->hHDF4 <= 0 )
2563
2830
                {
2564
2831
                    delete poDS;
2579
2846
#endif
2580
2847
                poDS->GetImageDimensions( szDimList );
2581
2848
 
 
2849
                int32 tilecode, tilerank;
 
2850
                if( GDtileinfo( hGD , poDS->pszFieldName , &tilecode, &tilerank, NULL ) == 0 )
 
2851
                {
 
2852
                    if ( tilecode == HDFE_TILE )
 
2853
                    {
 
2854
                        int32 *tiledims = NULL;
 
2855
                        tiledims = (int32 *) CPLCalloc( tilerank , sizeof( int32 ) );
 
2856
                        GDtileinfo( hGD , poDS->pszFieldName , &tilecode, &tilerank, tiledims );
 
2857
                        if ( ( tilerank == 2 ) && ( poDS->iRank == tilerank  ) )
 
2858
                        {
 
2859
                            poDS->nBlockPreferredXSize = tiledims[1];
 
2860
                            poDS->nBlockPreferredYSize = tiledims[0];
 
2861
                            poDS->bReadTile = true;
 
2862
#ifdef DEBUG
 
2863
                            CPLDebug( "HDF4_EOS:EOS_GRID:","tilerank in grid %s: %d",
 
2864
                                      poDS->pszFieldName , (int)tilerank );
 
2865
                            CPLDebug( "HDF4_EOS:EOS_GRID:","tiledimens in grid %s: %d,%d",
 
2866
                                      poDS->pszFieldName , (int)tiledims[0] , (int)tiledims[1] );
 
2867
#endif
 
2868
                        }
 
2869
#ifdef DEBUG
 
2870
                        else
 
2871
                        {
 
2872
                                CPLDebug( "HDF4_EOS:EOS_GRID:","tilerank in grid %s: %d not supported",
 
2873
                                          poDS->pszFieldName , (int)tilerank );
 
2874
                        }
 
2875
#endif
 
2876
                        CPLFree(tiledims);
 
2877
                    }
 
2878
                    else
 
2879
                    {
 
2880
#ifdef DEBUG
 
2881
                        CPLDebug( "HDF4_EOS:EOS_GRID:","tilecode == HDFE_NOTILE in grid %s: %d",
 
2882
                                poDS->pszFieldName ,
 
2883
                                (int)poDS->iRank );
 
2884
#endif
 
2885
                    }
 
2886
                }
 
2887
#ifdef DEBUG
 
2888
                else
 
2889
                {
 
2890
                    CPLDebug( "HDF4_EOS:EOS_GRID:","ERROR GDtileinfo %s", poDS->pszFieldName );
 
2891
                }
 
2892
#endif
 
2893
 
2582
2894
/* -------------------------------------------------------------------- */
2583
2895
/*      Fetch projection information                                    */
2584
2896
/* -------------------------------------------------------------------- */
2595
2907
#endif
2596
2908
                    poDS->oSRS.importFromUSGS( iProjCode, iZoneCode,
2597
2909
                                               adfProjParms, iSphereCode );
2598
 
                        
 
2910
 
2599
2911
                    if ( poDS->pszProjection )
2600
2912
                        CPLFree( poDS->pszProjection );
2601
2913
                    poDS->oSRS.exportToWkt( &poDS->pszProjection );
2679
2991
                CPLFree( pNoDataValue );
2680
2992
            }
2681
2993
            break;
2682
 
                
 
2994
 
2683
2995
            default:
2684
2996
              break;
2685
2997
          }
 
2998
 
 
2999
/* -------------------------------------------------------------------- */
 
3000
/*      Fetch unit type, scale, offset and description                  */
 
3001
/*      Should be similar among various HDF-EOS kinds.                  */
 
3002
/* -------------------------------------------------------------------- */
 
3003
          {
 
3004
              const char *pszTmp =
 
3005
                  CSLFetchNameValue( poDS->papszLocalMetadata,
 
3006
                                     "scale_factor" );
 
3007
              if ( pszTmp )
 
3008
              {
 
3009
                  dfScale = CPLAtof( pszTmp );
 
3010
                  bHaveScale = TRUE;
 
3011
              }
 
3012
 
 
3013
              pszTmp =
 
3014
                  CSLFetchNameValue( poDS->papszLocalMetadata, "add_offset" );
 
3015
              if ( pszTmp )
 
3016
              {
 
3017
                  dfOffset = CPLAtof( pszTmp );
 
3018
                  bHaveOffset = TRUE;
 
3019
              }
 
3020
 
 
3021
              pszUnits = CSLFetchNameValue( poDS->papszLocalMetadata,
 
3022
                                            "units" );
 
3023
              pszDescription = CSLFetchNameValue( poDS->papszLocalMetadata,
 
3024
                                            "long_name" );
 
3025
          }
2686
3026
      }
2687
3027
      break;
2688
3028
 
2697
3037
              poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_READ, 0 );
2698
3038
          else
2699
3039
              poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_WRITE, 0 );
2700
 
            
 
3040
 
2701
3041
          if( poDS->hHDF4 <= 0 )
2702
3042
          {
2703
3043
              delete poDS;
2710
3050
              delete poDS;
2711
3051
              return NULL;
2712
3052
          }
2713
 
                
 
3053
 
2714
3054
          if ( poDS->ReadGlobalAttributes( poDS->hSD ) != CE_None )
2715
3055
          {
2716
3056
              delete poDS;
2728
3068
          if (poDS->iDataset < 0 || poDS->iDataset >= nDatasets)
2729
3069
          {
2730
3070
              CPLError(CE_Failure, CPLE_AppDefined,
2731
 
                       "Subdataset index should be between 0 and %d", nDatasets - 1);
 
3071
                       "Subdataset index should be between 0 and %ld",
 
3072
                       (long int)nDatasets - 1);
2732
3073
              delete poDS;
2733
3074
              return NULL;
2734
3075
          }
2824
3165
          else
2825
3166
            poDS->nRasterYSize = 1;
2826
3167
 
2827
 
          // Special case projection info for NRL generated files. 
 
3168
          // Special case projection info for NRL generated files.
2828
3169
          const char *pszMapProjectionSystem =
2829
 
              CSLFetchNameValue(poDS->papszGlobalMetadata, 
 
3170
              CSLFetchNameValue(poDS->papszGlobalMetadata,
2830
3171
                                "mapProjectionSystem");
2831
 
          if( pszMapProjectionSystem != NULL 
 
3172
          if( pszMapProjectionSystem != NULL
2832
3173
              && EQUAL(pszMapProjectionSystem,"NRL(USGS)") )
2833
3174
          {
2834
3175
              poDS->CaptureNRLGeoTransform();
2835
3176
          }
2836
3177
 
2837
 
          // Special case for coastwatch hdf files. 
2838
 
          if( CSLFetchNameValue( poDS->papszGlobalMetadata, 
 
3178
          // Special case for coastwatch hdf files.
 
3179
          if( CSLFetchNameValue( poDS->papszGlobalMetadata,
2839
3180
                                 "gctp_sys" ) != NULL )
2840
3181
              poDS->CaptureCoastwatchGCTPInfo();
2841
3182
 
2842
3183
          // Special case for MODIS geolocation
2843
3184
          poDS->ProcessModisSDSGeolocation();
 
3185
 
 
3186
          // Special case for NASA/CCRS Landsat in HDF
 
3187
          poDS->CaptureL1GMTLInfo();
2844
3188
      }
2845
3189
      break;
2846
3190
 
2852
3196
            poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_READ, 0 );
2853
3197
        else
2854
3198
            poDS->hHDF4 = Hopen( poDS->pszFilename, DFACC_WRITE, 0 );
2855
 
            
 
3199
 
2856
3200
        if( poDS->hHDF4 <= 0 )
2857
3201
        {
2858
3202
            delete poDS;
2865
3209
            delete poDS;
2866
3210
            return NULL;
2867
3211
        }
2868
 
                
 
3212
 
2869
3213
        poDS->iGR = GRselect( poDS->hGR, poDS->iDataset );
2870
3214
        if ( GRgetiminfo( poDS->iGR, poDS->szName,
2871
3215
                          &poDS->iRank, &poDS->iNumType,
2884
3228
            char    szAttrName[H4_MAX_NC_NAME];
2885
3229
            GRattrinfo( poDS->iGR, iAttribute, szAttrName,
2886
3230
                        &iAttrNumType, &nValues );
2887
 
            poDS->papszLocalMetadata = 
 
3231
            poDS->papszLocalMetadata =
2888
3232
                poDS->TranslateHDF4Attributes( poDS->iGR, iAttribute,
2889
3233
                                               szAttrName, iAttrNumType,
2890
3234
                                               nValues,
2893
3237
        poDS->SetMetadata( poDS->papszLocalMetadata, "" );
2894
3238
        // Read colour table
2895
3239
        GDALColorEntry oEntry;
2896
 
                 
 
3240
 
2897
3241
        poDS->iPal = GRgetlutid ( poDS->iGR, poDS->iDataset );
2898
3242
        if ( poDS->iPal != -1 )
2899
3243
        {
2907
3251
                oEntry.c2 = poDS->aiPaletteData[i][1];
2908
3252
                oEntry.c3 = poDS->aiPaletteData[i][2];
2909
3253
                oEntry.c4 = 255;
2910
 
                        
 
3254
 
2911
3255
                poDS->poColorTable->SetColorEntry( i, &oEntry );
2912
3256
            }
2913
3257
        }
2920
3264
        delete poDS;
2921
3265
        return NULL;
2922
3266
    }
2923
 
    
 
3267
 
2924
3268
    poDS->nRasterXSize = poDS->aiDimSizes[poDS->iXDim];
2925
3269
    if (poDS->iYDim >= 0)
2926
3270
        poDS->nRasterYSize = poDS->aiDimSizes[poDS->iYDim];
2949
3293
/* -------------------------------------------------------------------- */
2950
3294
    for( i = 1; i <= poDS->nBandCount; i++ )
2951
3295
    {
2952
 
        HDF4ImageRasterBand *poBand = 
 
3296
        HDF4ImageRasterBand *poBand =
2953
3297
            new HDF4ImageRasterBand( poDS, i,
2954
3298
                                     poDS->GetDataType( poDS->iNumType ) );
2955
3299
        poDS->SetBand( i, poBand );
2956
3300
 
2957
3301
        if ( bNoDataSet )
2958
3302
            poBand->SetNoDataValue( dfNoData );
 
3303
        if ( bHaveScale )
 
3304
        {
 
3305
            poBand->bHaveScale = TRUE;
 
3306
            poBand->dfScale = dfScale;
 
3307
        }
 
3308
        if ( bHaveOffset )
 
3309
        {
 
3310
            poBand->bHaveOffset = TRUE;
 
3311
            poBand->dfOffset = dfOffset;
 
3312
        }
 
3313
        if ( pszUnits )
 
3314
            poBand->osUnitType =  pszUnits;
 
3315
        if ( pszDescription )
 
3316
            poBand->SetDescription( pszDescription );
2959
3317
    }
2960
3318
 
2961
3319
/* -------------------------------------------------------------------- */
2990
3348
                                             "TransformationMatrix")) )
2991
3349
          {
2992
3350
              int i = 0;
2993
 
              char *pszString = (char *) pszValue; 
 
3351
              char *pszString = (char *) pszValue;
2994
3352
              while ( *pszValue && i < 6 )
2995
3353
              {
2996
3354
                  poDS->adfGeoTransform[i++] = strtod(pszString, &pszString);
3032
3390
 
3033
3391
          // Read coordinate system and geotransform matrix
3034
3392
          poDS->oSRS.SetWellKnownGeogCS( "WGS84" );
3035
 
            
 
3393
 
3036
3394
          if ( EQUAL(CSLFetchNameValue(poDS->papszGlobalMetadata,
3037
3395
                                       "Map Projection"),
3038
3396
                     "Equidistant Cylindrical") )
3080
3438
              for( i = 1; i <= poDS->nBands; i++ )
3081
3439
              {
3082
3440
                  poDS->GetRasterBand(i)->SetNoDataValue(
3083
 
                      CPLAtof( CSLFetchNameValue(poDS->papszLocalMetadata, 
 
3441
                      CPLAtof( CSLFetchNameValue(poDS->papszLocalMetadata,
3084
3442
                                                 "missing_value") ) );
3085
3443
              }
3086
3444
          }
3087
3445
 
3088
3446
          // Coastwatch offset and scale.
3089
 
          if( CSLFetchNameValue( poDS->papszLocalMetadata, "scale_factor" ) 
 
3447
          if( CSLFetchNameValue( poDS->papszLocalMetadata, "scale_factor" )
3090
3448
              && CSLFetchNameValue( poDS->papszLocalMetadata, "add_offset" ) )
3091
3449
          {
3092
3450
              for( i = 1; i <= poDS->nBands; i++ )
3093
3451
              {
3094
 
                  HDF4ImageRasterBand *poBand = 
 
3452
                  HDF4ImageRasterBand *poBand =
3095
3453
                      (HDF4ImageRasterBand *) poDS->GetRasterBand(i);
3096
3454
 
3097
 
                  poBand->bHaveScaleAndOffset = TRUE;
3098
 
                  poBand->dfScale = 
3099
 
                      CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata, 
 
3455
                  poBand->bHaveScale = poBand->bHaveOffset = TRUE;
 
3456
                  poBand->dfScale =
 
3457
                      CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3100
3458
                                                  "scale_factor" ) );
3101
 
                  poBand->dfOffset = -1 * poBand->dfScale * 
3102
 
                      CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata, 
 
3459
                  poBand->dfOffset = -1 * poBand->dfScale *
 
3460
                      CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3103
3461
                                                  "add_offset" ) );
3104
3462
              }
3105
3463
          }
3107
3465
          // this is a modis level3 convention (data from ACT)
3108
3466
          // Eg data/hdf/act/modis/MODAM2004280160000.L3_NOAA_GMX
3109
3467
 
3110
 
          if( CSLFetchNameValue( poDS->papszLocalMetadata, 
3111
 
                                 "scalingSlope" ) 
3112
 
              && CSLFetchNameValue( poDS->papszLocalMetadata, 
 
3468
          if( CSLFetchNameValue( poDS->papszLocalMetadata,
 
3469
                                 "scalingSlope" )
 
3470
              && CSLFetchNameValue( poDS->papszLocalMetadata,
3113
3471
                                    "scalingIntercept" ) )
3114
3472
          {
3115
3473
              int i;
3116
3474
              CPLString osUnits;
3117
 
              
3118
 
              if( CSLFetchNameValue( poDS->papszLocalMetadata, 
 
3475
 
 
3476
              if( CSLFetchNameValue( poDS->papszLocalMetadata,
3119
3477
                                     "productUnits" ) )
3120
3478
              {
3121
 
                  osUnits = CSLFetchNameValue( poDS->papszLocalMetadata, 
 
3479
                  osUnits = CSLFetchNameValue( poDS->papszLocalMetadata,
3122
3480
                                               "productUnits" );
3123
3481
              }
3124
3482
 
3125
3483
              for( i = 1; i <= poDS->nBands; i++ )
3126
3484
              {
3127
 
                  HDF4ImageRasterBand *poBand = 
 
3485
                  HDF4ImageRasterBand *poBand =
3128
3486
                      (HDF4ImageRasterBand *) poDS->GetRasterBand(i);
3129
3487
 
3130
 
                  poBand->bHaveScaleAndOffset = TRUE;
3131
 
                  poBand->dfScale = 
3132
 
                      CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata, 
 
3488
                  poBand->bHaveScale = poBand->bHaveOffset = TRUE;
 
3489
                  poBand->dfScale =
 
3490
                      CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3133
3491
                                                  "scalingSlope" ) );
3134
 
                  poBand->dfOffset = 
3135
 
                      CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata, 
 
3492
                  poBand->dfOffset =
 
3493
                      CPLAtof( CSLFetchNameValue( poDS->papszLocalMetadata,
3136
3494
                                                  "scalingIntercept" ) );
3137
3495
 
3138
3496
                  poBand->osUnitType = osUnits;
3208
3566
    {
3209
3567
        CPLError( CE_Failure, CPLE_OpenFailed,
3210
3568
                  "Can't create HDF4 file %s", pszFilename );
 
3569
        delete poDS;
3211
3570
        return NULL;
3212
3571
    }
3213
3572
    poDS->iXDim = 1;