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