Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Polarimetric Workstation
4 : * Purpose: Radarsat 2 - XML Products (product.xml) driver
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2009-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "cpl_minixml.h"
31 : #include "gdal_frmts.h"
32 : #include "gdal_pam.h"
33 : #include "ogr_spatialref.h"
34 :
35 : typedef enum eCalibration_t
36 : {
37 : Sigma0 = 0,
38 : Gamma,
39 : Beta0,
40 : Uncalib,
41 : None
42 : } eCalibration;
43 :
44 : /*** Function to test for valid LUT files ***/
45 15 : static bool IsValidXMLFile(const char *pszPath, const char *pszLut)
46 : {
47 : /* Return true for valid xml file, false otherwise */
48 15 : char *pszLutFile = VSIStrdup(CPLFormFilename(pszPath, pszLut, nullptr));
49 :
50 15 : CPLXMLTreeCloser psLut(CPLParseXMLFile(pszLutFile));
51 :
52 15 : CPLFree(pszLutFile);
53 :
54 30 : return psLut.get() != nullptr;
55 : }
56 :
57 : /************************************************************************/
58 : /* ==================================================================== */
59 : /* RS2Dataset */
60 : /* ==================================================================== */
61 : /************************************************************************/
62 :
63 : class RS2Dataset final : public GDALPamDataset
64 : {
65 : CPLXMLNode *psProduct;
66 :
67 : int nGCPCount;
68 : GDAL_GCP *pasGCPList;
69 : OGRSpatialReference m_oSRS{};
70 : OGRSpatialReference m_oGCPSRS{};
71 : char **papszSubDatasets;
72 : double adfGeoTransform[6];
73 : bool bHaveGeoTransform;
74 :
75 : char **papszExtraFiles;
76 :
77 : protected:
78 : virtual int CloseDependentDatasets() override;
79 :
80 : public:
81 : RS2Dataset();
82 : virtual ~RS2Dataset();
83 :
84 : virtual int GetGCPCount() override;
85 : const OGRSpatialReference *GetGCPSpatialRef() const override;
86 : virtual const GDAL_GCP *GetGCPs() override;
87 :
88 : const OGRSpatialReference *GetSpatialRef() const override;
89 : virtual CPLErr GetGeoTransform(double *) override;
90 :
91 : virtual char **GetMetadataDomainList() override;
92 : virtual char **GetMetadata(const char *pszDomain = "") override;
93 : virtual char **GetFileList(void) override;
94 :
95 : static GDALDataset *Open(GDALOpenInfo *);
96 : static int Identify(GDALOpenInfo *);
97 :
98 : CPLXMLNode *GetProduct()
99 : {
100 : return psProduct;
101 : }
102 : };
103 :
104 : /************************************************************************/
105 : /* ==================================================================== */
106 : /* RS2RasterBand */
107 : /* ==================================================================== */
108 : /************************************************************************/
109 :
110 : class RS2RasterBand final : public GDALPamRasterBand
111 : {
112 : GDALDataset *poBandFile;
113 :
114 : public:
115 : RS2RasterBand(RS2Dataset *poDSIn, GDALDataType eDataTypeIn,
116 : const char *pszPole, GDALDataset *poBandFile);
117 : virtual ~RS2RasterBand();
118 :
119 : virtual CPLErr IReadBlock(int, int, void *) override;
120 :
121 : static GDALDataset *Open(GDALOpenInfo *);
122 : };
123 :
124 : /************************************************************************/
125 : /* RS2RasterBand */
126 : /************************************************************************/
127 :
128 8 : RS2RasterBand::RS2RasterBand(RS2Dataset *poDSIn, GDALDataType eDataTypeIn,
129 8 : const char *pszPole, GDALDataset *poBandFileIn)
130 8 : : poBandFile(poBandFileIn)
131 : {
132 8 : poDS = poDSIn;
133 :
134 8 : GDALRasterBand *poSrcBand = poBandFile->GetRasterBand(1);
135 :
136 8 : poSrcBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
137 :
138 8 : eDataType = eDataTypeIn;
139 :
140 8 : if (*pszPole != '\0')
141 8 : SetMetadataItem("POLARIMETRIC_INTERP", pszPole);
142 8 : }
143 :
144 : /************************************************************************/
145 : /* RSRasterBand() */
146 : /************************************************************************/
147 :
148 16 : RS2RasterBand::~RS2RasterBand()
149 :
150 : {
151 8 : if (poBandFile != nullptr)
152 8 : GDALClose(reinterpret_cast<GDALRasterBandH>(poBandFile));
153 16 : }
154 :
155 : /************************************************************************/
156 : /* IReadBlock() */
157 : /************************************************************************/
158 :
159 40 : CPLErr RS2RasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
160 :
161 : {
162 : /* -------------------------------------------------------------------- */
163 : /* If the last strip is partial, we need to avoid */
164 : /* over-requesting. We also need to initialize the extra part */
165 : /* of the block to zero. */
166 : /* -------------------------------------------------------------------- */
167 : int nRequestYSize;
168 40 : if ((nBlockYOff + 1) * nBlockYSize > nRasterYSize)
169 : {
170 0 : nRequestYSize = nRasterYSize - nBlockYOff * nBlockYSize;
171 0 : memset(pImage, 0,
172 0 : static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
173 0 : nBlockXSize * nBlockYSize);
174 : }
175 : else
176 : {
177 40 : nRequestYSize = nBlockYSize;
178 : }
179 :
180 : /*-------------------------------------------------------------------- */
181 : /* If the input imagery is tiled, also need to avoid over- */
182 : /* requesting in the X-direction. */
183 : /* ------------------------------------------------------------------- */
184 : int nRequestXSize;
185 40 : if ((nBlockXOff + 1) * nBlockXSize > nRasterXSize)
186 : {
187 0 : nRequestXSize = nRasterXSize - nBlockXOff * nBlockXSize;
188 0 : memset(pImage, 0,
189 0 : static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
190 0 : nBlockXSize * nBlockYSize);
191 : }
192 : else
193 : {
194 40 : nRequestXSize = nBlockXSize;
195 : }
196 40 : if (eDataType == GDT_CInt16 && poBandFile->GetRasterCount() == 2)
197 0 : return poBandFile->RasterIO(
198 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
199 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
200 0 : GDT_Int16, 2, nullptr, 4, nBlockXSize * 4, 2, nullptr);
201 :
202 : /* -------------------------------------------------------------------- */
203 : /* File has one sample marked as sample format void, a 32bits. */
204 : /* -------------------------------------------------------------------- */
205 40 : else if (eDataType == GDT_CInt16 && poBandFile->GetRasterCount() == 1)
206 : {
207 0 : CPLErr eErr = poBandFile->RasterIO(
208 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
209 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
210 0 : GDT_UInt32, 1, nullptr, 4, nBlockXSize * 4, 0, nullptr);
211 :
212 : #ifdef CPL_LSB
213 : /* First, undo the 32bit swap. */
214 0 : GDALSwapWords(pImage, 4, nBlockXSize * nBlockYSize, 4);
215 :
216 : /* Then apply 16 bit swap. */
217 0 : GDALSwapWords(pImage, 2, nBlockXSize * nBlockYSize * 2, 2);
218 : #endif
219 :
220 0 : return eErr;
221 : }
222 :
223 : /* -------------------------------------------------------------------- */
224 : /* The 16bit case is straight forward. The underlying file */
225 : /* looks like a 16bit unsigned data too. */
226 : /* -------------------------------------------------------------------- */
227 40 : else if (eDataType == GDT_UInt16)
228 0 : return poBandFile->RasterIO(
229 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
230 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
231 0 : GDT_UInt16, 1, nullptr, 2, nBlockXSize * 2, 0, nullptr);
232 40 : else if (eDataType == GDT_Byte)
233 : /* Ticket #2104: Support for ScanSAR products */
234 80 : return poBandFile->RasterIO(
235 40 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
236 : nRequestXSize, nRequestYSize, pImage, nRequestXSize, nRequestYSize,
237 40 : GDT_Byte, 1, nullptr, 1, nBlockXSize, 0, nullptr);
238 :
239 0 : CPLAssert(false);
240 : return CE_Failure;
241 : }
242 :
243 : /************************************************************************/
244 : /* ==================================================================== */
245 : /* RS2CalibRasterBand */
246 : /* ==================================================================== */
247 : /************************************************************************/
248 : /* Returns data that has been calibrated to sigma nought, gamma */
249 : /* or beta nought. */
250 : /************************************************************************/
251 :
252 : class RS2CalibRasterBand final : public GDALPamRasterBand
253 : {
254 : private:
255 : // eCalibration m_eCalib;
256 : GDALDataset *m_poBandDataset;
257 : GDALDataType m_eType; /* data type of data being ingested */
258 : float *m_nfTable;
259 : int m_nTableSize;
260 : float m_nfOffset;
261 : char *m_pszLUTFile;
262 :
263 : void ReadLUT();
264 :
265 : public:
266 : RS2CalibRasterBand(RS2Dataset *poDataset, const char *pszPolarization,
267 : GDALDataType eType, GDALDataset *poBandDataset,
268 : eCalibration eCalib, const char *pszLUT);
269 : ~RS2CalibRasterBand();
270 :
271 : CPLErr IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage) override;
272 : };
273 :
274 : /************************************************************************/
275 : /* ReadLUT() */
276 : /************************************************************************/
277 : /* Read the provided LUT in to m_ndTable */
278 : /************************************************************************/
279 2 : void RS2CalibRasterBand::ReadLUT()
280 : {
281 2 : CPLXMLNode *psLUT = CPLParseXMLFile(m_pszLUTFile);
282 :
283 2 : this->m_nfOffset = static_cast<float>(
284 2 : CPLAtof(CPLGetXMLValue(psLUT, "=lut.offset", "0.0")));
285 :
286 2 : char **papszLUTList = CSLTokenizeString2(
287 : CPLGetXMLValue(psLUT, "=lut.gains", ""), " ", CSLT_HONOURSTRINGS);
288 :
289 2 : m_nTableSize = CSLCount(papszLUTList);
290 :
291 2 : m_nfTable =
292 2 : reinterpret_cast<float *>(CPLMalloc(sizeof(float) * m_nTableSize));
293 :
294 514 : for (int i = 0; i < m_nTableSize; i++)
295 : {
296 512 : m_nfTable[i] = static_cast<float>(CPLAtof(papszLUTList[i]));
297 : }
298 :
299 2 : CPLDestroyXMLNode(psLUT);
300 :
301 2 : CSLDestroy(papszLUTList);
302 2 : }
303 :
304 : /************************************************************************/
305 : /* RS2CalibRasterBand() */
306 : /************************************************************************/
307 :
308 2 : RS2CalibRasterBand::RS2CalibRasterBand(
309 : RS2Dataset *poDataset, const char *pszPolarization, GDALDataType eType,
310 2 : GDALDataset *poBandDataset, eCalibration /* eCalib */, const char *pszLUT)
311 : : // m_eCalib(eCalib),
312 : m_poBandDataset(poBandDataset), m_eType(eType), m_nfTable(nullptr),
313 2 : m_nTableSize(0), m_nfOffset(0), m_pszLUTFile(VSIStrdup(pszLUT))
314 : {
315 2 : poDS = poDataset;
316 :
317 2 : if (*pszPolarization != '\0')
318 : {
319 2 : SetMetadataItem("POLARIMETRIC_INTERP", pszPolarization);
320 : }
321 :
322 2 : if (eType == GDT_CInt16)
323 0 : eDataType = GDT_CFloat32;
324 : else
325 2 : eDataType = GDT_Float32;
326 :
327 2 : GDALRasterBand *poRasterBand = poBandDataset->GetRasterBand(1);
328 2 : poRasterBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
329 :
330 2 : ReadLUT();
331 2 : }
332 :
333 : /************************************************************************/
334 : /* ~RS2CalibRasterBand() */
335 : /************************************************************************/
336 :
337 4 : RS2CalibRasterBand::~RS2CalibRasterBand()
338 : {
339 2 : CPLFree(m_nfTable);
340 2 : CPLFree(m_pszLUTFile);
341 2 : GDALClose(m_poBandDataset);
342 4 : }
343 :
344 : /************************************************************************/
345 : /* IReadBlock() */
346 : /************************************************************************/
347 :
348 20 : CPLErr RS2CalibRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
349 : void *pImage)
350 : {
351 : /* -------------------------------------------------------------------- */
352 : /* If the last strip is partial, we need to avoid */
353 : /* over-requesting. We also need to initialize the extra part */
354 : /* of the block to zero. */
355 : /* -------------------------------------------------------------------- */
356 : int nRequestYSize;
357 20 : if ((nBlockYOff + 1) * nBlockYSize > nRasterYSize)
358 : {
359 0 : nRequestYSize = nRasterYSize - nBlockYOff * nBlockYSize;
360 0 : memset(pImage, 0,
361 0 : static_cast<size_t>(GDALGetDataTypeSizeBytes(eDataType)) *
362 0 : nBlockXSize * nBlockYSize);
363 : }
364 : else
365 : {
366 20 : nRequestYSize = nBlockYSize;
367 : }
368 :
369 : CPLErr eErr;
370 20 : if (m_eType == GDT_CInt16)
371 : {
372 : /* read in complex values */
373 : GInt16 *pnImageTmp = reinterpret_cast<GInt16 *>(
374 0 : CPLMalloc(2 * nBlockXSize * nBlockYSize *
375 0 : GDALGetDataTypeSize(GDT_Int16) / 8));
376 0 : if (m_poBandDataset->GetRasterCount() == 2)
377 : {
378 0 : eErr = m_poBandDataset->RasterIO(
379 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
380 : nBlockXSize, nRequestYSize, pnImageTmp, nBlockXSize,
381 0 : nRequestYSize, GDT_Int16, 2, nullptr, 4, nBlockXSize * 4, 2,
382 : nullptr);
383 : }
384 : else
385 : {
386 0 : eErr = m_poBandDataset->RasterIO(
387 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
388 : nBlockXSize, nRequestYSize, pnImageTmp, nBlockXSize,
389 0 : nRequestYSize, GDT_UInt32, 1, nullptr, 4, nBlockXSize * 4, 0,
390 : nullptr);
391 :
392 : #ifdef CPL_LSB
393 : /* First, undo the 32bit swap. */
394 0 : GDALSwapWords(pImage, 4, nBlockXSize * nBlockYSize, 4);
395 :
396 : /* Then apply 16 bit swap. */
397 0 : GDALSwapWords(pImage, 2, nBlockXSize * nBlockYSize * 2, 2);
398 : #endif
399 : }
400 :
401 : /* calibrate the complex values */
402 0 : for (int i = 0; i < nBlockYSize; i++)
403 : {
404 0 : for (int j = 0; j < nBlockXSize; j++)
405 : {
406 : /* calculate pixel offset in memory*/
407 0 : int nPixOff = (2 * (i * nBlockXSize)) + (j * 2);
408 :
409 0 : reinterpret_cast<float *>(pImage)[nPixOff] =
410 0 : static_cast<float>(pnImageTmp[nPixOff]) /
411 0 : (m_nfTable[nBlockXOff + j]);
412 0 : reinterpret_cast<float *>(pImage)[nPixOff + 1] =
413 0 : static_cast<float>(pnImageTmp[nPixOff + 1]) /
414 0 : (m_nfTable[nBlockXOff + j]);
415 : }
416 : }
417 0 : CPLFree(pnImageTmp);
418 : }
419 20 : else if (m_eType == GDT_UInt16)
420 : {
421 : /* read in detected values */
422 0 : GUInt16 *pnImageTmp = reinterpret_cast<GUInt16 *>(CPLMalloc(
423 0 : nBlockXSize * nBlockYSize * GDALGetDataTypeSize(GDT_UInt16) / 8));
424 0 : eErr = m_poBandDataset->RasterIO(
425 0 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
426 : nBlockXSize, nRequestYSize, pnImageTmp, nBlockXSize, nRequestYSize,
427 0 : GDT_UInt16, 1, nullptr, 2, nBlockXSize * 2, 0, nullptr);
428 :
429 : /* iterate over detected values */
430 0 : for (int i = 0; i < nBlockYSize; i++)
431 : {
432 0 : for (int j = 0; j < nBlockXSize; j++)
433 : {
434 0 : int nPixOff = (i * nBlockXSize) + j;
435 :
436 0 : reinterpret_cast<float *>(pImage)[nPixOff] =
437 0 : ((static_cast<float>(pnImageTmp[nPixOff]) *
438 0 : static_cast<float>(pnImageTmp[nPixOff])) +
439 0 : m_nfOffset) /
440 0 : m_nfTable[nBlockXOff + j];
441 : }
442 : }
443 0 : CPLFree(pnImageTmp);
444 : } /* Ticket #2104: Support for ScanSAR products */
445 20 : else if (m_eType == GDT_Byte)
446 : {
447 40 : GByte *pnImageTmp = reinterpret_cast<GByte *>(CPLMalloc(
448 20 : nBlockXSize * nBlockYSize * GDALGetDataTypeSize(GDT_Byte) / 8));
449 40 : eErr = m_poBandDataset->RasterIO(
450 20 : GF_Read, nBlockXOff * nBlockXSize, nBlockYOff * nBlockYSize,
451 : nBlockXSize, nRequestYSize, pnImageTmp, nBlockXSize, nRequestYSize,
452 : GDT_Byte, 1, nullptr, 1, 1, 0, nullptr);
453 :
454 : /* iterate over detected values */
455 40 : for (int i = 0; i < nBlockYSize; i++)
456 : {
457 420 : for (int j = 0; j < nBlockXSize; j++)
458 : {
459 400 : int nPixOff = (i * nBlockXSize) + j;
460 :
461 400 : reinterpret_cast<float *>(pImage)[nPixOff] =
462 400 : ((pnImageTmp[nPixOff] * pnImageTmp[nPixOff]) + m_nfOffset) /
463 400 : m_nfTable[nBlockXOff + j];
464 : }
465 : }
466 20 : CPLFree(pnImageTmp);
467 : }
468 : else
469 : {
470 0 : CPLAssert(false);
471 : return CE_Failure;
472 : }
473 20 : return eErr;
474 : }
475 :
476 : /************************************************************************/
477 : /* ==================================================================== */
478 : /* RS2Dataset */
479 : /* ==================================================================== */
480 : /************************************************************************/
481 :
482 : /************************************************************************/
483 : /* RS2Dataset() */
484 : /************************************************************************/
485 :
486 5 : RS2Dataset::RS2Dataset()
487 : : psProduct(nullptr), nGCPCount(0), pasGCPList(nullptr),
488 : papszSubDatasets(nullptr), bHaveGeoTransform(FALSE),
489 5 : papszExtraFiles(nullptr)
490 : {
491 5 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
492 5 : m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
493 5 : adfGeoTransform[0] = 0.0;
494 5 : adfGeoTransform[1] = 1.0;
495 5 : adfGeoTransform[2] = 0.0;
496 5 : adfGeoTransform[3] = 0.0;
497 5 : adfGeoTransform[4] = 0.0;
498 5 : adfGeoTransform[5] = 1.0;
499 5 : }
500 :
501 : /************************************************************************/
502 : /* ~RS2Dataset() */
503 : /************************************************************************/
504 :
505 10 : RS2Dataset::~RS2Dataset()
506 :
507 : {
508 5 : RS2Dataset::FlushCache(true);
509 :
510 5 : CPLDestroyXMLNode(psProduct);
511 :
512 5 : if (nGCPCount > 0)
513 : {
514 5 : GDALDeinitGCPs(nGCPCount, pasGCPList);
515 5 : CPLFree(pasGCPList);
516 : }
517 :
518 5 : RS2Dataset::CloseDependentDatasets();
519 :
520 5 : CSLDestroy(papszSubDatasets);
521 5 : CSLDestroy(papszExtraFiles);
522 10 : }
523 :
524 : /************************************************************************/
525 : /* CloseDependentDatasets() */
526 : /************************************************************************/
527 :
528 11 : int RS2Dataset::CloseDependentDatasets()
529 : {
530 11 : int bHasDroppedRef = GDALPamDataset::CloseDependentDatasets();
531 :
532 11 : if (nBands != 0)
533 5 : bHasDroppedRef = TRUE;
534 :
535 21 : for (int iBand = 0; iBand < nBands; iBand++)
536 : {
537 10 : delete papoBands[iBand];
538 : }
539 11 : nBands = 0;
540 :
541 11 : return bHasDroppedRef;
542 : }
543 :
544 : /************************************************************************/
545 : /* GetFileList() */
546 : /************************************************************************/
547 :
548 2 : char **RS2Dataset::GetFileList()
549 :
550 : {
551 2 : char **papszFileList = GDALPamDataset::GetFileList();
552 :
553 2 : papszFileList = CSLInsertStrings(papszFileList, -1, papszExtraFiles);
554 :
555 2 : return papszFileList;
556 : }
557 :
558 : /************************************************************************/
559 : /* Identify() */
560 : /************************************************************************/
561 :
562 52814 : int RS2Dataset::Identify(GDALOpenInfo *poOpenInfo)
563 : {
564 : /* Check for the case where we're trying to read the calibrated data: */
565 52814 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "RADARSAT_2_CALIB:"))
566 : {
567 2 : return TRUE;
568 : }
569 :
570 : /* Check for directory access when there is a product.xml file in the
571 : directory. */
572 52812 : if (poOpenInfo->bIsDirectory)
573 : {
574 : CPLString osMDFilename =
575 724 : CPLFormCIFilename(poOpenInfo->pszFilename, "product.xml", nullptr);
576 :
577 : VSIStatBufL sStat;
578 362 : if (VSIStatL(osMDFilename, &sStat) == 0)
579 0 : return TRUE;
580 :
581 362 : return FALSE;
582 : }
583 :
584 : /* otherwise, do our normal stuff */
585 52450 : if (strlen(poOpenInfo->pszFilename) < 11 ||
586 51421 : !EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename) - 11,
587 : "product.xml"))
588 52442 : return FALSE;
589 :
590 8 : if (poOpenInfo->nHeaderBytes < 100)
591 0 : return FALSE;
592 :
593 8 : if (strstr((const char *)poOpenInfo->pabyHeader, "/rs2") == nullptr ||
594 8 : strstr((const char *)poOpenInfo->pabyHeader, "<product") == nullptr)
595 0 : return FALSE;
596 :
597 8 : return TRUE;
598 : }
599 :
600 : /************************************************************************/
601 : /* Open() */
602 : /************************************************************************/
603 :
604 5 : GDALDataset *RS2Dataset::Open(GDALOpenInfo *poOpenInfo)
605 :
606 : {
607 : /* -------------------------------------------------------------------- */
608 : /* Is this a RADARSAT-2 Product.xml definition? */
609 : /* -------------------------------------------------------------------- */
610 5 : if (!RS2Dataset::Identify(poOpenInfo))
611 : {
612 0 : return nullptr;
613 : }
614 :
615 : /* -------------------------------------------------------------------- */
616 : /* Get subdataset information, if relevant */
617 : /* -------------------------------------------------------------------- */
618 5 : const char *pszFilename = poOpenInfo->pszFilename;
619 5 : eCalibration eCalib = None;
620 :
621 5 : if (STARTS_WITH_CI(pszFilename, "RADARSAT_2_CALIB:"))
622 : {
623 1 : pszFilename += 17;
624 :
625 1 : if (STARTS_WITH_CI(pszFilename, "BETA0"))
626 1 : eCalib = Beta0;
627 0 : else if (STARTS_WITH_CI(pszFilename, "SIGMA0"))
628 0 : eCalib = Sigma0;
629 0 : else if (STARTS_WITH_CI(pszFilename, "GAMMA"))
630 0 : eCalib = Gamma;
631 0 : else if (STARTS_WITH_CI(pszFilename, "UNCALIB"))
632 0 : eCalib = Uncalib;
633 : else
634 0 : eCalib = None;
635 :
636 : /* advance the pointer to the actual filename */
637 6 : while (*pszFilename != '\0' && *pszFilename != ':')
638 5 : pszFilename++;
639 :
640 1 : if (*pszFilename == ':')
641 1 : pszFilename++;
642 :
643 : // need to redo the directory check:
644 : // the GDALOpenInfo check would have failed because of the calibration
645 : // string on the filename
646 : VSIStatBufL sStat;
647 1 : if (VSIStatL(pszFilename, &sStat) == 0)
648 1 : poOpenInfo->bIsDirectory = VSI_ISDIR(sStat.st_mode);
649 : }
650 :
651 10 : CPLString osMDFilename;
652 5 : if (poOpenInfo->bIsDirectory)
653 : {
654 0 : osMDFilename = CPLFormCIFilename(pszFilename, "product.xml", nullptr);
655 : }
656 : else
657 5 : osMDFilename = pszFilename;
658 :
659 : /* -------------------------------------------------------------------- */
660 : /* Ingest the Product.xml file. */
661 : /* -------------------------------------------------------------------- */
662 5 : CPLXMLNode *psProduct = CPLParseXMLFile(osMDFilename);
663 5 : if (psProduct == nullptr)
664 0 : return nullptr;
665 :
666 : /* -------------------------------------------------------------------- */
667 : /* Confirm the requested access is supported. */
668 : /* -------------------------------------------------------------------- */
669 5 : if (poOpenInfo->eAccess == GA_Update)
670 : {
671 0 : CPLDestroyXMLNode(psProduct);
672 0 : CPLError(CE_Failure, CPLE_NotSupported,
673 : "The RS2 driver does not support update access to existing"
674 : " datasets.\n");
675 0 : return nullptr;
676 : }
677 :
678 : CPLXMLNode *psImageAttributes =
679 5 : CPLGetXMLNode(psProduct, "=product.imageAttributes");
680 5 : if (psImageAttributes == nullptr)
681 : {
682 0 : CPLDestroyXMLNode(psProduct);
683 0 : CPLError(CE_Failure, CPLE_OpenFailed,
684 : "Failed to find <imageAttributes> in document.");
685 0 : return nullptr;
686 : }
687 :
688 : CPLXMLNode *psImageGenerationParameters =
689 5 : CPLGetXMLNode(psProduct, "=product.imageGenerationParameters");
690 5 : if (psImageGenerationParameters == nullptr)
691 : {
692 0 : CPLDestroyXMLNode(psProduct);
693 0 : CPLError(CE_Failure, CPLE_OpenFailed,
694 : "Failed to find <imageGenerationParameters> in document.");
695 0 : return nullptr;
696 : }
697 :
698 : /* -------------------------------------------------------------------- */
699 : /* Create the dataset. */
700 : /* -------------------------------------------------------------------- */
701 5 : RS2Dataset *poDS = new RS2Dataset();
702 :
703 5 : poDS->psProduct = psProduct;
704 :
705 : /* -------------------------------------------------------------------- */
706 : /* Get overall image information. */
707 : /* -------------------------------------------------------------------- */
708 5 : poDS->nRasterXSize = atoi(CPLGetXMLValue(
709 : psImageAttributes, "rasterAttributes.numberOfSamplesPerLine", "-1"));
710 5 : poDS->nRasterYSize = atoi(CPLGetXMLValue(
711 : psImageAttributes, "rasterAttributes.numberofLines", "-1"));
712 5 : if (poDS->nRasterXSize <= 1 || poDS->nRasterYSize <= 1)
713 : {
714 0 : CPLError(
715 : CE_Failure, CPLE_OpenFailed,
716 : "Non-sane raster dimensions provided in product.xml. If this is "
717 : "a valid RADARSAT-2 scene, please contact your data provider for "
718 : "a corrected dataset.");
719 0 : delete poDS;
720 0 : return nullptr;
721 : }
722 :
723 : /* -------------------------------------------------------------------- */
724 : /* Check product type, as to determine if there are LUTs for */
725 : /* calibration purposes. */
726 : /* -------------------------------------------------------------------- */
727 :
728 : const char *pszProductType =
729 5 : CPLGetXMLValue(psImageGenerationParameters,
730 : "generalProcessingInformation.productType", "UNK");
731 :
732 5 : poDS->SetMetadataItem("PRODUCT_TYPE", pszProductType);
733 :
734 : /* the following cases can be assumed to have no LUTs, as per
735 : * RN-RP-51-2713, but also common sense
736 : */
737 5 : bool bCanCalib = false;
738 5 : if (!(STARTS_WITH_CI(pszProductType, "UNK") ||
739 5 : STARTS_WITH_CI(pszProductType, "SSG") ||
740 5 : STARTS_WITH_CI(pszProductType, "SPG")))
741 : {
742 5 : bCanCalib = true;
743 : }
744 :
745 : /* -------------------------------------------------------------------- */
746 : /* Get dataType (so we can recognise complex data), and the */
747 : /* bitsPerSample. */
748 : /* -------------------------------------------------------------------- */
749 : const char *pszDataType =
750 5 : CPLGetXMLValue(psImageAttributes, "rasterAttributes.dataType", "");
751 5 : const int nBitsPerSample = atoi(CPLGetXMLValue(
752 : psImageAttributes, "rasterAttributes.bitsPerSample", ""));
753 :
754 : GDALDataType eDataType;
755 5 : if (nBitsPerSample == 16 && EQUAL(pszDataType, "Complex"))
756 0 : eDataType = GDT_CInt16;
757 5 : else if (nBitsPerSample == 16 && STARTS_WITH_CI(pszDataType, "Mag"))
758 0 : eDataType = GDT_UInt16;
759 5 : else if (nBitsPerSample == 8 && STARTS_WITH_CI(pszDataType, "Mag"))
760 5 : eDataType = GDT_Byte;
761 : else
762 : {
763 0 : delete poDS;
764 0 : CPLError(
765 : CE_Failure, CPLE_AppDefined,
766 : "dataType=%s, bitsPerSample=%d: not a supported configuration.",
767 : pszDataType, nBitsPerSample);
768 0 : return nullptr;
769 : }
770 :
771 : /* while we're at it, extract the pixel spacing information */
772 5 : const char *pszPixelSpacing = CPLGetXMLValue(
773 : psImageAttributes, "rasterAttributes.sampledPixelSpacing", "UNK");
774 5 : poDS->SetMetadataItem("PIXEL_SPACING", pszPixelSpacing);
775 :
776 5 : const char *pszLineSpacing = CPLGetXMLValue(
777 : psImageAttributes, "rasterAttributes.sampledLineSpacing", "UNK");
778 5 : poDS->SetMetadataItem("LINE_SPACING", pszLineSpacing);
779 :
780 : /* -------------------------------------------------------------------- */
781 : /* Open each of the data files as a complex band. */
782 : /* -------------------------------------------------------------------- */
783 10 : CPLString osBeta0LUT;
784 10 : CPLString osGammaLUT;
785 10 : CPLString osSigma0LUT;
786 :
787 5 : char *pszPath = CPLStrdup(CPLGetPath(osMDFilename));
788 5 : const int nFLen = static_cast<int>(osMDFilename.size());
789 :
790 5 : CPLXMLNode *psNode = psImageAttributes->psChild;
791 40 : for (; psNode != nullptr; psNode = psNode->psNext)
792 : {
793 35 : if (psNode->eType != CXT_Element ||
794 35 : !(EQUAL(psNode->pszValue, "fullResolutionImageData") ||
795 25 : EQUAL(psNode->pszValue, "lookupTable")))
796 10 : continue;
797 :
798 25 : if (EQUAL(psNode->pszValue, "lookupTable") && bCanCalib)
799 : {
800 : /* Determine which incidence angle correction this is */
801 : const char *pszLUTType =
802 15 : CPLGetXMLValue(psNode, "incidenceAngleCorrection", "");
803 15 : const char *pszLUTFile = CPLGetXMLValue(psNode, "", "");
804 : CPLString osLUTFilePath =
805 30 : CPLFormFilename(pszPath, pszLUTFile, nullptr);
806 :
807 15 : if (EQUAL(pszLUTType, ""))
808 0 : continue;
809 20 : else if (EQUAL(pszLUTType, "Beta Nought") &&
810 5 : IsValidXMLFile(pszPath, pszLUTFile))
811 : {
812 5 : poDS->papszExtraFiles =
813 5 : CSLAddString(poDS->papszExtraFiles, osLUTFilePath);
814 :
815 5 : const size_t nBufLen = nFLen + 27;
816 5 : char *pszBuf = reinterpret_cast<char *>(CPLMalloc(nBufLen));
817 5 : osBeta0LUT = pszLUTFile;
818 5 : poDS->SetMetadataItem("BETA_NOUGHT_LUT", pszLUTFile);
819 :
820 5 : snprintf(pszBuf, nBufLen, "RADARSAT_2_CALIB:BETA0:%s",
821 : osMDFilename.c_str());
822 5 : poDS->papszSubDatasets = CSLSetNameValue(
823 : poDS->papszSubDatasets, "SUBDATASET_3_NAME", pszBuf);
824 5 : poDS->papszSubDatasets =
825 5 : CSLSetNameValue(poDS->papszSubDatasets, "SUBDATASET_3_DESC",
826 : "Beta Nought calibrated");
827 5 : CPLFree(pszBuf);
828 : }
829 15 : else if (EQUAL(pszLUTType, "Sigma Nought") &&
830 5 : IsValidXMLFile(pszPath, pszLUTFile))
831 : {
832 5 : poDS->papszExtraFiles =
833 5 : CSLAddString(poDS->papszExtraFiles, osLUTFilePath);
834 :
835 5 : const size_t nBufLen = nFLen + 27;
836 5 : char *pszBuf = reinterpret_cast<char *>(CPLMalloc(nBufLen));
837 5 : osSigma0LUT = pszLUTFile;
838 5 : poDS->SetMetadataItem("SIGMA_NOUGHT_LUT", pszLUTFile);
839 :
840 5 : snprintf(pszBuf, nBufLen, "RADARSAT_2_CALIB:SIGMA0:%s",
841 : osMDFilename.c_str());
842 5 : poDS->papszSubDatasets = CSLSetNameValue(
843 : poDS->papszSubDatasets, "SUBDATASET_2_NAME", pszBuf);
844 5 : poDS->papszSubDatasets =
845 5 : CSLSetNameValue(poDS->papszSubDatasets, "SUBDATASET_2_DESC",
846 : "Sigma Nought calibrated");
847 5 : CPLFree(pszBuf);
848 : }
849 10 : else if (EQUAL(pszLUTType, "Gamma") &&
850 5 : IsValidXMLFile(pszPath, pszLUTFile))
851 : {
852 5 : poDS->papszExtraFiles =
853 5 : CSLAddString(poDS->papszExtraFiles, osLUTFilePath);
854 :
855 5 : const size_t nBufLen = nFLen + 27;
856 5 : char *pszBuf = reinterpret_cast<char *>(CPLMalloc(nBufLen));
857 5 : osGammaLUT = pszLUTFile;
858 5 : poDS->SetMetadataItem("GAMMA_LUT", pszLUTFile);
859 5 : snprintf(pszBuf, nBufLen, "RADARSAT_2_CALIB:GAMMA:%s",
860 : osMDFilename.c_str());
861 5 : poDS->papszSubDatasets = CSLSetNameValue(
862 : poDS->papszSubDatasets, "SUBDATASET_4_NAME", pszBuf);
863 5 : poDS->papszSubDatasets =
864 5 : CSLSetNameValue(poDS->papszSubDatasets, "SUBDATASET_4_DESC",
865 : "Gamma calibrated");
866 5 : CPLFree(pszBuf);
867 : }
868 15 : continue;
869 : }
870 :
871 : /* --------------------------------------------------------------------
872 : */
873 : /* Fetch filename. */
874 : /* --------------------------------------------------------------------
875 : */
876 10 : const char *pszBasename = CPLGetXMLValue(psNode, "", "");
877 10 : if (*pszBasename == '\0')
878 0 : continue;
879 :
880 : /* --------------------------------------------------------------------
881 : */
882 : /* Form full filename (path of product.xml + basename). */
883 : /* --------------------------------------------------------------------
884 : */
885 : char *pszFullname =
886 10 : CPLStrdup(CPLFormFilename(pszPath, pszBasename, nullptr));
887 :
888 : /* --------------------------------------------------------------------
889 : */
890 : /* Try and open the file. */
891 : /* --------------------------------------------------------------------
892 : */
893 : GDALDataset *poBandFile =
894 10 : GDALDataset::FromHandle(GDALOpen(pszFullname, GA_ReadOnly));
895 10 : if (poBandFile == nullptr)
896 : {
897 0 : CPLFree(pszFullname);
898 0 : continue;
899 : }
900 10 : if (poBandFile->GetRasterCount() == 0)
901 : {
902 0 : GDALClose(reinterpret_cast<GDALRasterBandH>(poBandFile));
903 0 : CPLFree(pszFullname);
904 0 : continue;
905 : }
906 :
907 10 : poDS->papszExtraFiles =
908 10 : CSLAddString(poDS->papszExtraFiles, pszFullname);
909 :
910 : /* --------------------------------------------------------------------
911 : */
912 : /* Create the band. */
913 : /* --------------------------------------------------------------------
914 : */
915 10 : if (eCalib == None || eCalib == Uncalib)
916 : {
917 : RS2RasterBand *poBand = new RS2RasterBand(
918 8 : poDS, eDataType, CPLGetXMLValue(psNode, "pole", ""),
919 8 : poBandFile);
920 :
921 8 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand);
922 : }
923 : else
924 : {
925 2 : const char *pszLUT = nullptr;
926 2 : switch (eCalib)
927 : {
928 0 : case Sigma0:
929 0 : pszLUT = osSigma0LUT;
930 0 : break;
931 2 : case Beta0:
932 2 : pszLUT = osBeta0LUT;
933 2 : break;
934 0 : case Gamma:
935 0 : pszLUT = osGammaLUT;
936 0 : break;
937 0 : default:
938 : /* we should bomb gracefully... */
939 0 : pszLUT = osSigma0LUT;
940 : }
941 : RS2CalibRasterBand *poBand = new RS2CalibRasterBand(
942 2 : poDS, CPLGetXMLValue(psNode, "pole", ""), eDataType, poBandFile,
943 2 : eCalib, CPLFormFilename(pszPath, pszLUT, nullptr));
944 2 : poDS->SetBand(poDS->GetRasterCount() + 1, poBand);
945 : }
946 :
947 10 : CPLFree(pszFullname);
948 : }
949 :
950 5 : if (poDS->papszSubDatasets != nullptr && eCalib == None)
951 : {
952 4 : const size_t nBufLen = nFLen + 28;
953 4 : char *pszBuf = reinterpret_cast<char *>(CPLMalloc(nBufLen));
954 4 : snprintf(pszBuf, nBufLen, "RADARSAT_2_CALIB:UNCALIB:%s",
955 : osMDFilename.c_str());
956 4 : poDS->papszSubDatasets = CSLSetNameValue(poDS->papszSubDatasets,
957 : "SUBDATASET_1_NAME", pszBuf);
958 4 : poDS->papszSubDatasets =
959 4 : CSLSetNameValue(poDS->papszSubDatasets, "SUBDATASET_1_DESC",
960 : "Uncalibrated digital numbers");
961 4 : CPLFree(pszBuf);
962 : }
963 1 : else if (poDS->papszSubDatasets != nullptr)
964 : {
965 1 : CSLDestroy(poDS->papszSubDatasets);
966 1 : poDS->papszSubDatasets = nullptr;
967 : }
968 :
969 : /* -------------------------------------------------------------------- */
970 : /* Set the appropriate MATRIX_REPRESENTATION. */
971 : /* -------------------------------------------------------------------- */
972 :
973 5 : if (poDS->GetRasterCount() == 4 &&
974 0 : (eDataType == GDT_CInt16 || eDataType == GDT_CFloat32))
975 : {
976 0 : poDS->SetMetadataItem("MATRIX_REPRESENTATION", "SCATTERING");
977 : }
978 :
979 : /* -------------------------------------------------------------------- */
980 : /* Collect a few useful metadata items */
981 : /* -------------------------------------------------------------------- */
982 :
983 : CPLXMLNode *psSourceAttrs =
984 5 : CPLGetXMLNode(psProduct, "=product.sourceAttributes");
985 5 : const char *pszItem = CPLGetXMLValue(psSourceAttrs, "satellite", "");
986 5 : poDS->SetMetadataItem("SATELLITE_IDENTIFIER", pszItem);
987 :
988 5 : pszItem = CPLGetXMLValue(psSourceAttrs, "sensor", "");
989 5 : poDS->SetMetadataItem("SENSOR_IDENTIFIER", pszItem);
990 :
991 5 : if (psSourceAttrs != nullptr)
992 : {
993 : /* Get beam mode mnemonic */
994 5 : pszItem = CPLGetXMLValue(psSourceAttrs, "beamModeMnemonic", "UNK");
995 5 : poDS->SetMetadataItem("BEAM_MODE", pszItem);
996 5 : pszItem = CPLGetXMLValue(psSourceAttrs, "rawDataStartTime", "UNK");
997 5 : poDS->SetMetadataItem("ACQUISITION_START_TIME", pszItem);
998 :
999 : pszItem =
1000 5 : CPLGetXMLValue(psSourceAttrs, "inputDatasetFacilityId", "UNK");
1001 5 : poDS->SetMetadataItem("FACILITY_IDENTIFIER", pszItem);
1002 :
1003 5 : pszItem = CPLGetXMLValue(
1004 : psSourceAttrs, "orbitAndAttitude.orbitInformation.passDirection",
1005 : "UNK");
1006 5 : poDS->SetMetadataItem("ORBIT_DIRECTION", pszItem);
1007 5 : pszItem = CPLGetXMLValue(
1008 : psSourceAttrs, "orbitAndAttitude.orbitInformation.orbitDataSource",
1009 : "UNK");
1010 5 : poDS->SetMetadataItem("ORBIT_DATA_SOURCE", pszItem);
1011 5 : pszItem = CPLGetXMLValue(
1012 : psSourceAttrs, "orbitAndAttitude.orbitInformation.orbitDataFile",
1013 : "UNK");
1014 5 : poDS->SetMetadataItem("ORBIT_DATA_FILE", pszItem);
1015 : }
1016 :
1017 : CPLXMLNode *psSarProcessingInformation =
1018 5 : CPLGetXMLNode(psProduct, "=product.imageGenerationParameters");
1019 :
1020 5 : if (psSarProcessingInformation != nullptr)
1021 : {
1022 : /* Get incidence angle information */
1023 5 : pszItem = CPLGetXMLValue(
1024 : psSarProcessingInformation,
1025 : "sarProcessingInformation.incidenceAngleNearRange", "UNK");
1026 5 : poDS->SetMetadataItem("NEAR_RANGE_INCIDENCE_ANGLE", pszItem);
1027 :
1028 5 : pszItem = CPLGetXMLValue(
1029 : psSarProcessingInformation,
1030 : "sarProcessingInformation.incidenceAngleFarRange", "UNK");
1031 5 : poDS->SetMetadataItem("FAR_RANGE_INCIDENCE_ANGLE", pszItem);
1032 :
1033 5 : pszItem = CPLGetXMLValue(psSarProcessingInformation,
1034 : "sarProcessingInformation.slantRangeNearEdge",
1035 : "UNK");
1036 5 : poDS->SetMetadataItem("SLANT_RANGE_NEAR_EDGE", pszItem);
1037 :
1038 5 : pszItem = CPLGetXMLValue(
1039 : psSarProcessingInformation,
1040 : "sarProcessingInformation.zeroDopplerTimeFirstLine", "UNK");
1041 5 : poDS->SetMetadataItem("FIRST_LINE_TIME", pszItem);
1042 :
1043 5 : pszItem = CPLGetXMLValue(
1044 : psSarProcessingInformation,
1045 : "sarProcessingInformation.zeroDopplerTimeLastLine", "UNK");
1046 5 : poDS->SetMetadataItem("LAST_LINE_TIME", pszItem);
1047 :
1048 : pszItem =
1049 5 : CPLGetXMLValue(psSarProcessingInformation,
1050 : "generalProcessingInformation.productType", "UNK");
1051 5 : poDS->SetMetadataItem("PRODUCT_TYPE", pszItem);
1052 :
1053 5 : pszItem = CPLGetXMLValue(
1054 : psSarProcessingInformation,
1055 : "generalProcessingInformation.processingFacility", "UNK");
1056 5 : poDS->SetMetadataItem("PROCESSING_FACILITY", pszItem);
1057 :
1058 5 : pszItem = CPLGetXMLValue(psSarProcessingInformation,
1059 : "generalProcessingInformation.processingTime",
1060 : "UNK");
1061 5 : poDS->SetMetadataItem("PROCESSING_TIME", pszItem);
1062 : }
1063 :
1064 : /*--------------------------------------------------------------------- */
1065 : /* Collect Map projection/Geotransform information, if present */
1066 : /* -------------------------------------------------------------------- */
1067 : CPLXMLNode *psMapProjection =
1068 5 : CPLGetXMLNode(psImageAttributes, "geographicInformation.mapProjection");
1069 :
1070 5 : if (psMapProjection != nullptr)
1071 : {
1072 : CPLXMLNode *psPos =
1073 0 : CPLGetXMLNode(psMapProjection, "positioningInformation");
1074 :
1075 : pszItem =
1076 0 : CPLGetXMLValue(psMapProjection, "mapProjectionDescriptor", "UNK");
1077 0 : poDS->SetMetadataItem("MAP_PROJECTION_DESCRIPTOR", pszItem);
1078 :
1079 : pszItem =
1080 0 : CPLGetXMLValue(psMapProjection, "mapProjectionOrientation", "UNK");
1081 0 : poDS->SetMetadataItem("MAP_PROJECTION_ORIENTATION", pszItem);
1082 :
1083 0 : pszItem = CPLGetXMLValue(psMapProjection, "resamplingKernel", "UNK");
1084 0 : poDS->SetMetadataItem("RESAMPLING_KERNEL", pszItem);
1085 :
1086 0 : pszItem = CPLGetXMLValue(psMapProjection, "satelliteHeading", "UNK");
1087 0 : poDS->SetMetadataItem("SATELLITE_HEADING", pszItem);
1088 :
1089 0 : if (psPos != nullptr)
1090 : {
1091 0 : const double tl_x = CPLStrtod(
1092 : CPLGetXMLValue(psPos, "upperLeftCorner.mapCoordinate.easting",
1093 : "0.0"),
1094 : nullptr);
1095 0 : const double tl_y = CPLStrtod(
1096 : CPLGetXMLValue(psPos, "upperLeftCorner.mapCoordinate.northing",
1097 : "0.0"),
1098 : nullptr);
1099 0 : const double bl_x = CPLStrtod(
1100 : CPLGetXMLValue(psPos, "lowerLeftCorner.mapCoordinate.easting",
1101 : "0.0"),
1102 : nullptr);
1103 0 : const double bl_y = CPLStrtod(
1104 : CPLGetXMLValue(psPos, "lowerLeftCorner.mapCoordinate.northing",
1105 : "0.0"),
1106 : nullptr);
1107 0 : const double tr_x = CPLStrtod(
1108 : CPLGetXMLValue(psPos, "upperRightCorner.mapCoordinate.easting",
1109 : "0.0"),
1110 : nullptr);
1111 0 : const double tr_y = CPLStrtod(
1112 : CPLGetXMLValue(psPos, "upperRightCorner.mapCoordinate.northing",
1113 : "0.0"),
1114 : nullptr);
1115 0 : poDS->adfGeoTransform[1] = (tr_x - tl_x) / (poDS->nRasterXSize - 1);
1116 0 : poDS->adfGeoTransform[4] = (tr_y - tl_y) / (poDS->nRasterXSize - 1);
1117 0 : poDS->adfGeoTransform[2] = (bl_x - tl_x) / (poDS->nRasterYSize - 1);
1118 0 : poDS->adfGeoTransform[5] = (bl_y - tl_y) / (poDS->nRasterYSize - 1);
1119 0 : poDS->adfGeoTransform[0] = (tl_x - 0.5 * poDS->adfGeoTransform[1] -
1120 0 : 0.5 * poDS->adfGeoTransform[2]);
1121 0 : poDS->adfGeoTransform[3] = (tl_y - 0.5 * poDS->adfGeoTransform[4] -
1122 0 : 0.5 * poDS->adfGeoTransform[5]);
1123 :
1124 : /* Use bottom right pixel to test geotransform */
1125 0 : const double br_x = CPLStrtod(
1126 : CPLGetXMLValue(psPos, "lowerRightCorner.mapCoordinate.easting",
1127 : "0.0"),
1128 : nullptr);
1129 0 : const double br_y = CPLStrtod(
1130 : CPLGetXMLValue(psPos, "lowerRightCorner.mapCoordinate.northing",
1131 : "0.0"),
1132 : nullptr);
1133 0 : const double testx =
1134 0 : poDS->adfGeoTransform[0] +
1135 0 : poDS->adfGeoTransform[1] * (poDS->nRasterXSize - 0.5) +
1136 0 : poDS->adfGeoTransform[2] * (poDS->nRasterYSize - 0.5);
1137 0 : const double testy =
1138 0 : poDS->adfGeoTransform[3] +
1139 0 : poDS->adfGeoTransform[4] * (poDS->nRasterXSize - 0.5) +
1140 0 : poDS->adfGeoTransform[5] * (poDS->nRasterYSize - 0.5);
1141 :
1142 : /* Give 1/4 pixel numerical error leeway in calculating location
1143 : based on affine transform */
1144 0 : if ((fabs(testx - br_x) >
1145 0 : fabs(0.25 *
1146 0 : (poDS->adfGeoTransform[1] + poDS->adfGeoTransform[2]))) ||
1147 0 : (fabs(testy - br_y) > fabs(0.25 * (poDS->adfGeoTransform[4] +
1148 0 : poDS->adfGeoTransform[5]))))
1149 : {
1150 0 : CPLError(CE_Warning, CPLE_AppDefined,
1151 : "Unexpected error in calculating affine transform: "
1152 : "corner coordinates inconsistent.");
1153 : }
1154 : else
1155 : {
1156 0 : poDS->bHaveGeoTransform = TRUE;
1157 : }
1158 : }
1159 : }
1160 :
1161 : /* -------------------------------------------------------------------- */
1162 : /* Collect Projection String Information */
1163 : /* -------------------------------------------------------------------- */
1164 : CPLXMLNode *psEllipsoid =
1165 5 : CPLGetXMLNode(psImageAttributes,
1166 : "geographicInformation.referenceEllipsoidParameters");
1167 :
1168 5 : if (psEllipsoid != nullptr)
1169 : {
1170 10 : OGRSpatialReference oLL, oPrj;
1171 5 : oLL.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1172 5 : oPrj.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1173 :
1174 : const char *pszEllipsoidName =
1175 5 : CPLGetXMLValue(psEllipsoid, "ellipsoidName", "");
1176 : double minor_axis =
1177 5 : CPLAtof(CPLGetXMLValue(psEllipsoid, "semiMinorAxis", "0.0"));
1178 : double major_axis =
1179 5 : CPLAtof(CPLGetXMLValue(psEllipsoid, "semiMajorAxis", "0.0"));
1180 :
1181 5 : if (EQUAL(pszEllipsoidName, "") || (minor_axis == 0.0) ||
1182 : (major_axis == 0.0))
1183 : {
1184 0 : CPLError(CE_Warning, CPLE_AppDefined,
1185 : "Warning- incomplete"
1186 : " ellipsoid information. Using wgs-84 parameters.\n");
1187 0 : oLL.SetWellKnownGeogCS("WGS84");
1188 0 : oPrj.SetWellKnownGeogCS("WGS84");
1189 : }
1190 5 : else if (EQUAL(pszEllipsoidName, "WGS84"))
1191 : {
1192 5 : oLL.SetWellKnownGeogCS("WGS84");
1193 5 : oPrj.SetWellKnownGeogCS("WGS84");
1194 : }
1195 : else
1196 : {
1197 0 : const double inv_flattening =
1198 0 : major_axis / (major_axis - minor_axis);
1199 0 : oLL.SetGeogCS("", "", pszEllipsoidName, major_axis, inv_flattening);
1200 0 : oPrj.SetGeogCS("", "", pszEllipsoidName, major_axis,
1201 : inv_flattening);
1202 : }
1203 :
1204 5 : if (psMapProjection != nullptr)
1205 : {
1206 : const char *pszProj =
1207 0 : CPLGetXMLValue(psMapProjection, "mapProjectionDescriptor", "");
1208 0 : bool bUseProjInfo = false;
1209 :
1210 : CPLXMLNode *psUtmParams =
1211 0 : CPLGetXMLNode(psMapProjection, "utmProjectionParameters");
1212 :
1213 : CPLXMLNode *psNspParams =
1214 0 : CPLGetXMLNode(psMapProjection, "nspProjectionParameters");
1215 :
1216 0 : if ((psUtmParams != nullptr) && poDS->bHaveGeoTransform)
1217 : {
1218 : /* double origEasting, origNorthing; */
1219 0 : bool bNorth = true;
1220 :
1221 : const int utmZone =
1222 0 : atoi(CPLGetXMLValue(psUtmParams, "utmZone", ""));
1223 : const char *pszHemisphere =
1224 0 : CPLGetXMLValue(psUtmParams, "hemisphere", "");
1225 : #if 0
1226 : origEasting = CPLStrtod(CPLGetXMLValue(
1227 : psUtmParams, "mapOriginFalseEasting", "0.0" ), nullptr);
1228 : origNorthing = CPLStrtod(CPLGetXMLValue(
1229 : psUtmParams, "mapOriginFalseNorthing", "0.0" ), nullptr);
1230 : #endif
1231 0 : if (STARTS_WITH_CI(pszHemisphere, "southern"))
1232 0 : bNorth = FALSE;
1233 :
1234 0 : if (STARTS_WITH_CI(pszProj, "UTM"))
1235 : {
1236 0 : oPrj.SetUTM(utmZone, bNorth);
1237 0 : bUseProjInfo = true;
1238 0 : }
1239 : }
1240 0 : else if ((psNspParams != nullptr) && poDS->bHaveGeoTransform)
1241 : {
1242 0 : const double origEasting = CPLStrtod(
1243 : CPLGetXMLValue(psNspParams, "mapOriginFalseEasting", "0.0"),
1244 : nullptr);
1245 : const double origNorthing =
1246 0 : CPLStrtod(CPLGetXMLValue(psNspParams,
1247 : "mapOriginFalseNorthing", "0.0"),
1248 : nullptr);
1249 0 : const double copLong = CPLStrtod(
1250 : CPLGetXMLValue(psNspParams, "centerOfProjectionLongitude",
1251 : "0.0"),
1252 : nullptr);
1253 0 : const double copLat = CPLStrtod(
1254 : CPLGetXMLValue(psNspParams, "centerOfProjectionLatitude",
1255 : "0.0"),
1256 : nullptr);
1257 0 : const double sP1 = CPLStrtod(
1258 : CPLGetXMLValue(psNspParams, "standardParallels1", "0.0"),
1259 : nullptr);
1260 0 : const double sP2 = CPLStrtod(
1261 : CPLGetXMLValue(psNspParams, "standardParallels2", "0.0"),
1262 : nullptr);
1263 :
1264 0 : if (STARTS_WITH_CI(pszProj, "ARC"))
1265 : {
1266 : /* Albers Conical Equal Area */
1267 0 : oPrj.SetACEA(sP1, sP2, copLat, copLong, origEasting,
1268 : origNorthing);
1269 0 : bUseProjInfo = true;
1270 : }
1271 0 : else if (STARTS_WITH_CI(pszProj, "LCC"))
1272 : {
1273 : /* Lambert Conformal Conic */
1274 0 : oPrj.SetLCC(sP1, sP2, copLat, copLong, origEasting,
1275 : origNorthing);
1276 0 : bUseProjInfo = true;
1277 : }
1278 0 : else if (STARTS_WITH_CI(pszProj, "STPL"))
1279 : {
1280 : /* State Plate
1281 : ASSUMPTIONS: "zone" in product.xml matches USGS
1282 : definition as expected by ogr spatial reference; NAD83
1283 : zones (versus NAD27) are assumed. */
1284 :
1285 : const int nSPZone =
1286 0 : atoi(CPLGetXMLValue(psNspParams, "zone", "1"));
1287 :
1288 0 : oPrj.SetStatePlane(nSPZone, TRUE, nullptr, 0.0);
1289 0 : bUseProjInfo = true;
1290 : }
1291 : }
1292 :
1293 0 : if (bUseProjInfo)
1294 : {
1295 0 : poDS->m_oSRS = std::move(oPrj);
1296 : }
1297 : else
1298 : {
1299 0 : CPLError(CE_Warning, CPLE_AppDefined,
1300 : "Unable to interpret "
1301 : "projection information; check mapProjection info in "
1302 : "product.xml!");
1303 : }
1304 : }
1305 :
1306 5 : poDS->m_oGCPSRS = std::move(oLL);
1307 : }
1308 :
1309 : /* -------------------------------------------------------------------- */
1310 : /* Collect GCPs. */
1311 : /* -------------------------------------------------------------------- */
1312 5 : CPLXMLNode *psGeoGrid = CPLGetXMLNode(
1313 : psImageAttributes, "geographicInformation.geolocationGrid");
1314 :
1315 5 : if (psGeoGrid != nullptr)
1316 : {
1317 : /* count GCPs */
1318 5 : poDS->nGCPCount = 0;
1319 :
1320 25 : for (psNode = psGeoGrid->psChild; psNode != nullptr;
1321 20 : psNode = psNode->psNext)
1322 : {
1323 20 : if (EQUAL(psNode->pszValue, "imageTiePoint"))
1324 20 : poDS->nGCPCount++;
1325 : }
1326 :
1327 5 : poDS->pasGCPList = reinterpret_cast<GDAL_GCP *>(
1328 5 : CPLCalloc(sizeof(GDAL_GCP), poDS->nGCPCount));
1329 :
1330 5 : poDS->nGCPCount = 0;
1331 :
1332 25 : for (psNode = psGeoGrid->psChild; psNode != nullptr;
1333 20 : psNode = psNode->psNext)
1334 : {
1335 20 : GDAL_GCP *psGCP = poDS->pasGCPList + poDS->nGCPCount;
1336 :
1337 20 : if (!EQUAL(psNode->pszValue, "imageTiePoint"))
1338 0 : continue;
1339 :
1340 20 : poDS->nGCPCount++;
1341 :
1342 : char szID[32];
1343 20 : snprintf(szID, sizeof(szID), "%d", poDS->nGCPCount);
1344 20 : psGCP->pszId = CPLStrdup(szID);
1345 20 : psGCP->pszInfo = CPLStrdup("");
1346 20 : psGCP->dfGCPPixel =
1347 20 : CPLAtof(CPLGetXMLValue(psNode, "imageCoordinate.pixel", "0")) +
1348 : 0.5;
1349 20 : psGCP->dfGCPLine =
1350 20 : CPLAtof(CPLGetXMLValue(psNode, "imageCoordinate.line", "0")) +
1351 : 0.5;
1352 20 : psGCP->dfGCPX = CPLAtof(
1353 : CPLGetXMLValue(psNode, "geodeticCoordinate.longitude", ""));
1354 20 : psGCP->dfGCPY = CPLAtof(
1355 : CPLGetXMLValue(psNode, "geodeticCoordinate.latitude", ""));
1356 20 : psGCP->dfGCPZ = CPLAtof(
1357 : CPLGetXMLValue(psNode, "geodeticCoordinate.height", ""));
1358 : }
1359 : }
1360 :
1361 5 : CPLFree(pszPath);
1362 :
1363 : /* -------------------------------------------------------------------- */
1364 : /* Collect RPC. */
1365 : /* -------------------------------------------------------------------- */
1366 5 : CPLXMLNode *psRationalFunctions = CPLGetXMLNode(
1367 : psImageAttributes, "geographicInformation.rationalFunctions");
1368 5 : if (psRationalFunctions != nullptr)
1369 : {
1370 5 : char **papszRPC = nullptr;
1371 : static const char *const apszXMLToGDALMapping[] = {
1372 : "biasError",
1373 : "ERR_BIAS",
1374 : "randomError",
1375 : "ERR_RAND",
1376 : //"lineFitQuality", "????",
1377 : //"pixelFitQuality", "????",
1378 : "lineOffset",
1379 : "LINE_OFF",
1380 : "pixelOffset",
1381 : "SAMP_OFF",
1382 : "latitudeOffset",
1383 : "LAT_OFF",
1384 : "longitudeOffset",
1385 : "LONG_OFF",
1386 : "heightOffset",
1387 : "HEIGHT_OFF",
1388 : "lineScale",
1389 : "LINE_SCALE",
1390 : "pixelScale",
1391 : "SAMP_SCALE",
1392 : "latitudeScale",
1393 : "LAT_SCALE",
1394 : "longitudeScale",
1395 : "LONG_SCALE",
1396 : "heightScale",
1397 : "HEIGHT_SCALE",
1398 : "lineNumeratorCoefficients",
1399 : "LINE_NUM_COEFF",
1400 : "lineDenominatorCoefficients",
1401 : "LINE_DEN_COEFF",
1402 : "pixelNumeratorCoefficients",
1403 : "SAMP_NUM_COEFF",
1404 : "pixelDenominatorCoefficients",
1405 : "SAMP_DEN_COEFF",
1406 : };
1407 85 : for (size_t i = 0; i < CPL_ARRAYSIZE(apszXMLToGDALMapping); i += 2)
1408 : {
1409 160 : const char *pszValue = CPLGetXMLValue(
1410 80 : psRationalFunctions, apszXMLToGDALMapping[i], nullptr);
1411 80 : if (pszValue)
1412 80 : papszRPC = CSLSetNameValue(
1413 80 : papszRPC, apszXMLToGDALMapping[i + 1], pszValue);
1414 : }
1415 5 : poDS->GDALDataset::SetMetadata(papszRPC, "RPC");
1416 5 : CSLDestroy(papszRPC);
1417 : }
1418 :
1419 : /* -------------------------------------------------------------------- */
1420 : /* Initialize any PAM information. */
1421 : /* -------------------------------------------------------------------- */
1422 5 : CPLString osDescription;
1423 :
1424 5 : switch (eCalib)
1425 : {
1426 0 : case Sigma0:
1427 : osDescription.Printf("RADARSAT_2_CALIB:SIGMA0:%s",
1428 0 : osMDFilename.c_str());
1429 0 : break;
1430 1 : case Beta0:
1431 : osDescription.Printf("RADARSAT_2_CALIB:BETA0:%s",
1432 1 : osMDFilename.c_str());
1433 1 : break;
1434 0 : case Gamma:
1435 : osDescription.Printf("RADARSAT_2_CALIB:GAMMA0:%s",
1436 0 : osMDFilename.c_str());
1437 0 : break;
1438 0 : case Uncalib:
1439 : osDescription.Printf("RADARSAT_2_CALIB:UNCALIB:%s",
1440 0 : osMDFilename.c_str());
1441 0 : break;
1442 4 : default:
1443 4 : osDescription = osMDFilename;
1444 : }
1445 :
1446 5 : if (eCalib != None)
1447 1 : poDS->papszExtraFiles =
1448 1 : CSLAddString(poDS->papszExtraFiles, osMDFilename);
1449 :
1450 : /* -------------------------------------------------------------------- */
1451 : /* Initialize any PAM information. */
1452 : /* -------------------------------------------------------------------- */
1453 5 : poDS->SetDescription(osDescription);
1454 :
1455 5 : poDS->SetPhysicalFilename(osMDFilename);
1456 5 : poDS->SetSubdatasetName(osDescription);
1457 :
1458 5 : poDS->TryLoadXML();
1459 :
1460 : /* -------------------------------------------------------------------- */
1461 : /* Check for overviews. */
1462 : /* -------------------------------------------------------------------- */
1463 5 : poDS->oOvManager.Initialize(poDS, ":::VIRTUAL:::");
1464 :
1465 5 : return poDS;
1466 : }
1467 :
1468 : /************************************************************************/
1469 : /* GetGCPCount() */
1470 : /************************************************************************/
1471 :
1472 0 : int RS2Dataset::GetGCPCount()
1473 :
1474 : {
1475 0 : return nGCPCount;
1476 : }
1477 :
1478 : /************************************************************************/
1479 : /* GetGCPSpatialRef() */
1480 : /************************************************************************/
1481 :
1482 0 : const OGRSpatialReference *RS2Dataset::GetGCPSpatialRef() const
1483 :
1484 : {
1485 0 : return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
1486 : }
1487 :
1488 : /************************************************************************/
1489 : /* GetGCPs() */
1490 : /************************************************************************/
1491 :
1492 0 : const GDAL_GCP *RS2Dataset::GetGCPs()
1493 :
1494 : {
1495 0 : return pasGCPList;
1496 : }
1497 :
1498 : /************************************************************************/
1499 : /* GetSpatialRef() */
1500 : /************************************************************************/
1501 :
1502 0 : const OGRSpatialReference *RS2Dataset::GetSpatialRef() const
1503 :
1504 : {
1505 0 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
1506 : }
1507 :
1508 : /************************************************************************/
1509 : /* GetGeoTransform() */
1510 : /************************************************************************/
1511 :
1512 0 : CPLErr RS2Dataset::GetGeoTransform(double *padfTransform)
1513 :
1514 : {
1515 0 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
1516 :
1517 0 : if (bHaveGeoTransform)
1518 0 : return CE_None;
1519 :
1520 0 : return CE_Failure;
1521 : }
1522 :
1523 : /************************************************************************/
1524 : /* GetMetadataDomainList() */
1525 : /************************************************************************/
1526 :
1527 0 : char **RS2Dataset::GetMetadataDomainList()
1528 : {
1529 0 : return BuildMetadataDomainList(GDALDataset::GetMetadataDomainList(), TRUE,
1530 0 : "SUBDATASETS", nullptr);
1531 : }
1532 :
1533 : /************************************************************************/
1534 : /* GetMetadata() */
1535 : /************************************************************************/
1536 :
1537 1 : char **RS2Dataset::GetMetadata(const char *pszDomain)
1538 :
1539 : {
1540 1 : if (pszDomain != nullptr && STARTS_WITH_CI(pszDomain, "SUBDATASETS") &&
1541 0 : papszSubDatasets != nullptr)
1542 0 : return papszSubDatasets;
1543 :
1544 1 : return GDALDataset::GetMetadata(pszDomain);
1545 : }
1546 :
1547 : /************************************************************************/
1548 : /* GDALRegister_RS2() */
1549 : /************************************************************************/
1550 :
1551 1512 : void GDALRegister_RS2()
1552 :
1553 : {
1554 1512 : if (GDALGetDriverByName("RS2") != nullptr)
1555 295 : return;
1556 :
1557 1217 : GDALDriver *poDriver = new GDALDriver();
1558 :
1559 1217 : poDriver->SetDescription("RS2");
1560 1217 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1561 1217 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "RadarSat 2 XML Product");
1562 1217 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/rs2.html");
1563 1217 : poDriver->SetMetadataItem(GDAL_DMD_SUBDATASETS, "YES");
1564 1217 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1565 :
1566 1217 : poDriver->pfnOpen = RS2Dataset::Open;
1567 1217 : poDriver->pfnIdentify = RS2Dataset::Identify;
1568 :
1569 1217 : GetGDALDriverManager()->RegisterDriver(poDriver);
1570 : }
|
