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