← Back to branch summary

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

« back to all changes in this revision

Viewing changes to frmts/tsx/tsxdataset.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/tsx/tsxdataset.cpp
1
1
/******************************************************************************
2
 
 * $Id: tsxdataset.cpp 17664 2009-09-21 21:16:45Z rouault $
 
2
 * $Id: tsxdataset.cpp 23569 2011-12-13 21:53:07Z rouault $
3
3
 *
4
4
 * Project:     TerraSAR-X XML Product Support
5
5
 * Purpose:     Support for TerraSAR-X XML Metadata files
31
31
 
32
32
#include "gdal_pam.h"
33
33
#include "cpl_minixml.h"
34
 
 
35
 
CPL_CVSID("$Id: tsxdataset.cpp 17664 2009-09-21 21:16:45Z rouault $");
 
34
#include "ogr_spatialref.h"
 
35
 
 
36
#define MAX_GCPS 5000    //this should be more than enough ground control points
 
37
 
 
38
CPL_CVSID("$Id: tsxdataset.cpp 23569 2011-12-13 21:53:07Z rouault $");
36
39
 
37
40
CPL_C_START
38
41
void GDALRegister_TSX(void);
40
43
 
41
44
 
42
45
enum ePolarization {
43
 
        HH=0,
44
 
        HV,
45
 
        VH,
46
 
        VV
 
46
    HH=0,
 
47
    HV,
 
48
    VH,
 
49
    VV
47
50
};
48
51
 
49
52
enum eProductType {
61
64
/* GetFilePath: return a relative path to a file within an XML node.
62
65
 * Returns Null on failure
63
66
 */
64
 
const char *GetFilePath(CPLXMLNode *psXMLNode, char **pszNodeType) {
65
 
        const char *pszDirectory, *pszFilename;
66
 
 
67
 
        pszDirectory = CPLGetXMLValue( psXMLNode, "file.location.path", "" );
68
 
        pszFilename = CPLGetXMLValue( psXMLNode, "file.location.filename", "" );
69
 
        *pszNodeType = strdup(CPLGetXMLValue (psXMLNode, "type", " " ));
70
 
 
71
 
        if (pszDirectory == NULL || pszFilename == NULL) {
72
 
                return NULL;
73
 
        }
74
 
 
75
 
        return strdup( CPLFormFilename( pszDirectory, pszFilename, "" ) );
 
67
const char *GetFilePath(CPLXMLNode *psXMLNode, const char **pszNodeType) {
 
68
    const char *pszDirectory, *pszFilename;
 
69
 
 
70
    pszDirectory = CPLGetXMLValue( psXMLNode, "file.location.path", "" );
 
71
    pszFilename = CPLGetXMLValue( psXMLNode, "file.location.filename", "" );
 
72
    *pszNodeType = CPLGetXMLValue (psXMLNode, "type", " " );
 
73
 
 
74
    if (pszDirectory == NULL || pszFilename == NULL) {
 
75
        return NULL;
 
76
    }
 
77
 
 
78
    return CPLFormFilename( pszDirectory, pszFilename, "" );
76
79
}
77
80
 
78
81
/************************************************************************/
79
82
/* ==================================================================== */
80
 
/*                                  TSXDataset                                          */
 
83
/*                                TSXDataset                                 */
81
84
/* ==================================================================== */
82
85
/************************************************************************/
83
86
 
87
90
 
88
91
    char *pszGCPProjection;
89
92
 
90
 
        char *pszGeorefFile;
91
 
        FILE *fp;
 
93
    char *pszProjection;
 
94
    double adfGeoTransform[6];
 
95
    bool bHaveGeoTransform;
 
96
 
 
97
    char *pszGeorefFile;
92
98
 
93
99
    eProductType nProduct;
94
100
public:
95
 
        TSXDataset();
96
 
        ~TSXDataset();
97
 
    
 
101
    TSXDataset();
 
102
    ~TSXDataset();
 
103
 
98
104
    virtual int GetGCPCount();
99
105
    virtual const char *GetGCPProjection();
100
106
    virtual const GDAL_GCP *GetGCPs();
101
107
 
 
108
    CPLErr GetGeoTransform( double* padfTransform);
 
109
    const char* GetProjectionRef();
 
110
 
102
111
    static GDALDataset *Open( GDALOpenInfo *poOpenInfo );
103
 
        static int Identify( GDALOpenInfo *poOpenInfo );
 
112
    static int Identify( GDALOpenInfo *poOpenInfo );
 
113
private:
 
114
    bool getGCPsFromGEOREF_XML(char *pszGeorefFilename);
104
115
};
105
116
 
106
117
/************************************************************************/
107
118
/* ==================================================================== */
108
 
/*                                          TSXRasterBand                           */
 
119
/*                                TSXRasterBand                           */
109
120
/* ==================================================================== */
110
121
/************************************************************************/
111
122
 
112
123
class TSXRasterBand : public GDALPamRasterBand {
113
 
        GDALDataset *poBand;
114
 
        ePolarization ePol;
 
124
    GDALDataset *poBand;
 
125
    ePolarization ePol;
115
126
public:
116
 
        TSXRasterBand( TSXDataset *poDSIn, GDALDataType eDataType, 
117
 
                ePolarization ePol, GDALDataset *poBand );
118
 
        virtual ~TSXRasterBand();
119
 
 
120
 
        virtual CPLErr IReadBlock( int nBlockXOff, int nBlockYOff, void *pImage );
121
 
 
122
 
        static GDALDataset *Open( GDALOpenInfo *poOpenInfo );
 
127
    TSXRasterBand( TSXDataset *poDSIn, GDALDataType eDataType,
 
128
        ePolarization ePol, GDALDataset *poBand );
 
129
    virtual ~TSXRasterBand();
 
130
 
 
131
    virtual CPLErr IReadBlock( int nBlockXOff, int nBlockYOff, void *pImage );
 
132
 
 
133
    static GDALDataset *Open( GDALOpenInfo *poOpenInfo );
123
134
};
124
135
 
125
136
/************************************************************************/
133
144
    this->eDataType = eDataType;
134
145
    this->ePol = ePol;
135
146
 
