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