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