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