136
 
        switch (ePol) {
137
 
                case HH:
138
 
                SetMetadataItem( "POLARIMETRIC_INTERP", "HH" );
139
 
                        break;
140
 
                case HV:
141
 
                SetMetadataItem( "POLARIMETRIC_INTERP", "HV" );
142
 
                        break;
143
 
                case VH:
144
 
                SetMetadataItem( "POLARIMETRIC_INTERP", "VH" );
145
 
                        break;
146
 
                case VV:
147
 
                SetMetadataItem( "POLARIMETRIC_INTERP", "VV" );
148
 
                        break;
149
 
        }
150
 
 
151
 
 
152
 
        /* now setup the actual raster reader */
153
 
        this->poBand = poBand;
154
 
 
155
 
        GDALRasterBand *poSrcBand = poBand->GetRasterBand( 1 );
156
 
        poSrcBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
 
147
    switch (ePol) {
 
148
        case HH:
 
149
            SetMetadataItem( "POLARIMETRIC_INTERP", "HH" );
 
150
            break;
 
151
        case HV:
 
152
            SetMetadataItem( "POLARIMETRIC_INTERP", "HV" );
 
153
            break;
 
154
        case VH:
 
155
            SetMetadataItem( "POLARIMETRIC_INTERP", "VH" );
 
156
            break;
 
157
        case VV:
 
158
            SetMetadataItem( "POLARIMETRIC_INTERP", "VV" );
 
159
            break;
 
160
    }
 
161
 
 
162
 
 
163
    /* now setup the actual raster reader */
 
164
    this->poBand = poBand;
 
165
 
 
166
    GDALRasterBand *poSrcBand = poBand->GetRasterBand( 1 );
 
167
    poSrcBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
157
168
}
158
169
 
159
170
/************************************************************************/
172
183
CPLErr TSXRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
173
184
                                  void * pImage )
174
185
{
175
 
        int nRequestYSize;
176
 
 
177
 
        /* Check if the last strip is partial so we can avoid over-requesting */
178
 
        if ( (nBlockYOff + 1) * nBlockYSize > nRasterYSize ) {
179
 
                nRequestYSize = nRasterYSize - nBlockYOff * nBlockYSize;
180
 
                memset( pImage, 0, (GDALGetDataTypeSize( eDataType ) / 8) *
181
 
                        nBlockXSize * nBlockYSize);
182
 
        }
183
 
        else {
184
 
                nRequestYSize = nBlockYSize;
185
 
        }
186
 
 
187
 
        /* Read Complex Data */
188
 
        if ( eDataType == GDT_CInt16 ) {
189
 
                return poBand->RasterIO( GF_Read, nBlockXOff * nBlockXSize,
190
 
                        nBlockYOff * nBlockYSize, nBlockXSize, nRequestYSize,
191
 
                        pImage, nBlockXSize, nRequestYSize, GDT_CInt16, 1, NULL, 4,
192
 
                        nBlockXSize * 4, 0 );
193
 
        }
194
 
        else { /* Detected Product */
195
 
                return poBand->RasterIO( GF_Read, nBlockXOff * nBlockXSize,
196
 
                        nBlockYOff * nBlockYSize, nBlockXSize, nRequestYSize,
197
 
                        pImage, nBlockXSize, nRequestYSize, GDT_UInt16, 1, NULL, 2,
198
 
                        nBlockXSize * 2, 0 );
199
 
        }
 
186
    int nRequestYSize;
 
187
 
 
188
    /* Check if the last strip is partial so we can avoid over-requesting */
 
189
    if ( (nBlockYOff + 1) * nBlockYSize > nRasterYSize ) {
 
190
        nRequestYSize = nRasterYSize - nBlockYOff * nBlockYSize;
 
191
        memset( pImage, 0, (GDALGetDataTypeSize( eDataType ) / 8) *
 
192
            nBlockXSize * nBlockYSize);
 
193
    }
 
194
    else {
 
195
        nRequestYSize = nBlockYSize;
 
196
    }
 
197
 
 
198
    /* Read Complex Data */
 
199
    if ( eDataType == GDT_CInt16 ) {
 
200
        return poBand->RasterIO( GF_Read, nBlockXOff * nBlockXSize,
 
201
            nBlockYOff * nBlockYSize, nBlockXSize, nRequestYSize,
 
202
            pImage, nBlockXSize, nRequestYSize, GDT_CInt16, 1, NULL, 4,
 
203
            nBlockXSize * 4, 0 );
 
204
    }
 
205
    else { /* Detected Product */
 
206
        return poBand->RasterIO( GF_Read, nBlockXOff * nBlockXSize,
 
207
            nBlockYOff * nBlockYSize, nBlockXSize, nRequestYSize,
 
208
            pImage, nBlockXSize, nRequestYSize, GDT_UInt16, 1, NULL, 2,
 
209
            nBlockXSize * 2, 0 );
 
210
    }
200
211
}
201
212
 
202
213
/************************************************************************/
203
214
/* ==================================================================== */
204
 
/*                                              TSXDataset                                              */
 
215
/*                                TSXDataset                                */
205
216
/* ==================================================================== */
206
217
/************************************************************************/
207
218
 
213
224
    nGCPCount = 0;
214
225
    pasGCPList = NULL;
215
226
    pszGCPProjection = CPLStrdup("");
 
227
    pszProjection = CPLStrdup("");
 
228
    adfGeoTransform[0] = 0.0;
 
229
    adfGeoTransform[1] = 1.0;
 
230
    adfGeoTransform[2] = 0.0;
 
231
    adfGeoTransform[3] = 0.0;
 
232
    adfGeoTransform[4] = 0.0;
 
233
    adfGeoTransform[5] = 1.0;
 
234
    bHaveGeoTransform = FALSE;
216
235
}
217
236
 
218
237
/************************************************************************/
221
240
 
222
241
TSXDataset::~TSXDataset() {
223
242
    FlushCache();
 
243
 
 
244
    CPLFree( pszProjection );
 
245
 
 
246
    CPLFree( pszGCPProjection );
 
247
    if( nGCPCount > 0 )
 
248
    {
 
249
        GDALDeinitGCPs( nGCPCount, pasGCPList );
 
250
        CPLFree( pasGCPList );
 
251
    }
224
252
}
225
253
 
