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 53597 : int TSXDataset::Identify(GDALOpenInfo *poOpenInfo)
250 : {
251 53597 : if (poOpenInfo->fpL == nullptr || poOpenInfo->nHeaderBytes < 260)
252 : {
253 49919 : if (poOpenInfo->bIsDirectory)
254 : {
255 0 : const CPLString osFilename = CPLFormCIFilenameSafe(
256 423 : poOpenInfo->pszFilename,
257 846 : 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 1269 : if (!(STARTS_WITH_CI(CPLGetBasenameSafe(osFilename).c_str(),
262 846 : "TSX1_SAR") ||
263 846 : STARTS_WITH_CI(CPLGetBasenameSafe(osFilename).c_str(),
264 : "TDX1_SAR") ||
265 846 : STARTS_WITH_CI(CPLGetBasenameSafe(osFilename).c_str(),
266 : "PAZ1_SAR")))
267 423 : return 0;
268 :
269 : VSIStatBufL sStat;
270 0 : if (VSIStatL(osFilename, &sStat) == 0)
271 0 : return 1;
272 : }
273 :
274 49496 : return 0;
275 : }
276 :
277 : /* Check if the filename contains TSX1_SAR (TerraSAR-X) or TDX1_SAR
278 : * (TanDEM-X) or PAZ1_SAR (PAZ) */
279 11034 : if (!(STARTS_WITH_CI(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str(),
280 7356 : "TSX1_SAR") ||
281 7356 : STARTS_WITH_CI(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str(),
282 : "TDX1_SAR") ||
283 7356 : STARTS_WITH_CI(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str(),
284 : "PAZ1_SAR")))
285 3678 : return 0;
286 :
287 : /* finally look for the <level1Product tag */
288 0 : if (!STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
289 : "<level1Product"))
290 0 : return 0;
291 :
292 0 : return 1;
293 : }
294 :
295 : /************************************************************************/
296 : /* getGCPsFromGEOREF_XML() */
297 : /*Reads georeferencing information from the TerraSAR-X GEOREF.XML file */
298 : /*and writes the information to the dataset's gcp list and projection */
299 : /*string. */
300 : /*Returns true on success. */
301 : /************************************************************************/
302 0 : bool TSXDataset::getGCPsFromGEOREF_XML(char *pszGeorefFilename)
303 : {
304 : // open GEOREF.xml
305 0 : CPLXMLNode *psGeorefData = CPLParseXMLFile(pszGeorefFilename);
306 0 : if (psGeorefData == nullptr)
307 0 : return false;
308 :
309 : // get the ellipsoid and semi-major, semi-minor axes
310 0 : OGRSpatialReference osr;
311 : CPLXMLNode *psSphere =
312 0 : CPLGetXMLNode(psGeorefData, "=geoReference.referenceFrames.sphere");
313 0 : if (psSphere != nullptr)
314 : {
315 : const char *pszEllipsoidName =
316 0 : CPLGetXMLValue(psSphere, "ellipsoidID", "");
317 : const double minor_axis =
318 0 : CPLAtof(CPLGetXMLValue(psSphere, "semiMinorAxis", "0.0"));
319 : const double major_axis =
320 0 : CPLAtof(CPLGetXMLValue(psSphere, "semiMajorAxis", "0.0"));
321 : // save datum parameters to the spatial reference
322 0 : if (EQUAL(pszEllipsoidName, "") || minor_axis == 0.0 ||
323 : major_axis == 0.0)
324 : {
325 0 : CPLError(CE_Warning, CPLE_AppDefined,
326 : "Warning- incomplete"
327 : " ellipsoid information. Using wgs-84 parameters.\n");
328 0 : osr.SetWellKnownGeogCS("WGS84");
329 : }
330 0 : else if (EQUAL(pszEllipsoidName, "WGS84"))
331 : {
332 0 : osr.SetWellKnownGeogCS("WGS84");
333 : }
334 : else
335 : {
336 0 : const double inv_flattening =
337 0 : major_axis / (major_axis - minor_axis);
338 0 : osr.SetGeogCS("", "", pszEllipsoidName, major_axis, inv_flattening);
339 : }
340 : }
341 :
342 : // get gcps
343 : CPLXMLNode *psGeolocationGrid =
344 0 : CPLGetXMLNode(psGeorefData, "=geoReference.geolocationGrid");
345 0 : if (psGeolocationGrid == nullptr)
346 : {
347 0 : CPLDestroyXMLNode(psGeorefData);
348 0 : return false;
349 : }
350 0 : nGCPCount = atoi(
351 : CPLGetXMLValue(psGeolocationGrid, "numberOfGridPoints.total", "0"));
352 : // count the gcps if the given count value is invalid
353 0 : CPLXMLNode *psNode = nullptr;
354 0 : if (nGCPCount <= 0)
355 : {
356 0 : for (psNode = psGeolocationGrid->psChild; psNode != nullptr;
357 0 : psNode = psNode->psNext)
358 0 : if (EQUAL(psNode->pszValue, "gridPoint"))
359 0 : nGCPCount++;
360 : }
361 : // if there are no gcps, fail
362 0 : if (nGCPCount <= 0)
363 : {
364 0 : CPLDestroyXMLNode(psGeorefData);
365 0 : return false;
366 : }
367 :
368 : // put some reasonable limits of the number of gcps
369 0 : if (nGCPCount > MAX_GCPS)
370 0 : nGCPCount = MAX_GCPS;
371 : // allocate memory for the gcps
372 0 : pasGCPList =
373 0 : reinterpret_cast<GDAL_GCP *>(CPLCalloc(sizeof(GDAL_GCP), nGCPCount));
374 :
375 : // loop through all gcps and set info
376 :
377 : // save the number allocated to ensure it does not run off the end of the
378 : // array
379 0 : const int gcps_allocated = nGCPCount;
380 0 : nGCPCount = 0; // reset to zero and count
381 : // do a check on the grid point to make sure it has lat,long row, and column
382 : // it seems that only SSC products contain row, col - how to map lat long
383 : // otherwise?? for now fail if row and col are not present - just check the
384 : // first and assume the rest are the same
385 0 : for (psNode = psGeolocationGrid->psChild; psNode != nullptr;
386 0 : psNode = psNode->psNext)
387 : {
388 0 : if (!EQUAL(psNode->pszValue, "gridPoint"))
389 0 : continue;
390 :
391 0 : if (!strcmp(CPLGetXMLValue(psNode, "col", "error"), "error") ||
392 0 : !strcmp(CPLGetXMLValue(psNode, "row", "error"), "error") ||
393 0 : !strcmp(CPLGetXMLValue(psNode, "lon", "error"), "error") ||
394 0 : !strcmp(CPLGetXMLValue(psNode, "lat", "error"), "error"))
395 : {
396 0 : CPLDestroyXMLNode(psGeorefData);
397 0 : return false;
398 : }
399 : }
400 0 : for (psNode = psGeolocationGrid->psChild; psNode != nullptr;
401 0 : psNode = psNode->psNext)
402 : {
403 : // break out if the end of the array has been reached
404 0 : if (nGCPCount >= gcps_allocated)
405 : {
406 0 : CPLError(CE_Warning, CPLE_AppDefined,
407 : "GDAL TSX driver: Truncating the number of GCPs.");
408 0 : break;
409 : }
410 :
411 0 : GDAL_GCP *psGCP = pasGCPList + nGCPCount;
412 :
413 0 : if (!EQUAL(psNode->pszValue, "gridPoint"))
414 0 : continue;
415 :
416 0 : nGCPCount++;
417 :
418 : char szID[32];
419 0 : snprintf(szID, sizeof(szID), "%d", nGCPCount);
420 0 : psGCP->pszId = CPLStrdup(szID);
421 0 : psGCP->pszInfo = CPLStrdup("");
422 0 : psGCP->dfGCPPixel = CPLAtof(CPLGetXMLValue(psNode, "col", "0"));
423 0 : psGCP->dfGCPLine = CPLAtof(CPLGetXMLValue(psNode, "row", "0"));
424 0 : psGCP->dfGCPX = CPLAtof(CPLGetXMLValue(psNode, "lon", ""));
425 0 : psGCP->dfGCPY = CPLAtof(CPLGetXMLValue(psNode, "lat", ""));
426 : // looks like height is in meters - should it be converted so xyz are
427 : // all on the same scale??
428 0 : psGCP->dfGCPZ = 0;
429 : // CPLAtof(CPLGetXMLValue(psNode,"height",""));
430 : }
431 :
432 0 : m_oGCPSRS = std::move(osr);
433 :
434 0 : CPLDestroyXMLNode(psGeorefData);
435 :
436 0 : return true;
437 : }
438 :
439 : /************************************************************************/
440 : /* Open() */
441 : /************************************************************************/
442 :
443 0 : GDALDataset *TSXDataset::Open(GDALOpenInfo *poOpenInfo)
444 : {
445 : /* -------------------------------------------------------------------- */
446 : /* Is this a TerraSAR-X product file? */
447 : /* -------------------------------------------------------------------- */
448 0 : if (!TSXDataset::Identify(poOpenInfo))
449 : {
450 0 : return nullptr; /* nope */
451 : }
452 :
453 : /* -------------------------------------------------------------------- */
454 : /* Confirm the requested access is supported. */
455 : /* -------------------------------------------------------------------- */
456 0 : if (poOpenInfo->eAccess == GA_Update)
457 : {
458 0 : ReportUpdateNotSupportedByDriver("TSX");
459 0 : return nullptr;
460 : }
461 :
462 0 : CPLString osFilename;
463 :
464 0 : if (poOpenInfo->bIsDirectory)
465 : {
466 0 : osFilename = CPLFormCIFilenameSafe(
467 0 : poOpenInfo->pszFilename, CPLGetFilename(poOpenInfo->pszFilename),
468 0 : "xml");
469 : }
470 : else
471 0 : osFilename = poOpenInfo->pszFilename;
472 :
473 : /* Ingest the XML */
474 0 : CPLXMLNode *psData = CPLParseXMLFile(osFilename);
475 0 : if (psData == nullptr)
476 0 : return nullptr;
477 :
478 : /* find the product components */
479 : CPLXMLNode *psComponents =
480 0 : CPLGetXMLNode(psData, "=level1Product.productComponents");
481 0 : if (psComponents == nullptr)
482 : {
483 0 : CPLError(CE_Failure, CPLE_OpenFailed,
484 : "Unable to find <productComponents> tag in file.\n");
485 0 : CPLDestroyXMLNode(psData);
486 0 : return nullptr;
487 : }
488 :
489 : /* find the product info tag */
490 : CPLXMLNode *psProductInfo =
491 0 : CPLGetXMLNode(psData, "=level1Product.productInfo");
492 0 : if (psProductInfo == nullptr)
493 : {
494 0 : CPLError(CE_Failure, CPLE_OpenFailed,
495 : "Unable to find <productInfo> tag in file.\n");
496 0 : CPLDestroyXMLNode(psData);
497 0 : return nullptr;
498 : }
499 :
500 : /* -------------------------------------------------------------------- */
501 : /* Create the dataset. */
502 : /* -------------------------------------------------------------------- */
503 :
504 0 : TSXDataset *poDS = new TSXDataset();
505 :
506 : /* -------------------------------------------------------------------- */
507 : /* Read in product info. */
508 : /* -------------------------------------------------------------------- */
509 :
510 0 : poDS->SetMetadataItem(
511 : "SCENE_CENTRE_TIME",
512 : CPLGetXMLValue(psProductInfo,
513 0 : "sceneInfo.sceneCenterCoord.azimuthTimeUTC", "unknown"));
514 0 : poDS->SetMetadataItem("OPERATIONAL_MODE",
515 : CPLGetXMLValue(psProductInfo,
516 : "generationInfo.groundOperationsType",
517 0 : "unknown"));
518 0 : poDS->SetMetadataItem(
519 : "ORBIT_CYCLE",
520 0 : CPLGetXMLValue(psProductInfo, "missionInfo.orbitCycle", "unknown"));
521 0 : poDS->SetMetadataItem(
522 : "ABSOLUTE_ORBIT",
523 0 : CPLGetXMLValue(psProductInfo, "missionInfo.absOrbit", "unknown"));
524 0 : poDS->SetMetadataItem(
525 : "ORBIT_DIRECTION",
526 0 : CPLGetXMLValue(psProductInfo, "missionInfo.orbitDirection", "unknown"));
527 0 : poDS->SetMetadataItem("IMAGING_MODE",
528 : CPLGetXMLValue(psProductInfo,
529 : "acquisitionInfo.imagingMode",
530 0 : "unknown"));
531 0 : poDS->SetMetadataItem("PRODUCT_VARIANT",
532 : CPLGetXMLValue(psProductInfo,
533 : "productVariantInfo.productVariant",
534 0 : "unknown"));
535 0 : char *pszDataType = CPLStrdup(CPLGetXMLValue(
536 : psProductInfo, "imageDataInfo.imageDataType", "unknown"));
537 0 : poDS->SetMetadataItem("IMAGE_TYPE", pszDataType);
538 :
539 : /* Get raster information */
540 0 : int nRows = atoi(CPLGetXMLValue(
541 : psProductInfo, "imageDataInfo.imageRaster.numberOfRows", ""));
542 0 : int nCols = atoi(CPLGetXMLValue(
543 : psProductInfo, "imageDataInfo.imageRaster.numberOfColumns", ""));
544 :
545 0 : poDS->nRasterXSize = nCols;
546 0 : poDS->nRasterYSize = nRows;
547 :
548 0 : poDS->SetMetadataItem("ROW_SPACING",
549 : CPLGetXMLValue(psProductInfo,
550 : "imageDataInfo.imageRaster.rowSpacing",
551 0 : "unknown"));
552 0 : poDS->SetMetadataItem(
553 : "COL_SPACING",
554 : CPLGetXMLValue(psProductInfo, "imageDataInfo.imageRaster.columnSpacing",
555 0 : "unknown"));
556 0 : poDS->SetMetadataItem(
557 : "COL_SPACING_UNITS",
558 : CPLGetXMLValue(psProductInfo,
559 : "imageDataInfo.imageRaster.columnSpacing.units",
560 0 : "unknown"));
561 :
562 : /* Get equivalent number of looks */
563 0 : poDS->SetMetadataItem(
564 : "AZIMUTH_LOOKS",
565 : CPLGetXMLValue(psProductInfo, "imageDataInfo.imageRaster.azimuthLooks",
566 0 : "unknown"));
567 0 : poDS->SetMetadataItem("RANGE_LOOKS",
568 : CPLGetXMLValue(psProductInfo,
569 : "imageDataInfo.imageRaster.rangeLooks",
570 0 : "unknown"));
571 :
572 0 : const char *pszProductVariant = CPLGetXMLValue(
573 : psProductInfo, "productVariantInfo.productVariant", "unknown");
574 :
575 0 : poDS->SetMetadataItem("PRODUCT_VARIANT", pszProductVariant);
576 :
577 : /* Determine what product variant this is */
578 0 : if (STARTS_WITH_CI(pszProductVariant, "SSC"))
579 0 : poDS->nProduct = eSSC;
580 0 : else if (STARTS_WITH_CI(pszProductVariant, "MGD"))
581 0 : poDS->nProduct = eMGD;
582 0 : else if (STARTS_WITH_CI(pszProductVariant, "EEC"))
583 0 : poDS->nProduct = eEEC;
584 0 : else if (STARTS_WITH_CI(pszProductVariant, "GEC"))
585 0 : poDS->nProduct = eGEC;
586 : else
587 0 : poDS->nProduct = eUnknown;
588 :
589 : /* Start reading in the product components */
590 0 : char *pszGeorefFile = nullptr;
591 0 : CPLErr geoTransformErr = CE_Failure;
592 0 : for (CPLXMLNode *psComponent = psComponents->psChild;
593 0 : psComponent != nullptr; psComponent = psComponent->psNext)
594 : {
595 0 : const char *pszType = nullptr;
596 : const std::string osPath =
597 0 : CPLFormFilenameSafe(CPLGetDirnameSafe(osFilename).c_str(),
598 0 : GetFilePath(psComponent, &pszType).c_str(), "");
599 0 : const char *pszPolLayer = CPLGetXMLValue(psComponent, "polLayer", " ");
600 :
601 0 : if (!STARTS_WITH_CI(pszType, " "))
602 : {
603 0 : if (STARTS_WITH_CI(pszType, "MAPPING_GRID"))
604 : {
605 : /* the mapping grid... save as a metadata item this path */
606 0 : poDS->SetMetadataItem("MAPPING_GRID", osPath.c_str());
607 : }
608 0 : else if (STARTS_WITH_CI(pszType, "GEOREF"))
609 : {
610 : /* save the path to the georef data for later use */
611 0 : CPLFree(pszGeorefFile);
612 0 : pszGeorefFile = CPLStrdup(osPath.c_str());
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 0 : GDALDataType eDataType = STARTS_WITH_CI(pszDataType, "COMPLEX")
638 0 : ? GDT_CInt16
639 : : 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, 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 : geoTransformErr =
658 0 : poBandData->GetGeoTransform(poDS->adfGeoTransform);
659 : }
660 : }
661 : }
662 :
663 : // now check if there is a geotransform
664 0 : if (!poDS->m_oSRS.IsEmpty() && geoTransformErr == CE_None)
665 : {
666 0 : poDS->bHaveGeoTransform = TRUE;
667 : }
668 : else
669 : {
670 0 : poDS->bHaveGeoTransform = FALSE;
671 0 : poDS->m_oSRS.Clear();
672 0 : poDS->adfGeoTransform[0] = 0.0;
673 0 : poDS->adfGeoTransform[1] = 1.0;
674 0 : poDS->adfGeoTransform[2] = 0.0;
675 0 : poDS->adfGeoTransform[3] = 0.0;
676 0 : poDS->adfGeoTransform[4] = 0.0;
677 0 : poDS->adfGeoTransform[5] = 1.0;
678 : }
679 :
680 0 : CPLFree(pszDataType);
681 :
682 : /* -------------------------------------------------------------------- */
683 : /* Check and set matrix representation. */
684 : /* -------------------------------------------------------------------- */
685 :
686 0 : if (poDS->GetRasterCount() == 4)
687 : {
688 0 : poDS->SetMetadataItem("MATRIX_REPRESENTATION", "SCATTERING");
689 : }
690 :
691 : /* -------------------------------------------------------------------- */
692 : /* Read the four corners and centre GCPs in */
693 : /* -------------------------------------------------------------------- */
694 :
695 : CPLXMLNode *psSceneInfo =
696 0 : CPLGetXMLNode(psData, "=level1Product.productInfo.sceneInfo");
697 0 : if (psSceneInfo != nullptr)
698 : {
699 : /* extract the GCPs from the provided file */
700 0 : bool success = false;
701 0 : if (pszGeorefFile != nullptr)
702 0 : success = poDS->getGCPsFromGEOREF_XML(pszGeorefFile);
703 :
704 : // if the gcp's cannot be extracted from the georef file, try to get the
705 : // corner coordinates for now just SSC because the others don't have
706 : // refColumn and refRow
707 0 : if (!success && poDS->nProduct == eSSC)
708 : {
709 0 : int nGCP = 0;
710 0 : double dfAvgHeight = CPLAtof(
711 : CPLGetXMLValue(psSceneInfo, "sceneAverageHeight", "0.0"));
712 :
713 : // count and allocate gcps - there should be five - 4 corners and a
714 : // centre
715 0 : poDS->nGCPCount = 0;
716 0 : CPLXMLNode *psNode = psSceneInfo->psChild;
717 0 : for (; psNode != nullptr; psNode = psNode->psNext)
718 : {
719 0 : if (!EQUAL(psNode->pszValue, "sceneCenterCoord") &&
720 0 : !EQUAL(psNode->pszValue, "sceneCornerCoord"))
721 0 : continue;
722 :
723 0 : poDS->nGCPCount++;
724 : }
725 0 : if (poDS->nGCPCount > 0)
726 : {
727 0 : poDS->pasGCPList =
728 0 : (GDAL_GCP *)CPLCalloc(sizeof(GDAL_GCP), poDS->nGCPCount);
729 :
730 : /* iterate over GCPs */
731 0 : for (psNode = psSceneInfo->psChild; psNode != nullptr;
732 0 : psNode = psNode->psNext)
733 : {
734 0 : GDAL_GCP *psGCP = poDS->pasGCPList + nGCP;
735 :
736 0 : if (!EQUAL(psNode->pszValue, "sceneCenterCoord") &&
737 0 : !EQUAL(psNode->pszValue, "sceneCornerCoord"))
738 0 : continue;
739 :
740 0 : psGCP->dfGCPPixel =
741 0 : CPLAtof(CPLGetXMLValue(psNode, "refColumn", "0.0"));
742 0 : psGCP->dfGCPLine =
743 0 : CPLAtof(CPLGetXMLValue(psNode, "refRow", "0.0"));
744 0 : psGCP->dfGCPX =
745 0 : CPLAtof(CPLGetXMLValue(psNode, "lon", "0.0"));
746 0 : psGCP->dfGCPY =
747 0 : CPLAtof(CPLGetXMLValue(psNode, "lat", "0.0"));
748 0 : psGCP->dfGCPZ = dfAvgHeight;
749 0 : psGCP->pszId = CPLStrdup(CPLSPrintf("%d", nGCP));
750 0 : psGCP->pszInfo = CPLStrdup("");
751 :
752 0 : nGCP++;
753 : }
754 :
755 : // set the projection string - the fields are lat/long - seems
756 : // to be WGS84 datum
757 0 : poDS->m_oGCPSRS.SetWellKnownGeogCS("WGS84");
758 : }
759 : }
760 :
761 : // gcps override geotransform - does it make sense to have both??
762 0 : if (poDS->nGCPCount > 0)
763 : {
764 0 : poDS->bHaveGeoTransform = FALSE;
765 0 : poDS->m_oSRS.Clear();
766 0 : poDS->adfGeoTransform[0] = 0.0;
767 0 : poDS->adfGeoTransform[1] = 1.0;
768 0 : poDS->adfGeoTransform[2] = 0.0;
769 0 : poDS->adfGeoTransform[3] = 0.0;
770 0 : poDS->adfGeoTransform[4] = 0.0;
771 0 : poDS->adfGeoTransform[5] = 1.0;
772 : }
773 : }
774 : else
775 : {
776 0 : CPLError(CE_Warning, CPLE_AppDefined,
777 : "Unable to find sceneInfo tag in XML document. "
778 : "Proceeding with caution.");
779 : }
780 :
781 0 : CPLFree(pszGeorefFile);
782 :
783 : /* -------------------------------------------------------------------- */
784 : /* Initialize any PAM information. */
785 : /* -------------------------------------------------------------------- */
786 0 : poDS->SetDescription(poOpenInfo->pszFilename);
787 0 : poDS->TryLoadXML();
788 :
789 : /* -------------------------------------------------------------------- */
790 : /* Check for overviews. */
791 : /* -------------------------------------------------------------------- */
792 0 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
793 :
794 0 : CPLDestroyXMLNode(psData);
795 :
796 0 : return poDS;
797 : }
798 :
799 : /************************************************************************/
800 : /* GetGCPCount() */
801 : /************************************************************************/
802 :
803 0 : int TSXDataset::GetGCPCount()
804 : {
805 0 : return nGCPCount;
806 : }
807 :
808 : /************************************************************************/
809 : /* GetGCPSpatialRef() */
810 : /************************************************************************/
811 :
812 0 : const OGRSpatialReference *TSXDataset::GetGCPSpatialRef() const
813 : {
814 0 : return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
815 : }
816 :
817 : /************************************************************************/
818 : /* GetGCPs() */
819 : /************************************************************************/
820 :
821 0 : const GDAL_GCP *TSXDataset::GetGCPs()
822 : {
823 0 : return pasGCPList;
824 : }
825 :
826 : /************************************************************************/
827 : /* GetSpatialRef() */
828 : /************************************************************************/
829 :
830 0 : const OGRSpatialReference *TSXDataset::GetSpatialRef() const
831 :
832 : {
833 0 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
834 : }
835 :
836 : /************************************************************************/
837 : /* GetGeotransform() */
838 : /************************************************************************/
839 0 : CPLErr TSXDataset::GetGeoTransform(double *padfTransform)
840 : {
841 0 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
842 :
843 0 : if (bHaveGeoTransform)
844 0 : return CE_None;
845 :
846 0 : return CE_Failure;
847 : }
848 :
849 : /************************************************************************/
850 : /* GDALRegister_TSX() */
851 : /************************************************************************/
852 :
853 1686 : void GDALRegister_TSX()
854 : {
855 1686 : if (GDALGetDriverByName("TSX") != nullptr)
856 302 : return;
857 :
858 1384 : GDALDriver *poDriver = new GDALDriver();
859 :
860 1384 : poDriver->SetDescription("TSX");
861 1384 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
862 1384 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "TerraSAR-X Product");
863 1384 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/tsx.html");
864 1384 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
865 :
866 1384 : poDriver->pfnOpen = TSXDataset::Open;
867 1384 : poDriver->pfnIdentify = TSXDataset::Identify;
868 :
869 1384 : GetGDALDriverManager()->RegisterDriver(poDriver);
870 : }
|
