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