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