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