1 : /******************************************************************************
2 : * $Id: tsxdataset.cpp 24254 2012-04-17 22:34:17Z rouault $
3 : *
4 : * Project: TerraSAR-X XML Product Support
5 : * Purpose: Support for TerraSAR-X XML Metadata files
6 : * Author: Philippe Vachon <philippe@cowpig.ca>
7 : * Description: This driver adds support for reading metadata and georef data
8 : * associated with TerraSAR-X products.
9 : *
10 : ******************************************************************************
11 : * Copyright (c) 2007, Philippe Vachon <philippe@cowpig.ca>
12 : *
13 : * Permission is hereby granted, free of charge, to any person obtaining a
14 : * copy of this software and associated documentation files (the "Software"),
15 : * to deal in the Software without restriction, including without limitation
16 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 : * and/or sell copies of the Software, and to permit persons to whom the
18 : * Software is furnished to do so, subject to the following conditions:
19 : *
20 : * The above copyright notice and this permission notice shall be included
21 : * in all copies or substantial portions of the Software.
22 : *
23 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 : * DEALINGS IN THE SOFTWARE.
30 : ****************************************************************************/
31 :
32 : #include "gdal_pam.h"
33 : #include "cpl_minixml.h"
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 24254 2012-04-17 22:34:17Z rouault $");
39 :
40 : CPL_C_START
41 : void GDALRegister_TSX(void);
42 : CPL_C_END
43 :
44 :
45 : enum ePolarization {
46 : HH=0,
47 : HV,
48 : VH,
49 : VV
50 : };
51 :
52 : enum eProductType {
53 : eSSC = 0,
54 : eMGD,
55 : eEEC,
56 : eGEC,
57 : eUnknown
58 : };
59 :
60 : /************************************************************************/
61 : /* Helper Functions */
62 : /************************************************************************/
63 :
64 : /* GetFilePath: return a relative path to a file within an XML node.
65 : * Returns Null on failure
66 : */
67 0 : const char *GetFilePath(CPLXMLNode *psXMLNode, const char **pszNodeType) {
68 : const char *pszDirectory, *pszFilename;
69 :
70 0 : pszDirectory = CPLGetXMLValue( psXMLNode, "file.location.path", "" );
71 0 : pszFilename = CPLGetXMLValue( psXMLNode, "file.location.filename", "" );
72 0 : *pszNodeType = CPLGetXMLValue (psXMLNode, "type", " " );
73 :
74 0 : if (pszDirectory == NULL || pszFilename == NULL) {
75 0 : return NULL;
76 : }
77 :
78 0 : return CPLFormFilename( pszDirectory, pszFilename, "" );
79 : }
80 :
81 : /************************************************************************/
82 : /* ==================================================================== */
83 : /* TSXDataset */
84 : /* ==================================================================== */
85 : /************************************************************************/
86 :
87 : class TSXDataset : public GDALPamDataset {
88 : int nGCPCount;
89 : GDAL_GCP *pasGCPList;
90 :
91 : char *pszGCPProjection;
92 :
93 : char *pszProjection;
94 : double adfGeoTransform[6];
95 : bool bHaveGeoTransform;
96 :
97 : char *pszGeorefFile;
98 :
99 : eProductType nProduct;
100 : public:
101 : TSXDataset();
102 : ~TSXDataset();
103 :
104 : virtual int GetGCPCount();
105 : virtual const char *GetGCPProjection();
106 : virtual const GDAL_GCP *GetGCPs();
107 :
108 : CPLErr GetGeoTransform( double* padfTransform);
109 : const char* GetProjectionRef();
110 :
111 : static GDALDataset *Open( GDALOpenInfo *poOpenInfo );
112 : static int Identify( GDALOpenInfo *poOpenInfo );
113 : private:
114 : bool getGCPsFromGEOREF_XML(char *pszGeorefFilename);
115 : };
116 :
117 : /************************************************************************/
118 : /* ==================================================================== */
119 : /* TSXRasterBand */
120 : /* ==================================================================== */
121 : /************************************************************************/
122 :
123 : class TSXRasterBand : public GDALPamRasterBand {
124 : GDALDataset *poBand;
125 : ePolarization ePol;
126 : public:
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 );
134 : };
135 :
136 : /************************************************************************/
137 : /* TSXRasterBand */
138 : /************************************************************************/
139 :
140 0 : TSXRasterBand::TSXRasterBand( TSXDataset *poDS, GDALDataType eDataType,
141 0 : ePolarization ePol, GDALDataset *poBand )
142 : {
143 0 : this->poDS = poDS;
144 0 : this->eDataType = eDataType;
145 0 : this->ePol = ePol;
146 :
147 0 : switch (ePol) {
148 : case HH:
149 0 : SetMetadataItem( "POLARIMETRIC_INTERP", "HH" );
150 0 : break;
151 : case HV:
152 0 : SetMetadataItem( "POLARIMETRIC_INTERP", "HV" );
153 0 : break;
154 : case VH:
155 0 : SetMetadataItem( "POLARIMETRIC_INTERP", "VH" );
156 0 : break;
157 : case VV:
158 0 : SetMetadataItem( "POLARIMETRIC_INTERP", "VV" );
159 : break;
160 : }
161 :
162 :
163 : /* now setup the actual raster reader */
164 0 : this->poBand = poBand;
165 :
166 0 : GDALRasterBand *poSrcBand = poBand->GetRasterBand( 1 );
167 0 : poSrcBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
168 0 : }
169 :
170 : /************************************************************************/
171 : /* TSXRasterBand() */
172 : /************************************************************************/
173 :
174 0 : TSXRasterBand::~TSXRasterBand() {
175 0 : if( poBand != NULL )
176 0 : GDALClose( (GDALRasterBandH) poBand );
177 0 : }
178 :
179 : /************************************************************************/
180 : /* IReadBlock() */
181 : /************************************************************************/
182 :
183 0 : CPLErr TSXRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
184 : void * pImage )
185 : {
186 : int nRequestYSize;
187 :
188 : /* Check if the last strip is partial so we can avoid over-requesting */
189 0 : if ( (nBlockYOff + 1) * nBlockYSize > nRasterYSize ) {
190 0 : nRequestYSize = nRasterYSize - nBlockYOff * nBlockYSize;
191 : memset( pImage, 0, (GDALGetDataTypeSize( eDataType ) / 8) *
192 0 : nBlockXSize * nBlockYSize);
193 : }
194 : else {
195 0 : nRequestYSize = nBlockYSize;
196 : }
197 :
198 : /* Read Complex Data */
199 0 : 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 0 : 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 0 : nBlockXSize * 2, 0 );
210 : }
211 : }
212 :
213 : /************************************************************************/
214 : /* ==================================================================== */
215 : /* TSXDataset */
216 : /* ==================================================================== */
217 : /************************************************************************/
218 :
219 : /************************************************************************/
220 : /* TSXDataset() */
221 : /************************************************************************/
222 :
223 0 : TSXDataset::TSXDataset() {
224 0 : nGCPCount = 0;
225 0 : pasGCPList = NULL;
226 0 : pszGCPProjection = CPLStrdup("");
227 0 : pszProjection = CPLStrdup("");
228 0 : adfGeoTransform[0] = 0.0;
229 0 : adfGeoTransform[1] = 1.0;
230 0 : adfGeoTransform[2] = 0.0;
231 0 : adfGeoTransform[3] = 0.0;
232 0 : adfGeoTransform[4] = 0.0;
233 0 : adfGeoTransform[5] = 1.0;
234 0 : bHaveGeoTransform = FALSE;
235 0 : }
236 :
237 : /************************************************************************/
238 : /* ~TSXDataset() */
239 : /************************************************************************/
240 :
241 0 : TSXDataset::~TSXDataset() {
242 0 : FlushCache();
243 :
244 0 : CPLFree( pszProjection );
245 :
246 0 : CPLFree( pszGCPProjection );
247 0 : if( nGCPCount > 0 )
248 : {
249 0 : GDALDeinitGCPs( nGCPCount, pasGCPList );
250 0 : CPLFree( pasGCPList );
251 : }
252 0 : }
253 :
254 : /************************************************************************/
255 : /* Identify() */
256 : /************************************************************************/
257 :
258 23374 : int TSXDataset::Identify( GDALOpenInfo *poOpenInfo )
259 : {
260 23374 : if (poOpenInfo->fp == NULL || poOpenInfo->nHeaderBytes < 260)
261 : {
262 22106 : if( poOpenInfo->bIsDirectory )
263 : {
264 : CPLString osFilename =
265 84 : CPLFormCIFilename( poOpenInfo->pszFilename, CPLGetFilename( poOpenInfo->pszFilename ), "xml" );
266 :
267 : /* Check if the filename contains TSX1_SAR (TerraSAR-X) or TDX1_SAR (TanDEM-X) */
268 84 : if (!(EQUALN(CPLGetBasename( osFilename ), "TSX1_SAR", 8) ||
269 : EQUALN(CPLGetBasename( osFilename ), "TDX1_SAR", 8)))
270 84 : return 0;
271 :
272 : VSIStatBufL sStat;
273 0 : if( VSIStatL( osFilename, &sStat ) == 0 )
274 0 : return 1;
275 : }
276 :
277 22022 : return 0;
278 : }
279 :
280 : /* Check if the filename contains TSX1_SAR (TerraSAR-X) or TDX1_SAR (TanDEM-X) */
281 1268 : if (!(EQUALN(CPLGetBasename( poOpenInfo->pszFilename ), "TSX1_SAR", 8) ||
282 : EQUALN(CPLGetBasename( poOpenInfo->pszFilename ), "TDX1_SAR", 8)))
283 1268 : return 0;
284 :
285 : /* finally look for the <level1Product tag */
286 0 : if (!EQUALN((char *)poOpenInfo->pabyHeader, "<level1Product", 14))
287 0 : return 0;
288 :
289 0 : 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 0 : bool TSXDataset::getGCPsFromGEOREF_XML(char *pszGeorefFilename)
300 : {
301 : //open GEOREF.xml
302 : CPLXMLNode *psGeorefData;
303 0 : psGeorefData = CPLParseXMLFile( pszGeorefFilename );
304 0 : if (psGeorefData==NULL)
305 0 : 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 0 : OGRSpatialReference osr;
312 0 : psSphere = CPLGetXMLNode( psGeorefData, "=geoReference.referenceFrames.sphere" );
313 0 : if (psSphere!=NULL)
314 : {
315 0 : pszEllipsoidName = CPLGetXMLValue( psSphere, "ellipsoidID", "" );
316 0 : minor_axis = atof(CPLGetXMLValue( psSphere, "semiMinorAxis", "0.0" ));
317 0 : major_axis = atof(CPLGetXMLValue( psSphere, "semiMajorAxis", "0.0" ));
318 : //save datum parameters to the spatial reference
319 0 : if ( EQUAL(pszEllipsoidName, "") || minor_axis==0.0 || major_axis==0.0 )
320 : {
321 : CPLError(CE_Warning,CPLE_AppDefined,"Warning- incomplete"
322 0 : " ellipsoid information. Using wgs-84 parameters.\n");
323 0 : osr.SetWellKnownGeogCS( "WGS84" );
324 : }
325 0 : else if ( EQUAL( pszEllipsoidName, "WGS84" ) ) {
326 0 : osr.SetWellKnownGeogCS( "WGS84" );
327 : }
328 : else {
329 0 : inv_flattening = major_axis/(major_axis - minor_axis);
330 0 : osr.SetGeogCS( "","",pszEllipsoidName, major_axis, inv_flattening);
331 : }
332 : }
333 :
334 : //get gcps
335 : CPLXMLNode *psNode;
336 0 : CPLXMLNode *psGeolocationGrid = CPLGetXMLNode( psGeorefData, "=geoReference.geolocationGrid" );
337 0 : if (psGeolocationGrid==NULL)
338 : {
339 0 : CPLDestroyXMLNode( psGeorefData );
340 0 : return false;
341 : }
342 0 : nGCPCount = atoi(CPLGetXMLValue( psGeolocationGrid, "numberOfGridPoints.total", "0" ));
343 : //count the gcps if the given count value is invalid
344 0 : if (nGCPCount<=0)
345 : {
346 0 : for( psNode = psGeolocationGrid->psChild; psNode != NULL; psNode = psNode->psNext )
347 0 : if( EQUAL(psNode->pszValue,"gridPoint") )
348 0 : nGCPCount++ ;
349 : }
350 : //if there are no gcps, fail
351 0 : if(nGCPCount<=0)
352 : {
353 0 : CPLDestroyXMLNode( psGeorefData );
354 0 : return false;
355 : }
356 :
357 : //put some reasonable limits of the number of gcps
358 0 : if (nGCPCount>MAX_GCPS )
359 0 : nGCPCount=MAX_GCPS;
360 : //allocate memory for the gcps
361 0 : pasGCPList = (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP),nGCPCount);
362 : //loop through all gcps and set info
363 0 : int gcps_allocated = nGCPCount; //save the number allocated to ensure it does not run off the end of the array
364 0 : 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 0 : for( psNode = psGeolocationGrid->psChild; psNode != NULL; psNode = psNode->psNext )
369 : {
370 0 : if( !EQUAL(psNode->pszValue,"gridPoint") )
371 0 : continue;
372 :
373 0 : 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 0 : CPLDestroyXMLNode( psGeorefData );
379 0 : return false;
380 : }
381 : }
382 0 : for( psNode = psGeolocationGrid->psChild; psNode != NULL; psNode = psNode->psNext )
383 : {
384 : //break out if the end of the array has been reached
385 0 : if (nGCPCount >= gcps_allocated)
386 : {
387 0 : CPLError(CE_Warning, CPLE_AppDefined, "GDAL TSX driver: Truncating the number of GCPs.");
388 0 : break;
389 : }
390 :
391 : char szID[32];
392 0 : GDAL_GCP *psGCP = pasGCPList + nGCPCount;
393 :
394 0 : if( !EQUAL(psNode->pszValue,"gridPoint") )
395 0 : continue;
396 :
397 0 : nGCPCount++ ;
398 :
399 0 : sprintf( szID, "%d", nGCPCount );
400 0 : psGCP->pszId = CPLStrdup( szID );
401 0 : psGCP->pszInfo = CPLStrdup("");
402 : psGCP->dfGCPPixel =
403 0 : atof(CPLGetXMLValue(psNode,"col","0"));
404 : psGCP->dfGCPLine =
405 0 : atof(CPLGetXMLValue(psNode,"row","0"));
406 : psGCP->dfGCPX =
407 0 : atof(CPLGetXMLValue(psNode,"lon",""));
408 : psGCP->dfGCPY =
409 0 : atof(CPLGetXMLValue(psNode,"lat",""));
410 : //looks like height is in meters - should it be converted so xyz are all on the same scale??
411 0 : psGCP->dfGCPZ = 0;
412 : //atof(CPLGetXMLValue(psNode,"height",""));
413 : }
414 :
415 0 : CPLFree(pszGCPProjection);
416 0 : osr.exportToWkt( &(pszGCPProjection) );
417 :
418 0 : CPLDestroyXMLNode( psGeorefData );
419 :
420 0 : return true;
421 : }
422 :
423 : /************************************************************************/
424 : /* Open() */
425 : /************************************************************************/
426 :
427 4464 : GDALDataset *TSXDataset::Open( GDALOpenInfo *poOpenInfo ) {
428 : /* -------------------------------------------------------------------- */
429 : /* Is this a TerraSAR-X product file? */
430 : /* -------------------------------------------------------------------- */
431 4464 : if (!TSXDataset::Identify( poOpenInfo ))
432 : {
433 4464 : return NULL; /* nope */
434 : }
435 :
436 : /* -------------------------------------------------------------------- */
437 : /* Confirm the requested access is supported. */
438 : /* -------------------------------------------------------------------- */
439 0 : if( poOpenInfo->eAccess == GA_Update )
440 : {
441 : CPLError( CE_Failure, CPLE_NotSupported,
442 : "The TSX driver does not support update access to existing"
443 0 : " datasets.\n" );
444 0 : return NULL;
445 : }
446 :
447 0 : CPLString osFilename;
448 :
449 0 : if( poOpenInfo->bIsDirectory )
450 : {
451 : osFilename =
452 0 : CPLFormCIFilename( poOpenInfo->pszFilename, CPLGetFilename( poOpenInfo->pszFilename ), "xml" );
453 : }
454 : else
455 0 : osFilename = poOpenInfo->pszFilename;
456 :
457 : /* Ingest the XML */
458 : CPLXMLNode *psData, *psComponents, *psProductInfo;
459 0 : psData = CPLParseXMLFile( osFilename );
460 0 : if (psData == NULL)
461 0 : return NULL;
462 :
463 : /* find the product components */
464 0 : psComponents = CPLGetXMLNode( psData, "=level1Product.productComponents" );
465 0 : if (psComponents == NULL) {
466 : CPLError( CE_Failure, CPLE_OpenFailed,
467 0 : "Unable to find <productComponents> tag in file.\n" );
468 0 : CPLDestroyXMLNode(psData);
469 0 : return NULL;
470 : }
471 :
472 : /* find the product info tag */
473 0 : psProductInfo = CPLGetXMLNode( psData, "=level1Product.productInfo" );
474 0 : if (psProductInfo == NULL) {
475 : CPLError( CE_Failure, CPLE_OpenFailed,
476 0 : "Unable to find <productInfo> tag in file.\n" );
477 0 : CPLDestroyXMLNode(psData);
478 0 : return NULL;
479 : }
480 :
481 : /* -------------------------------------------------------------------- */
482 : /* Create the dataset. */
483 : /* -------------------------------------------------------------------- */
484 :
485 0 : TSXDataset *poDS = new TSXDataset();
486 :
487 : /* -------------------------------------------------------------------- */
488 : /* Read in product info. */
489 : /* -------------------------------------------------------------------- */
490 :
491 : poDS->SetMetadataItem( "SCENE_CENTRE_TIME", CPLGetXMLValue( psProductInfo,
492 0 : "sceneInfo.sceneCenterCoord.azimuthTimeUTC", "unknown" ) );
493 : poDS->SetMetadataItem( "OPERATIONAL_MODE", CPLGetXMLValue( psProductInfo,
494 0 : "generationInfo.groundOperationsType", "unknown" ) );
495 : poDS->SetMetadataItem( "ORBIT_CYCLE", CPLGetXMLValue( psProductInfo,
496 0 : "missionInfo.orbitCycle", "unknown" ) );
497 : poDS->SetMetadataItem( "ABSOLUTE_ORBIT", CPLGetXMLValue( psProductInfo,
498 0 : "missionInfo.absOrbit", "unknown" ) );
499 : poDS->SetMetadataItem( "ORBIT_DIRECTION", CPLGetXMLValue( psProductInfo,
500 0 : "missionInfo.orbitDirection", "unknown" ) );
501 : poDS->SetMetadataItem( "IMAGING_MODE", CPLGetXMLValue( psProductInfo,
502 0 : "acquisitionInfo.imagingMode", "unknown" ) );
503 : poDS->SetMetadataItem( "PRODUCT_VARIANT", CPLGetXMLValue( psProductInfo,
504 0 : "productVariantInfo.productVariant", "unknown" ) );
505 : char *pszDataType = CPLStrdup( CPLGetXMLValue( psProductInfo,
506 0 : "imageDataInfo.imageDataType", "unknown" ) );
507 0 : poDS->SetMetadataItem( "IMAGE_TYPE", pszDataType );
508 :
509 : /* Get raster information */
510 : int nRows = atoi( CPLGetXMLValue( psProductInfo,
511 0 : "imageDataInfo.imageRaster.numberOfRows", "" ) );
512 : int nCols = atoi( CPLGetXMLValue( psProductInfo,
513 0 : "imageDataInfo.imageRaster.numberOfColumns", "" ) );
514 :
515 0 : poDS->nRasterXSize = nCols;
516 0 : poDS->nRasterYSize = nRows;
517 :
518 : poDS->SetMetadataItem( "ROW_SPACING", CPLGetXMLValue( psProductInfo,
519 0 : "imageDataInfo.imageRaster.rowSpacing", "unknown" ) );
520 : poDS->SetMetadataItem( "COL_SPACING", CPLGetXMLValue( psProductInfo,
521 0 : "imageDataInfo.imageRaster.columnSpacing", "unknown" ) );
522 : poDS->SetMetadataItem( "COL_SPACING_UNITS", CPLGetXMLValue( psProductInfo,
523 0 : "imageDataInfo.imageRaster.columnSpacing.units", "unknown" ) );
524 :
525 : /* Get equivalent number of looks */
526 : poDS->SetMetadataItem( "AZIMUTH_LOOKS", CPLGetXMLValue( psProductInfo,
527 0 : "imageDataInfo.imageRaster.azimuthLooks", "unknown" ) );
528 : poDS->SetMetadataItem( "RANGE_LOOKS", CPLGetXMLValue( psProductInfo,
529 0 : "imageDataInfo.imageRaster.rangeLooks", "unknown" ) );
530 :
531 : const char *pszProductVariant;
532 : pszProductVariant = CPLGetXMLValue( psProductInfo,
533 0 : "productVariantInfo.productVariant", "unknown" );
534 :
535 0 : poDS->SetMetadataItem( "PRODUCT_VARIANT", pszProductVariant );
536 :
537 : /* Determine what product variant this is */
538 0 : if (EQUALN(pszProductVariant,"SSC",3))
539 0 : poDS->nProduct = eSSC;
540 0 : else if (EQUALN(pszProductVariant,"MGD",3))
541 0 : poDS->nProduct = eMGD;
542 0 : else if (EQUALN(pszProductVariant,"EEC",3))
543 0 : poDS->nProduct = eEEC;
544 0 : else if (EQUALN(pszProductVariant,"GEC",3))
545 0 : poDS->nProduct = eGEC;
546 : else
547 0 : poDS->nProduct = eUnknown;
548 :
549 : /* Start reading in the product components */
550 : const char *pszPath;
551 0 : char *pszGeorefFile = NULL;
552 : CPLXMLNode *psComponent;
553 0 : CPLErr geoTransformErr=CE_Failure;
554 0 : for (psComponent = psComponents->psChild; psComponent != NULL;
555 : psComponent = psComponent->psNext)
556 : {
557 0 : const char *pszType = NULL;
558 : pszPath = CPLFormFilename(
559 : CPLGetDirname( osFilename ),
560 : GetFilePath(psComponent, &pszType),
561 0 : "" );
562 0 : const char *pszPolLayer = CPLGetXMLValue(psComponent, "polLayer", " ");
563 :
564 0 : if ( !EQUALN(pszType," ",1) ) {
565 0 : if (EQUALN(pszType, "MAPPING_GRID", 12) ) {
566 : /* the mapping grid... save as a metadata item this path */
567 0 : poDS->SetMetadataItem( "MAPPING_GRID", pszPath );
568 : }
569 0 : else if (EQUALN(pszType, "GEOREF", 6)) {
570 : /* save the path to the georef data for later use */
571 0 : pszGeorefFile = CPLStrdup( pszPath );
572 : }
573 : }
574 0 : else if( !EQUALN(pszPolLayer, " ", 1) &&
575 : EQUALN(psComponent->pszValue, "imageData", 9) ) {
576 : /* determine the polarization of this band */
577 : ePolarization ePol;
578 0 : if ( EQUALN(pszPolLayer, "HH", 2) ) {
579 0 : ePol = HH;
580 : }
581 0 : else if ( EQUALN(pszPolLayer, "HV" , 2) ) {
582 0 : ePol = HV;
583 : }
584 0 : else if ( EQUALN(pszPolLayer, "VH", 2) ) {
585 0 : ePol = VH;
586 : }
587 : else {
588 0 : ePol = VV;
589 : }
590 :
591 : GDALDataType eDataType = EQUALN(pszDataType, "COMPLEX", 7) ?
592 0 : GDT_CInt16 : GDT_UInt16;
593 :
594 : /* try opening the file that represents that band */
595 : TSXRasterBand *poBand;
596 : GDALDataset *poBandData;
597 :
598 0 : poBandData = (GDALDataset *) GDALOpen( pszPath, GA_ReadOnly );
599 0 : if ( poBandData != NULL ) {
600 : poBand = new TSXRasterBand( poDS, eDataType, ePol,
601 0 : poBandData );
602 0 : 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 0 : CPLFree( poDS->pszProjection );
608 0 : poDS->pszProjection = CPLStrdup(poBandData->GetProjectionRef());
609 0 : geoTransformErr = poBandData->GetGeoTransform(poDS->adfGeoTransform);
610 : }
611 : }
612 : }
613 :
614 : //now check if there is a geotransform
615 0 : if ( strcmp(poDS->pszProjection, "") && geoTransformErr==CE_None)
616 : {
617 0 : poDS->bHaveGeoTransform = TRUE;
618 : }
619 : else
620 : {
621 0 : poDS->bHaveGeoTransform = FALSE;
622 0 : CPLFree( poDS->pszProjection );
623 0 : poDS->pszProjection = CPLStrdup("");
624 0 : poDS->adfGeoTransform[0] = 0.0;
625 0 : poDS->adfGeoTransform[1] = 1.0;
626 0 : poDS->adfGeoTransform[2] = 0.0;
627 0 : poDS->adfGeoTransform[3] = 0.0;
628 0 : poDS->adfGeoTransform[4] = 0.0;
629 0 : poDS->adfGeoTransform[5] = 1.0;
630 : }
631 :
632 0 : CPLFree(pszDataType);
633 :
634 :
635 : /* -------------------------------------------------------------------- */
636 : /* Check and set matrix representation. */
637 : /* -------------------------------------------------------------------- */
638 :
639 0 : if (poDS->GetRasterCount() == 4) {
640 0 : poDS->SetMetadataItem( "MATRIX_REPRESENTATION", "SCATTERING" );
641 : }
642 :
643 : /* -------------------------------------------------------------------- */
644 : /* Read the four corners and centre GCPs in */
645 : /* -------------------------------------------------------------------- */
646 :
647 : CPLXMLNode *psSceneInfo = CPLGetXMLNode( psData,
648 0 : "=level1Product.productInfo.sceneInfo" );
649 0 : if (psSceneInfo != NULL)
650 : {
651 : /* extract the GCPs from the provided file */
652 0 : bool success = false;
653 0 : if (pszGeorefFile != NULL)
654 0 : 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 0 : if (!success && poDS->nProduct == eSSC)
659 : {
660 : CPLXMLNode *psNode;
661 0 : int nGCP = 0;
662 : double dfAvgHeight = atof(CPLGetXMLValue(psSceneInfo,
663 0 : "sceneAverageHeight", "0.0"));
664 :
665 : //count and allocate gcps - there should be five - 4 corners and a centre
666 0 : poDS->nGCPCount = 0;
667 0 : for (psNode = psSceneInfo->psChild; psNode != NULL; psNode = psNode->psNext )
668 : {
669 0 : if (!EQUAL(psNode->pszValue, "sceneCenterCoord") &&
670 : !EQUAL(psNode->pszValue, "sceneCornerCoord"))
671 0 : continue;
672 :
673 0 : poDS->nGCPCount++;
674 : }
675 0 : if (poDS->nGCPCount > 0)
676 : {
677 0 : poDS->pasGCPList = (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP), poDS->nGCPCount);
678 :
679 : /* iterate over GCPs */
680 0 : for (psNode = psSceneInfo->psChild; psNode != NULL; psNode = psNode->psNext )
681 : {
682 0 : GDAL_GCP *psGCP = poDS->pasGCPList + nGCP;
683 :
684 0 : if (!EQUAL(psNode->pszValue, "sceneCenterCoord") &&
685 : !EQUAL(psNode->pszValue, "sceneCornerCoord"))
686 0 : continue;
687 :
688 : psGCP->dfGCPPixel = atof(CPLGetXMLValue(psNode, "refColumn",
689 0 : "0.0"));
690 0 : psGCP->dfGCPLine = atof(CPLGetXMLValue(psNode, "refRow", "0.0"));
691 0 : psGCP->dfGCPX = atof(CPLGetXMLValue(psNode, "lon", "0.0"));
692 0 : psGCP->dfGCPY = atof(CPLGetXMLValue(psNode, "lat", "0.0"));
693 0 : psGCP->dfGCPZ = dfAvgHeight;
694 0 : psGCP->pszId = CPLStrdup( CPLSPrintf( "%d", nGCP ) );
695 0 : psGCP->pszInfo = CPLStrdup("");
696 :
697 0 : nGCP++;
698 : }
699 :
700 : //set the projection string - the fields are lat/long - seems to be WGS84 datum
701 0 : OGRSpatialReference osr;
702 0 : osr.SetWellKnownGeogCS( "WGS84" );
703 0 : CPLFree(poDS->pszGCPProjection);
704 0 : osr.exportToWkt( &(poDS->pszGCPProjection) );
705 : }
706 : }
707 :
708 : //gcps override geotransform - does it make sense to have both??
709 0 : if (poDS->nGCPCount>0)
710 : {
711 0 : poDS->bHaveGeoTransform = FALSE;
712 0 : CPLFree( poDS->pszProjection );
713 0 : poDS->pszProjection = CPLStrdup("");
714 0 : poDS->adfGeoTransform[0] = 0.0;
715 0 : poDS->adfGeoTransform[1] = 1.0;
716 0 : poDS->adfGeoTransform[2] = 0.0;
717 0 : poDS->adfGeoTransform[3] = 0.0;
718 0 : poDS->adfGeoTransform[4] = 0.0;
719 0 : poDS->adfGeoTransform[5] = 1.0;
720 : }
721 :
722 : }
723 : else {
724 : CPLError(CE_Warning, CPLE_AppDefined,
725 : "Unable to find sceneInfo tag in XML document. "
726 0 : "Proceeding with caution.");
727 : }
728 :
729 0 : CPLFree(pszGeorefFile);
730 :
731 : /* -------------------------------------------------------------------- */
732 : /* Initialize any PAM information. */
733 : /* -------------------------------------------------------------------- */
734 0 : poDS->SetDescription( poOpenInfo->pszFilename );
735 0 : poDS->TryLoadXML();
736 :
737 : /* -------------------------------------------------------------------- */
738 : /* Check for overviews. */
739 : /* -------------------------------------------------------------------- */
740 0 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
741 :
742 0 : CPLDestroyXMLNode(psData);
743 :
744 0 : return poDS;
745 : }
746 :
747 : /************************************************************************/
748 : /* GetGCPCount() */
749 : /************************************************************************/
750 :
751 0 : int TSXDataset::GetGCPCount() {
752 0 : return nGCPCount;
753 : }
754 :
755 : /************************************************************************/
756 : /* GetGCPProjection() */
757 : /************************************************************************/
758 :
759 0 : const char *TSXDataset::GetGCPProjection() {
760 0 : return pszGCPProjection;
761 : }
762 :
763 : /************************************************************************/
764 : /* GetGCPs() */
765 : /************************************************************************/
766 :
767 0 : const GDAL_GCP *TSXDataset::GetGCPs() {
768 0 : return pasGCPList;
769 : }
770 :
771 : /************************************************************************/
772 : /* GetProjectionRef() */
773 : /************************************************************************/
774 0 : const char *TSXDataset::GetProjectionRef()
775 : {
776 0 : return pszProjection;
777 : }
778 :
779 : /************************************************************************/
780 : /* GetGeotransform() */
781 : /************************************************************************/
782 0 : CPLErr TSXDataset::GetGeoTransform(double* padfTransform)
783 : {
784 0 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
785 :
786 0 : if (bHaveGeoTransform)
787 0 : return( CE_None );
788 :
789 0 : return( CE_Failure );
790 : }
791 :
792 : /************************************************************************/
793 : /* GDALRegister_TSX() */
794 : /************************************************************************/
795 :
796 1135 : void GDALRegister_TSX() {
797 : GDALDriver *poDriver;
798 :
799 1135 : if( GDALGetDriverByName( "TSX" ) == NULL )
800 : {
801 1093 : poDriver = new GDALDriver();
802 :
803 1093 : poDriver->SetDescription( "TSX" );
804 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
805 1093 : "TerraSAR-X Product" );
806 : /* poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "frmt_tsx.html" ); */
807 :
808 1093 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
809 :
810 1093 : poDriver->pfnOpen = TSXDataset::Open;
811 1093 : poDriver->pfnIdentify = TSXDataset::Identify;
812 :
813 1093 : GetGDALDriverManager()->RegisterDriver( poDriver );
814 : }
815 1135 : }
816 :
|
