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 : 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(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(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 : }
223 :
224 : /************************************************************************/
225 : /* ~TSXDataset() */
226 : /************************************************************************/
227 :
228 0 : TSXDataset::~TSXDataset()
229 : {
230 0 : FlushCache(true);
231 :
232 0 : if (nGCPCount > 0)
233 : {
234 0 : GDALDeinitGCPs(nGCPCount, pasGCPList);
235 0 : CPLFree(pasGCPList);
236 : }
237 0 : }
238 :
239 : /************************************************************************/
240 : /* Identify() */
241 : /************************************************************************/
242 :
243 58772 : int TSXDataset::Identify(GDALOpenInfo *poOpenInfo)
244 : {
245 58772 : if (poOpenInfo->fpL == nullptr || poOpenInfo->nHeaderBytes < 260)
246 : {
247 55059 : if (poOpenInfo->bIsDirectory)
248 : {
249 0 : const CPLString osFilename = CPLFormCIFilenameSafe(
250 612 : poOpenInfo->pszFilename,
251 1224 : CPLGetFilename(poOpenInfo->pszFilename), "xml");
252 :
253 : /* Check if the filename contains TSX1_SAR (TerraSAR-X) or TDX1_SAR
254 : * (TanDEM-X) or PAZ1_SAR (PAZ) */
255 1836 : if (!(STARTS_WITH_CI(CPLGetBasenameSafe(osFilename).c_str(),
256 1224 : "TSX1_SAR") ||
257 1224 : STARTS_WITH_CI(CPLGetBasenameSafe(osFilename).c_str(),
258 : "TDX1_SAR") ||
259 1224 : STARTS_WITH_CI(CPLGetBasenameSafe(osFilename).c_str(),
260 : "PAZ1_SAR")))
261 612 : return 0;
262 :
263 : VSIStatBufL sStat;
264 0 : if (VSIStatL(osFilename, &sStat) == 0)
265 0 : return 1;
266 : }
267 :
268 54447 : return 0;
269 : }
270 :
271 : /* Check if the filename contains TSX1_SAR (TerraSAR-X) or TDX1_SAR
272 : * (TanDEM-X) or PAZ1_SAR (PAZ) */
273 11139 : if (!(STARTS_WITH_CI(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str(),
274 7426 : "TSX1_SAR") ||
275 7426 : STARTS_WITH_CI(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str(),
276 : "TDX1_SAR") ||
277 7426 : STARTS_WITH_CI(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str(),
278 : "PAZ1_SAR")))
279 3713 : 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 : ReportUpdateNotSupportedByDriver("TSX");
453 0 : return nullptr;
454 : }
455 :
456 0 : CPLString osFilename;
457 :
458 0 : if (poOpenInfo->bIsDirectory)
459 : {
460 0 : osFilename = CPLFormCIFilenameSafe(
461 0 : poOpenInfo->pszFilename, CPLGetFilename(poOpenInfo->pszFilename),
462 0 : "xml");
463 : }
464 : else
465 0 : osFilename = poOpenInfo->pszFilename;
466 :
467 : /* Ingest the XML */
468 0 : CPLXMLNode *psData = CPLParseXMLFile(osFilename);
469 0 : if (psData == nullptr)
470 0 : return nullptr;
471 :
472 : /* find the product components */
473 : CPLXMLNode *psComponents =
474 0 : CPLGetXMLNode(psData, "=level1Product.productComponents");
475 0 : if (psComponents == nullptr)
476 : {
477 0 : CPLError(CE_Failure, CPLE_OpenFailed,
478 : "Unable to find <productComponents> tag in file.\n");
479 0 : CPLDestroyXMLNode(psData);
480 0 : return nullptr;
481 : }
482 :
483 : /* find the product info tag */
484 : CPLXMLNode *psProductInfo =
485 0 : CPLGetXMLNode(psData, "=level1Product.productInfo");
486 0 : if (psProductInfo == nullptr)
487 : {
488 0 : CPLError(CE_Failure, CPLE_OpenFailed,
489 : "Unable to find <productInfo> tag in file.\n");
490 0 : CPLDestroyXMLNode(psData);
491 0 : return nullptr;
492 : }
493 :
494 : /* -------------------------------------------------------------------- */
495 : /* Create the dataset. */
496 : /* -------------------------------------------------------------------- */
497 :
498 0 : TSXDataset *poDS = new TSXDataset();
499 :
500 : /* -------------------------------------------------------------------- */
501 : /* Read in product info. */
502 : /* -------------------------------------------------------------------- */
503 :
504 0 : poDS->SetMetadataItem(
505 : "SCENE_CENTRE_TIME",
506 : CPLGetXMLValue(psProductInfo,
507 0 : "sceneInfo.sceneCenterCoord.azimuthTimeUTC", "unknown"));
508 0 : poDS->SetMetadataItem("OPERATIONAL_MODE",
509 : CPLGetXMLValue(psProductInfo,
510 : "generationInfo.groundOperationsType",
511 0 : "unknown"));
512 0 : poDS->SetMetadataItem(
513 : "ORBIT_CYCLE",
514 0 : CPLGetXMLValue(psProductInfo, "missionInfo.orbitCycle", "unknown"));
515 0 : poDS->SetMetadataItem(
516 : "ABSOLUTE_ORBIT",
517 0 : CPLGetXMLValue(psProductInfo, "missionInfo.absOrbit", "unknown"));
518 0 : poDS->SetMetadataItem(
519 : "ORBIT_DIRECTION",
520 0 : CPLGetXMLValue(psProductInfo, "missionInfo.orbitDirection", "unknown"));
521 0 : poDS->SetMetadataItem("IMAGING_MODE",
522 : CPLGetXMLValue(psProductInfo,
523 : "acquisitionInfo.imagingMode",
524 0 : "unknown"));
525 0 : poDS->SetMetadataItem("PRODUCT_VARIANT",
526 : CPLGetXMLValue(psProductInfo,
527 : "productVariantInfo.productVariant",
528 0 : "unknown"));
529 0 : char *pszDataType = CPLStrdup(CPLGetXMLValue(
530 : psProductInfo, "imageDataInfo.imageDataType", "unknown"));
531 0 : poDS->SetMetadataItem("IMAGE_TYPE", pszDataType);
532 :
533 : /* Get raster information */
534 0 : int nRows = atoi(CPLGetXMLValue(
535 : psProductInfo, "imageDataInfo.imageRaster.numberOfRows", ""));
536 0 : int nCols = atoi(CPLGetXMLValue(
537 : psProductInfo, "imageDataInfo.imageRaster.numberOfColumns", ""));
538 :
539 0 : poDS->nRasterXSize = nCols;
540 0 : poDS->nRasterYSize = nRows;
541 :
542 0 : poDS->SetMetadataItem("ROW_SPACING",
543 : CPLGetXMLValue(psProductInfo,
544 : "imageDataInfo.imageRaster.rowSpacing",
545 0 : "unknown"));
546 0 : poDS->SetMetadataItem(
547 : "COL_SPACING",
548 : CPLGetXMLValue(psProductInfo, "imageDataInfo.imageRaster.columnSpacing",
549 0 : "unknown"));
550 0 : poDS->SetMetadataItem(
551 : "COL_SPACING_UNITS",
552 : CPLGetXMLValue(psProductInfo,
553 : "imageDataInfo.imageRaster.columnSpacing.units",
554 0 : "unknown"));
555 :
556 : /* Get equivalent number of looks */
557 0 : poDS->SetMetadataItem(
558 : "AZIMUTH_LOOKS",
559 : CPLGetXMLValue(psProductInfo, "imageDataInfo.imageRaster.azimuthLooks",
560 0 : "unknown"));
561 0 : poDS->SetMetadataItem("RANGE_LOOKS",
562 : CPLGetXMLValue(psProductInfo,
563 : "imageDataInfo.imageRaster.rangeLooks",
564 0 : "unknown"));
565 :
566 0 : const char *pszProductVariant = CPLGetXMLValue(
567 : psProductInfo, "productVariantInfo.productVariant", "unknown");
568 :
569 0 : poDS->SetMetadataItem("PRODUCT_VARIANT", pszProductVariant);
570 :
571 : /* Determine what product variant this is */
572 0 : if (STARTS_WITH_CI(pszProductVariant, "SSC"))
573 0 : poDS->nProduct = eSSC;
574 0 : else if (STARTS_WITH_CI(pszProductVariant, "MGD"))
575 0 : poDS->nProduct = eMGD;
576 0 : else if (STARTS_WITH_CI(pszProductVariant, "EEC"))
577 0 : poDS->nProduct = eEEC;
578 0 : else if (STARTS_WITH_CI(pszProductVariant, "GEC"))
579 0 : poDS->nProduct = eGEC;
580 : else
581 0 : poDS->nProduct = eUnknown;
582 :
583 : /* Start reading in the product components */
584 0 : char *pszGeorefFile = nullptr;
585 0 : CPLErr geoTransformErr = CE_Failure;
586 0 : for (CPLXMLNode *psComponent = psComponents->psChild;
587 0 : psComponent != nullptr; psComponent = psComponent->psNext)
588 : {
589 0 : const char *pszType = nullptr;
590 : const std::string osPath =
591 0 : CPLFormFilenameSafe(CPLGetDirnameSafe(osFilename).c_str(),
592 0 : GetFilePath(psComponent, &pszType).c_str(), "");
593 0 : const char *pszPolLayer = CPLGetXMLValue(psComponent, "polLayer", " ");
594 :
595 0 : if (!STARTS_WITH_CI(pszType, " "))
596 : {
597 0 : if (STARTS_WITH_CI(pszType, "MAPPING_GRID"))
598 : {
599 : /* the mapping grid... save as a metadata item this path */
600 0 : poDS->SetMetadataItem("MAPPING_GRID", osPath.c_str());
601 : }
602 0 : else if (STARTS_WITH_CI(pszType, "GEOREF"))
603 : {
604 : /* save the path to the georef data for later use */
605 0 : CPLFree(pszGeorefFile);
606 0 : pszGeorefFile = CPLStrdup(osPath.c_str());
607 : }
608 : }
609 0 : else if (!STARTS_WITH_CI(pszPolLayer, " ") &&
610 0 : STARTS_WITH_CI(psComponent->pszValue, "imageData"))
611 : {
612 : /* determine the polarization of this band */
613 : ePolarization ePol;
614 0 : if (STARTS_WITH_CI(pszPolLayer, "HH"))
615 : {
616 0 : ePol = HH;
617 : }
618 0 : else if (STARTS_WITH_CI(pszPolLayer, "HV"))
619 : {
620 0 : ePol = HV;
621 : }
622 0 : else if (STARTS_WITH_CI(pszPolLayer, "VH"))
623 : {
624 0 : ePol = VH;
625 : }
626 : else
627 : {
628 0 : ePol = VV;
629 : }
630 :
631 0 : GDALDataType eDataType = STARTS_WITH_CI(pszDataType, "COMPLEX")
632 0 : ? GDT_CInt16
633 : : GDT_UInt16;
634 :
635 : /* try opening the file that represents that band */
636 : GDALDataset *poBandData =
637 0 : GDALDataset::FromHandle(GDALOpen(osPath.c_str(), GA_ReadOnly));
638 0 : if (poBandData != nullptr)
639 : {
640 : TSXRasterBand *poBand =
641 0 : new TSXRasterBand(poDS, eDataType, ePol, poBandData);
642 0 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand);
643 :
644 : // copy georeferencing info from the band
645 : // need error checking??
646 : // it will just save the info from the last band
647 0 : const auto poSrcSRS = poBandData->GetSpatialRef();
648 0 : if (poSrcSRS)
649 0 : poDS->m_oSRS = *poSrcSRS;
650 :
651 0 : geoTransformErr = poBandData->GetGeoTransform(poDS->m_gt);
652 : }
653 : }
654 : }
655 :
656 : // now check if there is a geotransform
657 0 : if (!poDS->m_oSRS.IsEmpty() && geoTransformErr == CE_None)
658 : {
659 0 : poDS->bHaveGeoTransform = TRUE;
660 : }
661 : else
662 : {
663 0 : poDS->bHaveGeoTransform = FALSE;
664 0 : poDS->m_oSRS.Clear();
665 0 : poDS->m_gt = GDALGeoTransform();
666 : }
667 :
668 0 : CPLFree(pszDataType);
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 : CPLXMLNode *psSceneInfo =
684 0 : CPLGetXMLNode(psData, "=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 (pszGeorefFile != nullptr)
690 0 : success = poDS->getGCPsFromGEOREF_XML(pszGeorefFile);
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 : 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 0 : CPLFree(pszGeorefFile);
765 :
766 : /* -------------------------------------------------------------------- */
767 : /* Initialize any PAM information. */
768 : /* -------------------------------------------------------------------- */
769 0 : poDS->SetDescription(poOpenInfo->pszFilename);
770 0 : poDS->TryLoadXML();
771 :
772 : /* -------------------------------------------------------------------- */
773 : /* Check for overviews. */
774 : /* -------------------------------------------------------------------- */
775 0 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
776 :
777 0 : CPLDestroyXMLNode(psData);
778 :
779 0 : return poDS;
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 1935 : void GDALRegister_TSX()
837 : {
838 1935 : if (GDALGetDriverByName("TSX") != nullptr)
839 282 : return;
840 :
841 1653 : GDALDriver *poDriver = new GDALDriver();
842 :
843 1653 : poDriver->SetDescription("TSX");
844 1653 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
845 1653 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "TerraSAR-X Product");
846 1653 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/tsx.html");
847 1653 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
848 :
849 1653 : poDriver->pfnOpen = TSXDataset::Open;
850 1653 : poDriver->pfnIdentify = TSXDataset::Identify;
851 :
852 1653 : GetGDALDriverManager()->RegisterDriver(poDriver);
853 : }
|