226
254
/************************************************************************/
227
255
/*                              Identify()                              */
228
256
/************************************************************************/
229
257
 
230
 
int TSXDataset::Identify( GDALOpenInfo *poOpenInfo ) {
231
 
        if (poOpenInfo->fp == NULL || poOpenInfo->nHeaderBytes < 260) 
232
 
                return 0;
233
 
 
234
 
        /* Check if the filename contains TSX1_SAR */
235
 
        if (!EQUALN(CPLGetBasename( poOpenInfo->pszFilename ), "TSX1_SAR", 8))
236
 
                return 0;
237
 
 
238
 
        /* finally look for the <level1Product tag */
239
 
        if (!EQUALN((char *)poOpenInfo->pabyHeader, "<level1Product", 14)) 
240
 
                return 0;
241
 
 
242
 
        return 1;
 
258
int TSXDataset::Identify( GDALOpenInfo *poOpenInfo )
 
259
{
 
260
    if (poOpenInfo->fp == NULL || poOpenInfo->nHeaderBytes < 260)
 
261
    {
 
262
        if( poOpenInfo->bIsDirectory )
 
263
        {
 
264
            CPLString osFilename =
 
265
                CPLFormCIFilename( poOpenInfo->pszFilename, CPLGetFilename( poOpenInfo->pszFilename ), "xml" );
 
266
 
 
267
            /* Check if the filename contains TSX1_SAR (TerraSAR-X) or TDX1_SAR (TanDEM-X) */
 
268
            if (!(EQUALN(CPLGetBasename( osFilename ), "TSX1_SAR", 8) ||
 
269
                  EQUALN(CPLGetBasename( osFilename ), "TDX1_SAR", 8)))
 
270
                return 0;
 
271
 
 
272
            VSIStatBufL sStat;
 
273
            if( VSIStatL( osFilename, &sStat ) == 0 )
 
274
                return 1;
 
275
 
 
276
            return 0;
 
277
        }
 
278
    }
 
279
 
 
280
    /* Check if the filename contains TSX1_SAR (TerraSAR-X) or TDX1_SAR (TanDEM-X) */
 
281
    if (!(EQUALN(CPLGetBasename( poOpenInfo->pszFilename ), "TSX1_SAR", 8) ||
 
282
          EQUALN(CPLGetBasename( poOpenInfo->pszFilename ), "TDX1_SAR", 8)))
 
283
        return 0;
 
284
 
 
285
    /* finally look for the <level1Product tag */
 
286
    if (!EQUALN((char *)poOpenInfo->pabyHeader, "<level1Product", 14))
 
287
        return 0;
 
288
 
 
289
    return 1;
 
290
}
 
291
 
 
292
/************************************************************************/
 
293
/*                                getGCPsFromGEOREF_XML()               */
 
294
/*Reads georeferencing information from the TerraSAR-X GEOREF.XML file    */
 
295
/*and writes the information to the dataset's gcp list and projection     */
 
296
/*string.                                                                */
 
297
/*Returns true on success.                                                */
 
298
/************************************************************************/
 
299
bool TSXDataset::getGCPsFromGEOREF_XML(char *pszGeorefFilename)
 
300
{
 
301
    //open GEOREF.xml
 
302
    CPLXMLNode *psGeorefData;
 
303
    psGeorefData = CPLParseXMLFile( pszGeorefFilename );
 
304
    if (psGeorefData==NULL)
 
305
        return false;
 
306
 
 
307
    //get the ellipsoid and semi-major, semi-minor axes
 
308
    CPLXMLNode *psSphere;
 
309
    const char *pszEllipsoidName;
 
310
    double minor_axis, major_axis, inv_flattening;
 
311
    OGRSpatialReference osr;
 
312
    psSphere = CPLGetXMLNode( psGeorefData, "=geoReference.referenceFrames.sphere" );
 
313
    if (psSphere!=NULL)
 
314
    {
 
315
        pszEllipsoidName = CPLGetXMLValue( psSphere, "ellipsoidID", "" );
 
316
        minor_axis = atof(CPLGetXMLValue( psSphere, "semiMinorAxis", "0.0" ));
 
317
        major_axis = atof(CPLGetXMLValue( psSphere, "semiMajorAxis", "0.0" ));
 
318
        //save datum parameters to the spatial reference
 
319
        if ( EQUAL(pszEllipsoidName, "") || minor_axis==0.0 || major_axis==0.0 )
 
320
        {
 
321
            CPLError(CE_Warning,CPLE_AppDefined,"Warning- incomplete"
 
322
                " ellipsoid information.  Using wgs-84 parameters.\n");
 
323
            osr.SetWellKnownGeogCS( "WGS84" );
 
324
        }
 
325
        else if ( EQUAL( pszEllipsoidName, "WGS84" ) ) {
 
326
            osr.SetWellKnownGeogCS( "WGS84" );
 
327
        }
 
328
        else {
 
329
            inv_flattening = major_axis/(major_axis - minor_axis);
 
330
            osr.SetGeogCS( "","",pszEllipsoidName, major_axis, inv_flattening);
 
331
        }
 
332
    }
 
333
 
 
334
    //get gcps
 
335
    CPLXMLNode *psNode;
 
336
    CPLXMLNode *psGeolocationGrid = CPLGetXMLNode( psGeorefData, "=geoReference.geolocationGrid" );
 
337
    if (psGeolocationGrid==NULL)
 
338
    {
 
339
        CPLDestroyXMLNode( psGeorefData );
 
340
        return false;
 
341
    }
 
342
    nGCPCount = atoi(CPLGetXMLValue( psGeolocationGrid, "numberOfGridPoints.total", "0" ));
 
343
    //count the gcps if the given count value is invalid
 
344
    if (nGCPCount<=0)
 
345
    {
 
346
        for( psNode = psGeolocationGrid->psChild; psNode != NULL; psNode = psNode->psNext )
 
347
            if( EQUAL(psNode->pszValue,"gridPoint") )
 
348
                nGCPCount++ ;
 
349
    }
 
350
    //if there are no gcps, fail
 
351
    if(nGCPCount<=0)
 
352
    {
 
353
        CPLDestroyXMLNode( psGeorefData );
 
354
        return false;
 
355
    }
 
356
 
 
357
    //put some reasonable limits of the number of gcps
 
358
    if (nGCPCount>MAX_GCPS )
 
359
        nGCPCount=MAX_GCPS;
 
360
    //allocate memory for the gcps
 
361
    pasGCPList = (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP),nGCPCount);
 
362
    //loop through all gcps and set info
 
363
    int gcps_allocated = nGCPCount;    //save the number allocated to ensure it does not run off the end of the array
 
364
    nGCPCount=0;    //reset to zero and count
 
365
    //do a check on the grid point to make sure it has lat,long row, and column
 
366
    //it seems that only SSC products contain row, col - how to map lat long otherwise??
 
367
    //for now fail if row and col are not present - just check the first and assume the rest are the same
 
368
    for( psNode = psGeolocationGrid->psChild; psNode != NULL; psNode = psNode->psNext )
 
369
    {
 
370
         if( !EQUAL(psNode->pszValue,"gridPoint") )
 
371
             continue;
 
372
 
 
373
         if (    !strcmp(CPLGetXMLValue(psNode,"col","error"), "error") ||
 
374
                 !strcmp(CPLGetXMLValue(psNode,"row","error"), "error") ||
 
375
                 !strcmp(CPLGetXMLValue(psNode,"lon","error"), "error") ||
 
376
                 !strcmp(CPLGetXMLValue(psNode,"lat","error"), "error"))
 
377
        {
 
378
            CPLDestroyXMLNode( psGeorefData );
 
379
            return false;
 
380
        }
 
381
    }
 
382
    for( psNode = psGeolocationGrid->psChild; psNode != NULL; psNode = psNode->psNext )
 
383
    {
 
384
        //break out if the end of the array has been reached
 
385
        if (nGCPCount >= gcps_allocated)
 
386
        {
 
387
            CPLError(CE_Warning, CPLE_AppDefined, "GDAL TSX driver: Truncating the number of GCPs.");
 
388
            break;
 
389
        }
 
390
 
 
391
         char    szID[32];
 
392
         GDAL_GCP   *psGCP = pasGCPList + nGCPCount;
 
393
 
 
394
         if( !EQUAL(psNode->pszValue,"gridPoint") )
 
395
             continue;
 
396
 
 
397
         nGCPCount++ ;
 
398
 
 
399
         sprintf( szID, "%d", nGCPCount );
 
400
         psGCP->pszId = CPLStrdup( szID );
 
401
         psGCP->pszInfo = CPLStrdup("");
 
402
         psGCP->dfGCPPixel =
 
403
             atof(CPLGetXMLValue(psNode,"col","0"));
 
404
         psGCP->dfGCPLine =
 
405
             atof(CPLGetXMLValue(psNode,"row","0"));
 
406
         psGCP->dfGCPX =
 
407
             atof(CPLGetXMLValue(psNode,"lon",""));
 
408
         psGCP->dfGCPY =
 
409
             atof(CPLGetXMLValue(psNode,"lat",""));
 
410
         //looks like height is in meters - should it be converted so xyz are all on the same scale??
 
411
         psGCP->dfGCPZ = 0;
 
412
             //atof(CPLGetXMLValue(psNode,"height",""));
 
413
    }
 
414
 
 
415
    CPLFree(pszGCPProjection);
 
416
    osr.exportToWkt( &(pszGCPProjection) );
 
417
 
 
418
    CPLDestroyXMLNode( psGeorefData );
 
419
 
 
420
    return true;
243
421
}
244
422
 
