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