245
423
/************************************************************************/
250
428
/* -------------------------------------------------------------------- */
251
429
/*      Is this a TerraSAR-X product file?                              */
252
430
/* -------------------------------------------------------------------- */
253
 
        if (!TSXDataset::Identify( poOpenInfo )) {
254
 
                return NULL; /* nope */
255
 
        }
256
 
    
 
431
    if (!TSXDataset::Identify( poOpenInfo ))
 
432
    {
 
433
        return NULL; /* nope */
 
434
    }
 
435
 
257
436
/* -------------------------------------------------------------------- */
258
437
/*      Confirm the requested access is supported.                      */
259
438
/* -------------------------------------------------------------------- */
260
439
    if( poOpenInfo->eAccess == GA_Update )
261
440
    {
262
 
        CPLError( CE_Failure, CPLE_NotSupported, 
 
441
        CPLError( CE_Failure, CPLE_NotSupported,
263
442
                  "The TSX driver does not support update access to existing"
264
443
                  " datasets.\n" );
265
444
        return NULL;
266
445
    }
267
 
    
268
 
        /* Ingest the XML */
269
 
        CPLXMLNode *psData, *psComponents, *psProductInfo;
270
 
        psData = CPLParseXMLFile( poOpenInfo->pszFilename );
271
 
 
272
 
        /* find the product components */
273
 
        psComponents = CPLGetXMLNode( psData, "=level1Product.productComponents" );
274
 
        if (psComponents == NULL) {
275
 
                CPLError( CE_Failure, CPLE_OpenFailed, 
276
 
                        "Unable to find <productComponents> tag in file.\n" );
277
 
                return NULL;
278
 
        }
279
 
 
280
 
        /* find the product info tag */
281
 
        psProductInfo = CPLGetXMLNode( psData, "=level1Product.productInfo" );
282
 
        if (psComponents == NULL) {
283
 
                CPLError( CE_Failure, CPLE_OpenFailed,
284
 
                        "Unable to find <productInfo> tag in file.\n" );
285
 
                return NULL;
286
 
        }
 
446
 
 
447
    CPLString osFilename;
 
448
 
 
449
    if( poOpenInfo->bIsDirectory )
 
450
    {
 
451
        osFilename =
 
452
                CPLFormCIFilename( poOpenInfo->pszFilename, CPLGetFilename( poOpenInfo->pszFilename ), "xml" );
 
453
    }
 
454
    else
 
455
        osFilename = poOpenInfo->pszFilename;
 
456
 
 
457
    /* Ingest the XML */
 
458
    CPLXMLNode *psData, *psComponents, *psProductInfo;
 
459
    psData = CPLParseXMLFile( osFilename );
 
460
    if (psData == NULL)
 
461
        return NULL;
 
462
 
 
463
    /* find the product components */
 
464
    psComponents = CPLGetXMLNode( psData, "=level1Product.productComponents" );
 
465
    if (psComponents == NULL) {
 
466
        CPLError( CE_Failure, CPLE_OpenFailed,
 
467
            "Unable to find <productComponents> tag in file.\n" );
 
468
        CPLDestroyXMLNode(psData);
 
469
        return NULL;
 
470
    }
 
471
 
 
472
    /* find the product info tag */
 
473
    psProductInfo = CPLGetXMLNode( psData, "=level1Product.productInfo" );
 
474
    if (psProductInfo == NULL) {
 
475
        CPLError( CE_Failure, CPLE_OpenFailed,
 
476
            "Unable to find <productInfo> tag in file.\n" );
 
477
        CPLDestroyXMLNode(psData);
 
478
        return NULL;
 
479
    }
287
480
 
288
481
/* -------------------------------------------------------------------- */
289
482
/*      Create the dataset.                                             */
290
483
/* -------------------------------------------------------------------- */
291
 
        
 
484
 
292
485
    TSXDataset *poDS = new TSXDataset();
293
 
        poDS->fp = poOpenInfo->fp;
294
 
        poOpenInfo->fp = NULL;
295
486
 
296
487
/* -------------------------------------------------------------------- */
297
488
/*      Read in product info.                                           */
299
490
 
300
491
    poDS->SetMetadataItem( "SCENE_CENTRE_TIME", CPLGetXMLValue( psProductInfo,
301
492
        "sceneInfo.sceneCenterCoord.azimuthTimeUTC", "unknown" ) );
302
 
        poDS->SetMetadataItem( "OPERATIONAL_MODE", CPLGetXMLValue( psProductInfo, 
303
 
                "generationInfo.groundOperationsType", "unknown" ) );
304
 
        poDS->SetMetadataItem( "ORBIT_CYCLE", CPLGetXMLValue( psProductInfo,
305
 
                "missionInfo.orbitCycle", "unknown" ) );
306
 
        poDS->SetMetadataItem( "ABSOLUTE_ORBIT", CPLGetXMLValue( psProductInfo,
307
 
                "missionInfo.absOrbit", "unknown" ) );
308
 
        poDS->SetMetadataItem( "ORBIT_DIRECTION", CPLGetXMLValue( psProductInfo,
309
 
                "missionInfo.orbitDirection", "unknown" ) );
310
 
        poDS->SetMetadataItem( "IMAGING_MODE", CPLGetXMLValue( psProductInfo,
311
 
                "acquisitionInfo.imagingMode", "unknown" ) );
312
 
        poDS->SetMetadataItem( "PRODUCT_VARIANT", CPLGetXMLValue( psProductInfo,
313
 
                "productVariantInfo.productVariant", "unknown" ) ); 
314
 
        char *pszDataType = strdup( CPLGetXMLValue( psProductInfo,
315
 
                "imageDataInfo.imageDataType", "unknown" ) );
316
 
        poDS->SetMetadataItem( "IMAGE_TYPE", pszDataType ); 
317
 
        
318
 
        /* Get raster information */
319
 
        int nRows = atoi( CPLGetXMLValue( psProductInfo,
320
 
                "imageDataInfo.imageRaster.numberOfRows", "" ) );
321
 
        int nCols = atoi( CPLGetXMLValue( psProductInfo,
322
 
                "imageDataInfo.imageRaster.numberOfColumns", "" ) );
323
 
 
324
 
        poDS->nRasterXSize = nCols;
325
 
        poDS->nRasterYSize = nRows;
326
 
 
327
 
        poDS->SetMetadataItem( "ROW_SPACING", CPLGetXMLValue( psProductInfo,
328
 
                "imageDataInfo.imageRaster.rowSpacing", "unknown" ) );
329
 
        poDS->SetMetadataItem( "COL_SPACING", CPLGetXMLValue( psProductInfo,
330
 
                "imageDataInfo.imageRaster.columnSpacing", "unknown" ) );
 
493
    poDS->SetMetadataItem( "OPERATIONAL_MODE", CPLGetXMLValue( psProductInfo,
 
494
        "generationInfo.groundOperationsType", "unknown" ) );
 
495
    poDS->SetMetadataItem( "ORBIT_CYCLE", CPLGetXMLValue( psProductInfo,
 
496
        "missionInfo.orbitCycle", "unknown" ) );
 
497
    poDS->SetMetadataItem( "ABSOLUTE_ORBIT", CPLGetXMLValue( psProductInfo,
 
498
        "missionInfo.absOrbit", "unknown" ) );
 
499
    poDS->SetMetadataItem( "ORBIT_DIRECTION", CPLGetXMLValue( psProductInfo,
 
500
        "missionInfo.orbitDirection", "unknown" ) );
 
501
    poDS->SetMetadataItem( "IMAGING_MODE", CPLGetXMLValue( psProductInfo,
 
502
        "acquisitionInfo.imagingMode", "unknown" ) );
 
503
    poDS->SetMetadataItem( "PRODUCT_VARIANT", CPLGetXMLValue( psProductInfo,
 
504
        "productVariantInfo.productVariant", "unknown" ) );
 
505
    char *pszDataType = CPLStrdup( CPLGetXMLValue( psProductInfo,
 
506
        "imageDataInfo.imageDataType", "unknown" ) );
 
507
    poDS->SetMetadataItem( "IMAGE_TYPE", pszDataType );
 
508
 
 
509
    /* Get raster information */
 
510
    int nRows = atoi( CPLGetXMLValue( psProductInfo,
 
511
        "imageDataInfo.imageRaster.numberOfRows", "" ) );
 
512
    int nCols = atoi( CPLGetXMLValue( psProductInfo,
 
513
        "imageDataInfo.imageRaster.numberOfColumns", "" ) );
 
514
 
 
515
    poDS->nRasterXSize = nCols;
 
516
    poDS->nRasterYSize = nRows;
 
517
 
 
518
    poDS->SetMetadataItem( "ROW_SPACING", CPLGetXMLValue( psProductInfo,
 
519
        "imageDataInfo.imageRaster.rowSpacing", "unknown" ) );
 
520
    poDS->SetMetadataItem( "COL_SPACING", CPLGetXMLValue( psProductInfo,
 
521
        "imageDataInfo.imageRaster.columnSpacing", "unknown" ) );
331
522
    poDS->SetMetadataItem( "COL_SPACING_UNITS", CPLGetXMLValue( psProductInfo,
332
523
        "imageDataInfo.imageRaster.columnSpacing.units", "unknown" ) );
333
524
 
334
 
        /* Get equivalent number of looks */
335
 
        poDS->SetMetadataItem( "AZIMUTH_LOOKS", CPLGetXMLValue( psProductInfo,
336
 
                "imageDataInfo.imageRaster.azimuthLooks", "unknown" ) );
337
 
        poDS->SetMetadataItem( "RANGE_LOOKS", CPLGetXMLValue( psProductInfo,
338
 
                "imageDataInfo.imageRaster.rangeLooks", "unknown" ) );
 
525
    /* Get equivalent number of looks */
 
526
    poDS->SetMetadataItem( "AZIMUTH_LOOKS", CPLGetXMLValue( psProductInfo,
 
527
        "imageDataInfo.imageRaster.azimuthLooks", "unknown" ) );
 
528
    poDS->SetMetadataItem( "RANGE_LOOKS", CPLGetXMLValue( psProductInfo,
 
529
        "imageDataInfo.imageRaster.rangeLooks", "unknown" ) );
339
530
 
340
531
    const char *pszProductVariant;
341
 
    pszProductVariant = CPLGetXMLValue( psProductInfo, 
 
532
    pszProductVariant = CPLGetXMLValue( psProductInfo,
342
533
        "productVariantInfo.productVariant", "unknown" );
343
534
 
344
535
    poDS->SetMetadataItem( "PRODUCT_VARIANT", pszProductVariant );
355
546
    else
356
547
        poDS->nProduct = eUnknown;
357
548
 
358
 
        /* Start reading in the product components */
359
 
        const char *pszPath;
360
 
        char *pszGeorefFile = NULL;
361
 
        CPLXMLNode *psComponent;
362
 
        for (psComponent = psComponents->psChild; psComponent != NULL;
363
 
                 psComponent = psComponent->psNext)
364
 
        {
365
 
                char *pszType;
366
 
                pszPath = CPLFormFilename( 
367
 
                                CPLGetDirname( poOpenInfo->pszFilename ),
368
 
                                GetFilePath(psComponent, &pszType), 
369
 
                                "" );
370
 
                const char *pszPolLayer = CPLGetXMLValue(psComponent, "polLayer", " ");
371
 
 
372
 
                if ( !EQUALN(pszType," ",1) ) {
373
 
                        if (EQUALN(pszType, "MAPPING_GRID", 12) ) {
374
 
                                /* the mapping grid... save as a metadata item this path */
375
 
                                poDS->SetMetadataItem( "MAPPING_GRID", pszPath );       
376
 
                        }
377
 
                        else if (EQUALN(pszType, "GEOREF", 6)) {
378
 
                                /* save the path to the georef data for later use */
379
 
                                pszGeorefFile = strdup( pszPath );
380
 
                        }
381
 
                        CPLFree(pszType);
382
 
                }
383
 
                else if( !EQUALN(pszPolLayer, " ", 1) && 
384
 
                        EQUALN(psComponent->pszValue, "imageData", 9) ) {
385
 
                        /* determine the polarization of this band */
386
 
                        ePolarization ePol;
387
 
                        if ( EQUALN(pszPolLayer, "HH", 2) ) {
388
 
                                ePol = HH;
389
 
                        }
390
 
                        else if ( EQUALN(pszPolLayer, "HV" , 2) ) {
391
 
                                ePol = HV;
392
 
                        }
393
 
                        else if ( EQUALN(pszPolLayer, "VH", 2) ) {
394
 
                                ePol = VH;
395
 
                        }
396
 
                        else {
397
 
                                ePol = VV;
398
 
                        }
399
 
 
400
 
                        GDALDataType eDataType = EQUALN(pszDataType, "COMPLEX", 7) ?
401
 
                                GDT_CInt16 : GDT_UInt16;
402
 
 
403
 
                        /* try opening the file that represents that band */
404
 
                        TSXRasterBand *poBand;
405
 
                        GDALDataset *poBandData;
406
 
 
407
 
                        poBandData = (GDALDataset *) GDALOpen( pszPath, GA_ReadOnly );
408
 
                        if ( poBandData != NULL ) {
409
 
                                poBand = new TSXRasterBand( poDS, eDataType, ePol, 
410
 
                                        poBandData );
411
 
                                poDS->SetBand( poDS->GetRasterCount() + 1, poBand );
412
 
                        }
413
 
                }
414
 
        }
415
 
 
416
 
        CPLFree(pszDataType);
 
549
    /* Start reading in the product components */
 
550
    const char *pszPath;
 
551
    char *pszGeorefFile = NULL;
 
552
    CPLXMLNode *psComponent;
 
553
    CPLErr geoTransformErr=CE_Failure;
 
554
    for (psComponent = psComponents->psChild; psComponent != NULL;
 
555
         psComponent = psComponent->psNext)
 
556
    {
 
557
        const char *pszType = NULL;
 
558
        pszPath = CPLFormFilename(
 
559
                CPLGetDirname( osFilename ),
 
560
                GetFilePath(psComponent, &pszType),
 
561
                "" );
 
562
        const char *pszPolLayer = CPLGetXMLValue(psComponent, "polLayer", " ");
 
563
 
 
564
        if ( !EQUALN(pszType," ",1) ) {
 
565
            if (EQUALN(pszType, "MAPPING_GRID", 12) ) {
 
566
                /* the mapping grid... save as a metadata item this path */
 
567
                poDS->SetMetadataItem( "MAPPING_GRID", pszPath );
 
568
            }
 
569
            else if (EQUALN(pszType, "GEOREF", 6)) {
 
570
                /* save the path to the georef data for later use */
 
571
                pszGeorefFile = CPLStrdup( pszPath );
 
572
            }
 
573
        }
 
574
        else if( !EQUALN(pszPolLayer, " ", 1) &&
 
575
            EQUALN(psComponent->pszValue, "imageData", 9) ) {
 
576
            /* determine the polarization of this band */
 
577
            ePolarization ePol;
 
578
            if ( EQUALN(pszPolLayer, "HH", 2) ) {
 
579
                ePol = HH;
 
580
            }
 
581
            else if ( EQUALN(pszPolLayer, "HV" , 2) ) {
 
582
                ePol = HV;
 
583
            }
 
584
            else if ( EQUALN(pszPolLayer, "VH", 2) ) {
 
585
                ePol = VH;
 
586
            }
 
587
            else {
 
588
                ePol = VV;
 
589
            }
 
590
 
 
591
            GDALDataType eDataType = EQUALN(pszDataType, "COMPLEX", 7) ?
 
592
                GDT_CInt16 : GDT_UInt16;
 
593
 
 
594
            /* try opening the file that represents that band */
 
595
            TSXRasterBand *poBand;
 
596
            GDALDataset *poBandData;
 
597
 
 
598
            poBandData = (GDALDataset *) GDALOpen( pszPath, GA_ReadOnly );
 
599
            if ( poBandData != NULL ) {
 
600
                poBand = new TSXRasterBand( poDS, eDataType, ePol,
 
601
                    poBandData );
 
602
                poDS->SetBand( poDS->GetRasterCount() + 1, poBand );
 
603
 
 
604
                //copy georeferencing info from the band
 
605
                //need error checking??
 
606
                //it will just save the info from the last band
 
607
                CPLFree( poDS->pszProjection );
 
608
                poDS->pszProjection = CPLStrdup(poBandData->GetProjectionRef());
 
609
                geoTransformErr = poBandData->GetGeoTransform(poDS->adfGeoTransform);
 
610
            }
 
611
        }
 
612
    }
 
613
 
 
614
    //now check if there is a geotransform
 
615
    if ( strcmp(poDS->pszProjection, "") && geoTransformErr==CE_None)
 
616
    {
 
617
        poDS->bHaveGeoTransform = TRUE;
 
618
    }
 
619
    else
 
620
    {
 
621
        poDS->bHaveGeoTransform = FALSE;
 
622
        CPLFree( poDS->pszProjection );
 
623
        poDS->pszProjection = CPLStrdup("");
 
624
        poDS->adfGeoTransform[0] = 0.0;
 
625
        poDS->adfGeoTransform[1] = 1.0;
 
626
        poDS->adfGeoTransform[2] = 0.0;
 
627
        poDS->adfGeoTransform[3] = 0.0;
 
628
        poDS->adfGeoTransform[4] = 0.0;
 
629
        poDS->adfGeoTransform[5] = 1.0;
 
630
    }
 
631
 
 
632
    CPLFree(pszDataType);
417
633
 
418
634
 
419
635
/* -------------------------------------------------------------------- */
420
636
/*      Check and set matrix representation.                            */
421
637
/* -------------------------------------------------------------------- */
422
638
 
423
 
        if (poDS->GetRasterCount() == 4) {
424
 
                poDS->SetMetadataItem( "MATRIX_REPRESENTATION", "SCATTERING" );
425
 
        }
 
639
    if (poDS->GetRasterCount() == 4) {
 
640
        poDS->SetMetadataItem( "MATRIX_REPRESENTATION", "SCATTERING" );
 
641
    }
426
642
 
427
643
/* -------------------------------------------------------------------- */
428
644
/*      Read the four corners and centre GCPs in                        */
429
645
/* -------------------------------------------------------------------- */
430
646
 
431
 
    CPLXMLNode *psSceneInfo = CPLGetXMLNode( psData, 
 
647
    CPLXMLNode *psSceneInfo = CPLGetXMLNode( psData,
432
648
        "=level1Product.productInfo.sceneInfo" );
433
 
    /* for SSC products */
434
 
    if (poDS->nProduct == eSSC && psSceneInfo != NULL) {
435
 
        CPLXMLNode *psNode;
436
 
        int nGCP = 0;
437
 
        double dfAvgHeight = atof(CPLGetXMLValue(psSceneInfo, 
438
 
            "sceneAverageHeight", "0.0"));
439
 
        char szID[3];
440
 
 
441
 
        poDS->nGCPCount = 5; /* 5 GCPs provided */
442
 
        poDS->pasGCPList = (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP), 
443
 
            poDS->nGCPCount);
444
 
 
445
 
        /* iterate over GCPs */
446
 
        for (psNode = psSceneInfo->psChild; psNode != NULL; 
447
 
             psNode = psNode->psNext )
448
 
        {
449
 
            GDAL_GCP *psGCP = poDS->pasGCPList + nGCP;
450
 
 
451
 
            if (!EQUAL(psNode->pszValue, "sceneCenterCoord") && 
452
 
                !EQUAL(psNode->pszValue, "sceneCornerCoord"))
453
 
                continue;
454
 
 
455
 
            CPLSPrintf( szID, "%d", nGCP );
456
 
            
457
 
            psGCP->dfGCPPixel = atof(CPLGetXMLValue(psNode, "refColumn", 
458
 
                "0.0"));
459
 
            psGCP->dfGCPLine = atof(CPLGetXMLValue(psNode, "refRow", "0.0"));
460
 
            psGCP->dfGCPX = atof(CPLGetXMLValue(psNode, "lon", "0.0"));
461
 
            psGCP->dfGCPY = atof(CPLGetXMLValue(psNode, "lat", "0.0"));
462
 
            psGCP->dfGCPZ = dfAvgHeight;
463
 
            psGCP->pszId = CPLStrdup( szID );
464
 
            psGCP->pszInfo = CPLStrdup("");
465
 
 
466
 
            nGCP++;
467
 
        }
468
 
    }
469
 
    else if (psSceneInfo != NULL) {
 
649
    if (psSceneInfo != NULL)
 
650
    {
470
651
        /* extract the GCPs from the provided file */
471
 
 
472
 
        /* TODO */
 
652
        bool success = false;
 
653
        if (pszGeorefFile != NULL)
 
654
            success = poDS->getGCPsFromGEOREF_XML(pszGeorefFile);
 
655
 
 
656
        //if the gcp's cannot be extracted from the georef file, try to get the corner coordinates
 
657
        //for now just SSC because the others don't have refColumn and refRow
 
658
        if (!success && poDS->nProduct == eSSC)
 
659
        {
 
660
            CPLXMLNode *psNode;
 
661
            int nGCP = 0;
 
662
            double dfAvgHeight = atof(CPLGetXMLValue(psSceneInfo,
 
663
                "sceneAverageHeight", "0.0"));
 
664
 
 
665
            //count and allocate gcps - there should be five - 4 corners and a centre
 
666
            poDS->nGCPCount = 0;
 
667
            for (psNode = psSceneInfo->psChild; psNode != NULL; psNode = psNode->psNext )
 
668
            {
 
669
                if (!EQUAL(psNode->pszValue, "sceneCenterCoord") &&
 
670
                    !EQUAL(psNode->pszValue, "sceneCornerCoord"))
 
671
                    continue;
 
672
 
 
673
                poDS->nGCPCount++;
 
674
            }
 
675
            if (poDS->nGCPCount > 0)
 
676
            {
 
677
                poDS->pasGCPList = (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP), poDS->nGCPCount);
 
678
 
 
679
                /* iterate over GCPs */
 
680
                for (psNode = psSceneInfo->psChild; psNode != NULL; psNode = psNode->psNext )
 
681
                {
 
682
                    GDAL_GCP *psGCP = poDS->pasGCPList + nGCP;
 
683
 
 
684
                    if (!EQUAL(psNode->pszValue, "sceneCenterCoord") &&
 
685
                        !EQUAL(psNode->pszValue, "sceneCornerCoord"))
 
686
                        continue;
 
687
 
 
688
                    psGCP->dfGCPPixel = atof(CPLGetXMLValue(psNode, "refColumn",
 
689
                        "0.0"));
 
690
                    psGCP->dfGCPLine = atof(CPLGetXMLValue(psNode, "refRow", "0.0"));
 
691
                    psGCP->dfGCPX = atof(CPLGetXMLValue(psNode, "lon", "0.0"));
 
692
                    psGCP->dfGCPY = atof(CPLGetXMLValue(psNode, "lat", "0.0"));
 
693
                    psGCP->dfGCPZ = dfAvgHeight;
 
694
                    psGCP->pszId = CPLStrdup( CPLSPrintf( "%d", nGCP ) );
 
695
                    psGCP->pszInfo = CPLStrdup("");
 
696
 
 
697
                    nGCP++;
 
698
                }
 
699
 
 
700
                //set the projection string - the fields are lat/long - seems to be WGS84 datum
 
701
                OGRSpatialReference osr;
 
702
                osr.SetWellKnownGeogCS( "WGS84" );
 
703
                CPLFree(poDS->pszGCPProjection);
 
704
                osr.exportToWkt( &(poDS->pszGCPProjection) );
 
705
            }
 
706
        }
 
707
 
 
708
        //gcps override geotransform - does it make sense to have both??
 
709
        if (poDS->nGCPCount>0)
 
710
        {
 
711
            poDS->bHaveGeoTransform = FALSE;
 
712
            CPLFree( poDS->pszProjection );
 
713
            poDS->pszProjection = CPLStrdup("");
 
714
            poDS->adfGeoTransform[0] = 0.0;
 
715
            poDS->adfGeoTransform[1] = 1.0;
 
716
            poDS->adfGeoTransform[2] = 0.0;
 
717
            poDS->adfGeoTransform[3] = 0.0;
 
718
            poDS->adfGeoTransform[4] = 0.0;
 
719
            poDS->adfGeoTransform[5] = 1.0;
 
720
        }
 
721
 
473
722
    }
474
723
    else {
475
 
        CPLError(CE_Warning, CPLE_AppDefined, 
476
 
            "Unable to find sceneInfo tag in XML document. " 
 
724
        CPLError(CE_Warning, CPLE_AppDefined,
 
725
            "Unable to find sceneInfo tag in XML document. "
477
726
            "Proceeding with caution.");
478
727
    }
479
728
 
 
729
    CPLFree(pszGeorefFile);
 
730
 
480
731
/* -------------------------------------------------------------------- */
481
732
/*      Initialize any PAM information.                                 */
482
733
/* -------------------------------------------------------------------- */
518
769
}
519
770
 
520
771
/************************************************************************/
521
 
/*                         GDALRegister_TSX()                          */
 
772
/*                          GetProjectionRef()                          */
 
773
/************************************************************************/
 
774
const char *TSXDataset::GetProjectionRef()
 
775
{
 
776
    return pszProjection;
 
777
}
 
778
 
 
779
/************************************************************************/
 
780
/*                               GetGeotransform()                      */
 
781
/************************************************************************/
 
782
CPLErr TSXDataset::GetGeoTransform(double* padfTransform)
 
783
{
 
784
    memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
 
785
 
 
786
    if (bHaveGeoTransform)
 
787
        return( CE_None );
 
788
 
 
789
    return( CE_Failure );
 
790
}
 
791
 
 
792
/************************************************************************/
 
793
/*                         GDALRegister_TSX()                           */
522
794
/************************************************************************/
523
795
 
524
796
void GDALRegister_TSX() {
525
 
    GDALDriver  *poDriver;
 
797
    GDALDriver    *poDriver;
526
798
 
527
799
    if( GDALGetDriverByName( "TSX" ) == NULL )
528
800
    {
529
801
        poDriver = new GDALDriver();
530
 
        
 
802
 
531
803
        poDriver->SetDescription( "TSX" );
532
 
        poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, 
 
804
        poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
533
805
                                   "TerraSAR-X Product" );
534
806
/*        poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "frmt_tsx.html" ); */
535
807
 
 
808
        poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
 
809
 
536
810
        poDriver->pfnOpen = TSXDataset::Open;
537
 
                poDriver->pfnIdentify = TSXDataset::Identify;
 
811
        poDriver->pfnIdentify = TSXDataset::Identify;
538
812
 
539
813
        GetGDALDriverManager()->RegisterDriver( poDriver );
540
814
    }