1 : /******************************************************************************
2 : * $Id: rs2dataset.cpp 23570 2011-12-13 22:23:32Z rouault $
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 23570 2011-12-13 22:23:32Z rouault $");
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 12 : 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 12 : bool retVal = false;
56 :
57 12 : pszLutFile = VSIStrdup(CPLFormFilename(pszPath, pszLut, NULL));
58 :
59 12 : psLut = CPLParseXMLFile(pszLutFile);
60 :
61 12 : if (psLut != NULL)
62 : {
63 12 : CPLDestroyXMLNode(psLut);
64 12 : retVal = true;
65 : }
66 :
67 12 : CPLFree(pszLutFile);
68 :
69 12 : 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 4 : RS2RasterBand::RS2RasterBand( RS2Dataset *poDSIn,
144 : GDALDataType eDataTypeIn,
145 : const char *pszPole,
146 4 : GDALDataset *poBandFileIn )
147 :
148 : {
149 : GDALRasterBand *poSrcBand;
150 :
151 4 : poDS = poDSIn;
152 4 : poBandFile = poBandFileIn;
153 :
154 4 : poSrcBand = poBandFile->GetRasterBand( 1 );
155 :
156 4 : poSrcBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
157 :
158 4 : eDataType = eDataTypeIn;
159 :
160 4 : if( *pszPole != '\0' )
161 4 : SetMetadataItem( "POLARIMETRIC_INTERP", pszPole );
162 4 : }
163 :
164 : /************************************************************************/
165 : /* RSRasterBand() */
166 : /************************************************************************/
167 :
168 4 : RS2RasterBand::~RS2RasterBand()
169 :
170 : {
171 4 : if( poBandFile != NULL )
172 4 : GDALClose( (GDALRasterBandH) poBandFile );
173 4 : }
174 :
175 : /************************************************************************/
176 : /* IReadBlock() */
177 : /************************************************************************/
178 :
179 40 : 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 40 : 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 40 : 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 40 : 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 40 : nRequestXSize = nBlockXSize;
215 : }
216 40 : 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 40 : 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 40 : 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 40 : 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 40 : 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 4 : void RS2CalibRasterBand::ReadLUT() {
318 : CPLXMLNode *psLUT;
319 : char **papszLUTList;
320 :
321 4 : psLUT = CPLParseXMLFile(m_pszLUTFile);
322 :
323 : this->m_nfOffset = (float) CPLAtof(CPLGetXMLValue(psLUT, "=lut.offset",
324 4 : "0.0"));
325 :
326 : papszLUTList = CSLTokenizeString2( CPLGetXMLValue(psLUT,
327 4 : "=lut.gains", ""), " ", CSLT_HONOURSTRINGS);
328 :
329 4 : this->m_nTableSize = CSLCount(papszLUTList);
330 :
331 4 : this->m_nfTable = (float *)CPLMalloc(sizeof(float) * this->m_nTableSize);
332 :
333 1028 : for (int i = 0; i < this->m_nTableSize; i++) {
334 1024 : m_nfTable[i] = (float) CPLAtof(papszLUTList[i]);
335 : }
336 :
337 4 : CPLDestroyXMLNode(psLUT);
338 :
339 4 : CSLDestroy(papszLUTList);
340 4 : }
341 :
342 : /************************************************************************/
343 : /* RS2CalibRasterBand() */
344 : /************************************************************************/
345 :
346 4 : RS2CalibRasterBand::RS2CalibRasterBand( RS2Dataset *poDataset,
347 : const char *pszPolarization, GDALDataType eType, GDALDataset *poBandDataset,
348 4 : eCalibration eCalib, const char *pszLUT )
349 : {
350 4 : this->poDS = poDataset;
351 :
352 4 : if (*pszPolarization != '\0') {
353 4 : this->SetMetadataItem( "POLARIMETRIC_INTERP", pszPolarization );
354 : }
355 :
356 4 : this->m_eType = eType;
357 4 : this->m_poBandDataset = poBandDataset;
358 4 : this->m_eCalib = eCalib;
359 4 : this->m_pszLUTFile = VSIStrdup(pszLUT);
360 :
361 4 : this->m_nfTable = NULL;
362 4 : this->m_nTableSize = 0;
363 :
364 4 : if (eType == GDT_CInt16)
365 0 : this->eDataType = GDT_CFloat32;
366 : else
367 4 : this->eDataType = GDT_Float32;
368 :
369 4 : GDALRasterBand *poRasterBand = poBandDataset->GetRasterBand( 1 );
370 4 : poRasterBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
371 :
372 :
373 4 : this->ReadLUT();
374 4 : }
375 :
376 : /************************************************************************/
377 : /* ~RS2CalibRasterBand() */
378 : /************************************************************************/
379 :
380 4 : RS2CalibRasterBand::~RS2CalibRasterBand() {
381 4 : if (this->m_nfTable != NULL)
382 4 : CPLFree(this->m_nfTable);
383 :
384 4 : if (this->m_pszLUTFile != NULL)
385 4 : CPLFree(this->m_pszLUTFile);
386 :
387 4 : if (this->m_poBandDataset != NULL)
388 4 : GDALClose( this->m_poBandDataset );
389 4 : }
390 :
391 : /************************************************************************/
392 : /* IReadBlock() */
393 : /************************************************************************/
394 :
395 40 : 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 40 : 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 40 : nRequestYSize = nBlockYSize;
415 : }
416 :
417 40 : 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 40 : 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 40 : else if (this->m_eType == GDT_Byte) {
491 : GByte *pnImageTmp;
492 : pnImageTmp = (GByte *)CPLMalloc(nBlockXSize * nBlockYSize *
493 40 : 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 40 : 1, NULL, 1, 1, 0);
501 :
502 : /* iterate over detected values */
503 80 : for (int i = 0; i < nBlockYSize; i++) {
504 840 : for (int j = 0; j < nBlockXSize; j++) {
505 800 : int nPixOff = (i * nBlockXSize) + j;
506 :
507 800 : ((float *)pImage)[nPixOff] = ((pnImageTmp[nPixOff] *
508 800 : pnImageTmp[nPixOff]) +
509 1600 : this->m_nfOffset)/m_nfTable[nBlockXOff + j];
510 : }
511 : }
512 40 : CPLFree(pnImageTmp);
513 : }
514 : else {
515 0 : CPLAssert( FALSE );
516 0 : return CE_Failure;
517 : }
518 40 : return eErr;
519 : }
520 :
521 :
522 : /************************************************************************/
523 : /* ==================================================================== */
524 : /* RS2Dataset */
525 : /* ==================================================================== */
526 : /************************************************************************/
527 :
528 : /************************************************************************/
529 : /* RS2Dataset() */
530 : /************************************************************************/
531 :
532 4 : RS2Dataset::RS2Dataset()
533 : {
534 4 : psProduct = NULL;
535 :
536 4 : nGCPCount = 0;
537 4 : pasGCPList = NULL;
538 4 : pszGCPProjection = CPLStrdup("");
539 4 : pszProjection = CPLStrdup("");
540 4 : adfGeoTransform[0] = 0.0;
541 4 : adfGeoTransform[1] = 1.0;
542 4 : adfGeoTransform[2] = 0.0;
543 4 : adfGeoTransform[3] = 0.0;
544 4 : adfGeoTransform[4] = 0.0;
545 4 : adfGeoTransform[5] = 1.0;
546 4 : bHaveGeoTransform = FALSE;
547 :
548 4 : papszSubDatasets = NULL;
549 4 : papszExtraFiles = NULL;
550 4 : }
551 :
552 : /************************************************************************/
553 : /* ~RS2Dataset() */
554 : /************************************************************************/
555 :
556 4 : RS2Dataset::~RS2Dataset()
557 :
558 : {
559 4 : FlushCache();
560 :
561 4 : CPLDestroyXMLNode( psProduct );
562 4 : CPLFree( pszProjection );
563 :
564 4 : CPLFree( pszGCPProjection );
565 4 : if( nGCPCount > 0 )
566 : {
567 4 : GDALDeinitGCPs( nGCPCount, pasGCPList );
568 4 : CPLFree( pasGCPList );
569 : }
570 :
571 4 : CloseDependentDatasets();
572 :
573 4 : CSLDestroy( papszSubDatasets );
574 4 : CSLDestroy( papszExtraFiles );
575 4 : }
576 :
577 : /************************************************************************/
578 : /* CloseDependentDatasets() */
579 : /************************************************************************/
580 :
581 4 : int RS2Dataset::CloseDependentDatasets()
582 : {
583 4 : int bHasDroppedRef = GDALPamDataset::CloseDependentDatasets();
584 :
585 4 : if (nBands != 0)
586 4 : bHasDroppedRef = TRUE;
587 :
588 12 : for( int iBand = 0; iBand < nBands; iBand++ )
589 : {
590 8 : delete papoBands[iBand];
591 : }
592 4 : nBands = 0;
593 :
594 4 : 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 25704 : int RS2Dataset::Identify( GDALOpenInfo *poOpenInfo )
616 : {
617 :
618 : /* Check for the case where we're trying to read the calibrated data: */
619 25704 : if (EQUALN("RADARSAT_2_CALIB:",poOpenInfo->pszFilename,17)) {
620 2 : return 1;
621 : }
622 :
623 : /* Check for directory access when there is a product.xml file in the
624 : directory. */
625 25702 : if( poOpenInfo->bIsDirectory )
626 : {
627 : VSIStatBufL sStat;
628 :
629 : CPLString osMDFilename =
630 84 : CPLFormCIFilename( poOpenInfo->pszFilename, "product.xml", NULL );
631 :
632 84 : if( VSIStatL( osMDFilename, &sStat ) == 0 )
633 0 : return TRUE;
634 : else
635 84 : return FALSE;
636 : }
637 :
638 : /* otherwise, do our normal stuff */
639 25618 : if( strlen(poOpenInfo->pszFilename) < 11
640 : || !EQUAL(poOpenInfo->pszFilename + strlen(poOpenInfo->pszFilename)-11,
641 : "product.xml") )
642 25616 : return 0;
643 :
644 2 : if( poOpenInfo->nHeaderBytes < 100 )
645 0 : return 0;
646 :
647 2 : if( strstr((const char *) poOpenInfo->pabyHeader, "/rs2" ) == NULL
648 : || strstr((const char *) poOpenInfo->pabyHeader, "<product" ) == NULL)
649 0 : return 0;
650 :
651 2 : return 1;
652 :
653 : }
654 :
655 :
656 : /************************************************************************/
657 : /* Open() */
658 : /************************************************************************/
659 :
660 6792 : GDALDataset *RS2Dataset::Open( GDALOpenInfo * poOpenInfo )
661 :
662 : {
663 : /* -------------------------------------------------------------------- */
664 : /* Is this a RADARSAT-2 Product.xml definition? */
665 : /* -------------------------------------------------------------------- */
666 6792 : if ( !RS2Dataset::Identify( poOpenInfo ) ) {
667 6788 : return NULL;
668 : }
669 :
670 : /* -------------------------------------------------------------------- */
671 : /* Get subdataset information, if relevant */
672 : /* -------------------------------------------------------------------- */
673 4 : CPLString osMDFilename;
674 4 : const char *pszFilename = poOpenInfo->pszFilename;
675 4 : eCalibration eCalib = None;
676 :
677 4 : if (EQUALN("RADARSAT_2_CALIB:",pszFilename,17)) {
678 2 : pszFilename += 17;
679 :
680 2 : if (EQUALN("BETA0",pszFilename,5))
681 2 : 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 14 : while ( *pszFilename != '\0' && *pszFilename != ':' )
693 10 : pszFilename++;
694 :
695 2 : if (*pszFilename == ':')
696 2 : 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 2 : if( VSIStatL( pszFilename, &sStat ) == 0 )
702 2 : poOpenInfo->bIsDirectory = VSI_ISDIR( sStat.st_mode );
703 : }
704 :
705 4 : if( poOpenInfo->bIsDirectory )
706 : {
707 : osMDFilename =
708 0 : CPLFormCIFilename( pszFilename, "product.xml", NULL );
709 : }
710 : else
711 4 : osMDFilename = pszFilename;
712 :
713 : /* -------------------------------------------------------------------- */
714 : /* Ingest the Product.xml file. */
715 : /* -------------------------------------------------------------------- */
716 : CPLXMLNode *psProduct, *psImageAttributes, *psImageGenerationParameters;
717 :
718 4 : psProduct = CPLParseXMLFile( osMDFilename );
719 4 : if( psProduct == NULL )
720 0 : return NULL;
721 :
722 : /* -------------------------------------------------------------------- */
723 : /* Confirm the requested access is supported. */
724 : /* -------------------------------------------------------------------- */
725 4 : 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 4 : psImageAttributes = CPLGetXMLNode(psProduct, "=product.imageAttributes" );
735 4 : 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 4 : "=product.imageGenerationParameters" );
745 4 : 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 4 : RS2Dataset *poDS = new RS2Dataset();
756 :
757 4 : poDS->psProduct = psProduct;
758 :
759 : /* -------------------------------------------------------------------- */
760 : /* Get overall image information. */
761 : /* -------------------------------------------------------------------- */
762 : poDS->nRasterXSize =
763 : atoi(CPLGetXMLValue( psImageAttributes,
764 : "rasterAttributes.numberOfSamplesPerLine",
765 8 : "-1" ));
766 : poDS->nRasterYSize =
767 : atoi(CPLGetXMLValue( psImageAttributes,
768 : "rasterAttributes.numberofLines",
769 4 : "-1" ));
770 4 : 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 4 : char *pszBeta0LUT = NULL;
785 4 : char *pszGammaLUT = NULL;
786 4 : char *pszSigma0LUT = NULL;
787 4 : int bCanCalib = 0;
788 :
789 : const char *pszProductType = CPLGetXMLValue( psImageGenerationParameters,
790 : "generalProcessingInformation.productType",
791 4 : "UNK" );
792 :
793 4 : 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 4 : if (!(EQUALN(pszProductType, "UNK", 3) ||
799 : EQUALN(pszProductType, "SSG", 3) ||
800 : EQUALN(pszProductType, "SPG", 3)))
801 : {
802 4 : 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 4 : "" );
814 : int nBitsPerSample =
815 : atoi( CPLGetXMLValue( psImageAttributes,
816 4 : "rasterAttributes.bitsPerSample", "" ) );
817 :
818 4 : if( nBitsPerSample == 16 && EQUAL(pszDataType,"Complex") )
819 0 : eDataType = GDT_CInt16;
820 4 : else if( nBitsPerSample == 16 && EQUALN(pszDataType,"Mag",3) )
821 0 : eDataType = GDT_UInt16;
822 8 : else if( nBitsPerSample == 8 && EQUALN(pszDataType,"Mag",3) )
823 4 : 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 4 : "rasterAttributes.sampledPixelSpacing", "UNK" );
836 4 : poDS->SetMetadataItem( "PIXEL_SPACING", pszPixelSpacing );
837 :
838 : const char *pszLineSpacing = CPLGetXMLValue( psImageAttributes,
839 4 : "rasterAttributes.sampledLineSpacing", "UNK" );
840 4 : poDS->SetMetadataItem( "LINE_SPACING", pszLineSpacing );
841 :
842 : /* -------------------------------------------------------------------- */
843 : /* Open each of the data files as a complex band. */
844 : /* -------------------------------------------------------------------- */
845 : CPLXMLNode *psNode;
846 4 : char *pszPath = CPLStrdup(CPLGetPath( osMDFilename ));
847 : char *pszBuf;
848 4 : int nFLen = strlen(osMDFilename);
849 :
850 32 : for( psNode = psImageAttributes->psChild;
851 : psNode != NULL;
852 : psNode = psNode->psNext )
853 : {
854 : const char *pszBasename;
855 28 : if( psNode->eType != CXT_Element
856 : || !(EQUAL(psNode->pszValue,"fullResolutionImageData")
857 : || EQUAL(psNode->pszValue,"lookupTable")) )
858 8 : continue;
859 :
860 20 : if ( EQUAL(psNode->pszValue, "lookupTable") && bCanCalib ) {
861 : /* Determine which incidence angle correction this is */
862 : const char *pszLUTType = CPLGetXMLValue( psNode,
863 12 : "incidenceAngleCorrection", "" );
864 12 : const char *pszLUTFile = CPLGetXMLValue( psNode, "", "" );
865 : CPLString osLUTFilePath = CPLFormFilename( pszPath, pszLUTFile,
866 12 : NULL );
867 :
868 12 : if (EQUAL(pszLUTType, ""))
869 0 : continue;
870 12 : else if (EQUAL(pszLUTType, "Beta Nought") &&
871 : IsValidXMLFile(pszPath,pszLUTFile))
872 : {
873 : poDS->papszExtraFiles =
874 4 : CSLAddString( poDS->papszExtraFiles, osLUTFilePath );
875 :
876 4 : pszBuf = (char *)CPLMalloc(nFLen + 27);
877 4 : pszBeta0LUT = VSIStrdup( pszLUTFile );
878 4 : poDS->SetMetadataItem( "BETA_NOUGHT_LUT", pszLUTFile );
879 :
880 : sprintf(pszBuf, "RADARSAT_2_CALIB:BETA0:%s",
881 4 : osMDFilename.c_str() );
882 : poDS->papszSubDatasets = CSLSetNameValue(
883 4 : poDS->papszSubDatasets, "SUBDATASET_3_NAME", pszBuf );
884 : poDS->papszSubDatasets = CSLSetNameValue(
885 : poDS->papszSubDatasets, "SUBDATASET_3_DESC",
886 4 : "Beta Nought calibrated" );
887 4 : CPLFree(pszBuf);
888 : }
889 8 : else if (EQUAL(pszLUTType, "Sigma Nought") &&
890 : IsValidXMLFile(pszPath,pszLUTFile))
891 : {
892 : poDS->papszExtraFiles =
893 4 : CSLAddString( poDS->papszExtraFiles, osLUTFilePath );
894 :
895 4 : pszBuf = (char *)CPLMalloc(nFLen + 27);
896 4 : pszSigma0LUT = VSIStrdup( pszLUTFile );
897 4 : poDS->SetMetadataItem( "SIGMA_NOUGHT_LUT", pszLUTFile );
898 :
899 : sprintf(pszBuf, "RADARSAT_2_CALIB:SIGMA0:%s",
900 4 : osMDFilename.c_str() );
901 : poDS->papszSubDatasets = CSLSetNameValue(
902 4 : poDS->papszSubDatasets, "SUBDATASET_2_NAME", pszBuf );
903 : poDS->papszSubDatasets = CSLSetNameValue(
904 : poDS->papszSubDatasets, "SUBDATASET_2_DESC",
905 4 : "Sigma Nought calibrated" );
906 4 : CPLFree(pszBuf);
907 : }
908 4 : else if (EQUAL(pszLUTType, "Gamma") &&
909 : IsValidXMLFile(pszPath,pszLUTFile))
910 : {
911 : poDS->papszExtraFiles =
912 4 : CSLAddString( poDS->papszExtraFiles, osLUTFilePath );
913 :
914 4 : pszBuf = (char *)CPLMalloc(nFLen + 27);
915 4 : pszGammaLUT = VSIStrdup( pszLUTFile );
916 4 : poDS->SetMetadataItem( "GAMMA_LUT", pszLUTFile );
917 : sprintf(pszBuf, "RADARSAT_2_CALIB:GAMMA:%s",
918 4 : osMDFilename.c_str());
919 : poDS->papszSubDatasets = CSLSetNameValue(
920 4 : poDS->papszSubDatasets, "SUBDATASET_4_NAME", pszBuf );
921 : poDS->papszSubDatasets = CSLSetNameValue(
922 : poDS->papszSubDatasets, "SUBDATASET_4_DESC",
923 4 : "Gamma calibrated" );
924 4 : CPLFree(pszBuf);
925 : }
926 12 : continue;
927 : }
928 :
929 : /* -------------------------------------------------------------------- */
930 : /* Fetch filename. */
931 : /* -------------------------------------------------------------------- */
932 8 : pszBasename = CPLGetXMLValue( psNode, "", "" );
933 8 : if( *pszBasename == '\0' )
934 0 : continue;
935 :
936 : /* -------------------------------------------------------------------- */
937 : /* Form full filename (path of product.xml + basename). */
938 : /* -------------------------------------------------------------------- */
939 : char *pszFullname =
940 8 : CPLStrdup(CPLFormFilename( pszPath, pszBasename, NULL ));
941 :
942 : /* -------------------------------------------------------------------- */
943 : /* Try and open the file. */
944 : /* -------------------------------------------------------------------- */
945 : GDALDataset *poBandFile;
946 :
947 8 : poBandFile = (GDALDataset *) GDALOpen( pszFullname, GA_ReadOnly );
948 8 : if( poBandFile == NULL )
949 : {
950 0 : CPLFree(pszFullname);
951 0 : continue;
952 : }
953 8 : 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 8 : pszFullname );
962 :
963 : /* -------------------------------------------------------------------- */
964 : /* Create the band. */
965 : /* -------------------------------------------------------------------- */
966 12 : if (eCalib == None || eCalib == Uncalib) {
967 : RS2RasterBand *poBand;
968 : poBand = new RS2RasterBand( poDS, eDataType,
969 : CPLGetXMLValue( psNode, "pole", "" ),
970 4 : poBandFile );
971 :
972 8 : poDS->SetBand( poDS->GetRasterCount() + 1, poBand );
973 : }
974 : else {
975 : const char *pszLUT;
976 4 : switch (eCalib) {
977 : case Sigma0:
978 0 : pszLUT = pszSigma0LUT;
979 0 : break;
980 : case Beta0:
981 4 : pszLUT = pszBeta0LUT;
982 4 : 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 4 : CPLFormFilename(pszPath, pszLUT, NULL) );
994 8 : poDS->SetBand( poDS->GetRasterCount() + 1, poBand );
995 : }
996 :
997 8 : CPLFree( pszFullname );
998 : }
999 :
1000 6 : if (poDS->papszSubDatasets != NULL && eCalib == None) {
1001 : char *pszBuf;
1002 2 : pszBuf = (char *)CPLMalloc(nFLen + 28);
1003 : sprintf(pszBuf, "RADARSAT_2_CALIB:UNCALIB:%s",
1004 2 : osMDFilename.c_str() );
1005 : poDS->papszSubDatasets = CSLSetNameValue( poDS->papszSubDatasets,
1006 2 : "SUBDATASET_1_NAME", pszBuf );
1007 : poDS->papszSubDatasets = CSLSetNameValue( poDS->papszSubDatasets,
1008 2 : "SUBDATASET_1_DESC", "Uncalibrated digital numbers" );
1009 2 : CPLFree(pszBuf);
1010 : }
1011 2 : else if (poDS->papszSubDatasets != NULL) {
1012 2 : CSLDestroy( poDS->papszSubDatasets );
1013 2 : poDS->papszSubDatasets = NULL;
1014 : }
1015 :
1016 : /* -------------------------------------------------------------------- */
1017 : /* Set the appropriate MATRIX_REPRESENTATION. */
1018 : /* -------------------------------------------------------------------- */
1019 :
1020 4 : 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 4 : "=product.sourceAttributes");
1032 : const char *pszItem;
1033 :
1034 : pszItem = CPLGetXMLValue( psSourceAttrs,
1035 4 : "satellite", "" );
1036 4 : poDS->SetMetadataItem( "SATELLITE_IDENTIFIER", pszItem );
1037 :
1038 : pszItem = CPLGetXMLValue( psSourceAttrs,
1039 4 : "sensor", "" );
1040 4 : poDS->SetMetadataItem( "SENSOR_IDENTIFIER", pszItem );
1041 :
1042 4 : if (psSourceAttrs != NULL) {
1043 : /* Get beam mode mnemonic */
1044 4 : pszItem = CPLGetXMLValue( psSourceAttrs, "beamModeMnemonic", "UNK" );
1045 4 : poDS->SetMetadataItem( "BEAM_MODE", pszItem );
1046 4 : pszItem = CPLGetXMLValue( psSourceAttrs, "rawDataStartTime", "UNK" );
1047 4 : poDS->SetMetadataItem( "ACQUISITION_START_TIME", pszItem );
1048 : }
1049 :
1050 : CPLXMLNode *psSarProcessingInformation =
1051 4 : CPLGetXMLNode( psProduct, "=product.imageGenerationParameters" );
1052 :
1053 4 : if (psSarProcessingInformation != NULL) {
1054 : /* Get incidence angle information */
1055 : pszItem = CPLGetXMLValue( psSarProcessingInformation,
1056 4 : "sarProcessingInformation.incidenceAngleNearRange", "UNK" );
1057 4 : poDS->SetMetadataItem( "NEAR_RANGE_INCIDENCE_ANGLE", pszItem );
1058 :
1059 : pszItem = CPLGetXMLValue( psSarProcessingInformation,
1060 4 : "sarProcessingInformation.incidenceAngleFarRange", "UNK" );
1061 4 : poDS->SetMetadataItem( "FAR_RANGE_INCIDENCE_ANGLE", pszItem );
1062 :
1063 : pszItem = CPLGetXMLValue( psSarProcessingInformation,
1064 4 : "sarProcessingInformation.slantRangeNearEdge", "UNK" );
1065 4 : poDS->SetMetadataItem( "SLANT_RANGE_NEAR_EDGE", pszItem );
1066 : }
1067 :
1068 : /*--------------------------------------------------------------------- */
1069 : /* Collect Map projection/Geotransform information, if present */
1070 : /* -------------------------------------------------------------------- */
1071 : CPLXMLNode *psMapProjection =
1072 : CPLGetXMLNode( psImageAttributes,
1073 4 : "geographicInformation.mapProjection" );
1074 :
1075 4 : if ( psMapProjection != NULL ) {
1076 : CPLXMLNode *psPos =
1077 0 : CPLGetXMLNode( psMapProjection, "positioningInformation" );
1078 :
1079 : pszItem = CPLGetXMLValue( psMapProjection,
1080 0 : "mapProjectionDescriptor", "UNK" );
1081 0 : poDS->SetMetadataItem( "MAP_PROJECTION_DESCRIPTOR", pszItem );
1082 :
1083 : pszItem = CPLGetXMLValue( psMapProjection,
1084 0 : "mapProjectionOrientation", "UNK" );
1085 0 : poDS->SetMetadataItem( "MAP_PROJECTION_ORIENTATION", pszItem );
1086 :
1087 : pszItem = CPLGetXMLValue( psMapProjection,
1088 0 : "resamplingKernel", "UNK" );
1089 0 : poDS->SetMetadataItem( "RESAMPLING_KERNEL", pszItem );
1090 :
1091 : pszItem = CPLGetXMLValue( psMapProjection,
1092 0 : "satelliteHeading", "UNK" );
1093 0 : poDS->SetMetadataItem( "SATELLITE_HEADING", pszItem );
1094 :
1095 0 : if (psPos != NULL) {
1096 : double testx, testy, br_x, br_y, tl_x, tl_y, tr_x, tr_y,
1097 : bl_x, bl_y;
1098 :
1099 : tl_x = strtod(CPLGetXMLValue( psPos,
1100 0 : "upperLeftCorner.mapCoordinate.easting", "0.0" ), NULL);
1101 : tl_y = strtod(CPLGetXMLValue( psPos,
1102 0 : "upperLeftCorner.mapCoordinate.northing", "0.0" ), NULL);
1103 : bl_x = strtod(CPLGetXMLValue( psPos,
1104 0 : "lowerLeftCorner.mapCoordinate.easting", "0.0" ), NULL);
1105 : bl_y = strtod(CPLGetXMLValue( psPos,
1106 0 : "lowerLeftCorner.mapCoordinate.northing", "0.0" ), NULL);
1107 : tr_x = strtod(CPLGetXMLValue( psPos,
1108 0 : "upperRightCorner.mapCoordinate.easting", "0.0" ), NULL);
1109 : tr_y = strtod(CPLGetXMLValue( psPos,
1110 0 : "upperRightCorner.mapCoordinate.northing", "0.0" ), NULL);
1111 0 : poDS->adfGeoTransform[1] = (tr_x - tl_x)/(poDS->nRasterXSize - 1);
1112 0 : poDS->adfGeoTransform[4] = (tr_y - tl_y)/(poDS->nRasterXSize - 1);
1113 0 : poDS->adfGeoTransform[2] = (bl_x - tl_x)/(poDS->nRasterYSize - 1);
1114 0 : poDS->adfGeoTransform[5] = (bl_y - tl_y)/(poDS->nRasterYSize - 1);
1115 0 : poDS->adfGeoTransform[0] = (tl_x - 0.5*poDS->adfGeoTransform[1]
1116 0 : - 0.5*poDS->adfGeoTransform[2]);
1117 0 : poDS->adfGeoTransform[3] = (tl_y - 0.5*poDS->adfGeoTransform[4]
1118 0 : - 0.5*poDS->adfGeoTransform[5]);
1119 :
1120 : /* Use bottom right pixel to test geotransform */
1121 : br_x = strtod(CPLGetXMLValue( psPos,
1122 0 : "lowerRightCorner.mapCoordinate.easting", "0.0" ), NULL);
1123 : br_y = strtod(CPLGetXMLValue( psPos,
1124 0 : "lowerRightCorner.mapCoordinate.northing", "0.0" ), NULL);
1125 0 : testx = poDS->adfGeoTransform[0] + poDS->adfGeoTransform[1] *
1126 0 : (poDS->nRasterXSize - 0.5) + poDS->adfGeoTransform[2] *
1127 0 : (poDS->nRasterYSize - 0.5);
1128 0 : testy = poDS->adfGeoTransform[3] + poDS->adfGeoTransform[4] *
1129 0 : (poDS->nRasterXSize - 0.5) + poDS->adfGeoTransform[5] *
1130 0 : (poDS->nRasterYSize - 0.5);
1131 :
1132 : /* Give 1/4 pixel numerical error leeway in calculating location
1133 : based on affine transform */
1134 0 : if ( (fabs(testx - br_x) >
1135 0 : fabs(0.25*(poDS->adfGeoTransform[1]+poDS->adfGeoTransform[2])))
1136 0 : || (fabs(testy - br_y) > fabs(0.25*(poDS->adfGeoTransform[4] +
1137 0 : poDS->adfGeoTransform[5]))) )
1138 : {
1139 : CPLError( CE_Warning, CPLE_AppDefined,
1140 : "Unexpected error in calculating affine transform: "
1141 0 : "corner coordinates inconsistent.");
1142 : }
1143 : else
1144 0 : poDS->bHaveGeoTransform = TRUE;
1145 :
1146 : }
1147 :
1148 : }
1149 :
1150 :
1151 : /* -------------------------------------------------------------------- */
1152 : /* Collect Projection String Information */
1153 : /* -------------------------------------------------------------------- */
1154 : CPLXMLNode *psEllipsoid =
1155 : CPLGetXMLNode( psImageAttributes,
1156 4 : "geographicInformation.referenceEllipsoidParameters" );
1157 :
1158 4 : if ( psEllipsoid != NULL ) {
1159 : const char *pszEllipsoidName;
1160 : double minor_axis, major_axis, inv_flattening;
1161 4 : OGRSpatialReference oLL, oPrj;
1162 :
1163 4 : pszEllipsoidName = CPLGetXMLValue( psEllipsoid, "ellipsoidName", "" );
1164 : minor_axis = atof(CPLGetXMLValue( psEllipsoid, "semiMinorAxis",
1165 4 : "0.0" ));
1166 : major_axis = atof(CPLGetXMLValue( psEllipsoid, "semiMajorAxis",
1167 4 : "0.0" ));
1168 :
1169 4 : if ( EQUAL(pszEllipsoidName, "") || ( minor_axis == 0.0 ) ||
1170 : ( major_axis == 0.0 ) )
1171 : {
1172 : CPLError(CE_Warning,CPLE_AppDefined,"Warning- incomplete"
1173 0 : " ellipsoid information. Using wgs-84 parameters.\n");
1174 0 : oLL.SetWellKnownGeogCS( "WGS84" );
1175 0 : oPrj.SetWellKnownGeogCS( "WGS84" );
1176 : }
1177 4 : else if ( EQUAL( pszEllipsoidName, "WGS84" ) ) {
1178 4 : oLL.SetWellKnownGeogCS( "WGS84" );
1179 4 : oPrj.SetWellKnownGeogCS( "WGS84" );
1180 : }
1181 : else {
1182 0 : inv_flattening = major_axis/(major_axis - minor_axis);
1183 : oLL.SetGeogCS( "","",pszEllipsoidName, major_axis,
1184 0 : inv_flattening);
1185 : oPrj.SetGeogCS( "","",pszEllipsoidName, major_axis,
1186 0 : inv_flattening);
1187 : }
1188 :
1189 4 : if ( psMapProjection != NULL ) {
1190 : const char *pszProj = CPLGetXMLValue(
1191 0 : psMapProjection, "mapProjectionDescriptor", "" );
1192 0 : bool bUseProjInfo = FALSE;
1193 :
1194 : CPLXMLNode *psUtmParams =
1195 : CPLGetXMLNode( psMapProjection,
1196 0 : "utmProjectionParameters" );
1197 :
1198 : CPLXMLNode *psNspParams =
1199 : CPLGetXMLNode( psMapProjection,
1200 0 : "nspProjectionParameters" );
1201 :
1202 0 : if ((psUtmParams != NULL) && poDS->bHaveGeoTransform ) {
1203 : const char *pszHemisphere;
1204 : int utmZone;
1205 : double origEasting, origNorthing;
1206 0 : bool bNorth = TRUE;
1207 :
1208 0 : utmZone = atoi(CPLGetXMLValue( psUtmParams, "utmZone", "" ));
1209 : pszHemisphere = CPLGetXMLValue( psUtmParams,
1210 0 : "hemisphere", "" );
1211 : origEasting = strtod(CPLGetXMLValue( psUtmParams,
1212 0 : "mapOriginFalseEasting", "0.0" ), NULL);
1213 : origNorthing = strtod(CPLGetXMLValue( psUtmParams,
1214 0 : "mapOriginFalseNorthing", "0.0" ), NULL);
1215 :
1216 0 : if ( EQUALN(pszHemisphere,"southern",8) )
1217 0 : bNorth = FALSE;
1218 :
1219 0 : if (EQUALN(pszProj,"UTM",3)) {
1220 0 : oPrj.SetUTM(utmZone, bNorth);
1221 0 : bUseProjInfo = TRUE;
1222 : }
1223 : }
1224 0 : else if ((psNspParams != NULL) && poDS->bHaveGeoTransform) {
1225 : double origEasting, origNorthing, copLong, copLat, sP1, sP2;
1226 :
1227 : origEasting = strtod(CPLGetXMLValue( psNspParams,
1228 0 : "mapOriginFalseEasting", "0.0" ), NULL);
1229 : origNorthing = strtod(CPLGetXMLValue( psNspParams,
1230 0 : "mapOriginFalseNorthing", "0.0" ), NULL);
1231 : copLong = strtod(CPLGetXMLValue( psNspParams,
1232 0 : "centerOfProjectionLongitude", "0.0" ), NULL);
1233 : copLat = strtod(CPLGetXMLValue( psNspParams,
1234 0 : "centerOfProjectionLatitude", "0.0" ), NULL);
1235 : sP1 = strtod(CPLGetXMLValue( psNspParams,
1236 0 : "standardParallels1", "0.0" ), NULL);
1237 : sP2 = strtod(CPLGetXMLValue( psNspParams,
1238 0 : "standardParallels2", "0.0" ), NULL);
1239 :
1240 0 : if (EQUALN(pszProj,"ARC",3)) {
1241 : /* Albers Conical Equal Area */
1242 : oPrj.SetACEA(sP1, sP2, copLat, copLong, origEasting,
1243 0 : origNorthing);
1244 0 : bUseProjInfo = TRUE;
1245 : }
1246 0 : else if (EQUALN(pszProj,"LCC",3)) {
1247 : /* Lambert Conformal Conic */
1248 : oPrj.SetLCC(sP1, sP2, copLat, copLong, origEasting,
1249 0 : origNorthing);
1250 0 : bUseProjInfo = TRUE;
1251 : }
1252 0 : else if (EQUALN(pszProj,"STPL",3)) {
1253 : /* State Plate
1254 : ASSUMPTIONS: "zone" in product.xml matches USGS
1255 : definition as expected by ogr spatial reference; NAD83
1256 : zones (versus NAD27) are assumed. */
1257 :
1258 : int nSPZone = atoi(CPLGetXMLValue( psNspParams,
1259 0 : "zone", "1" ));
1260 :
1261 0 : oPrj.SetStatePlane( nSPZone, TRUE, NULL, 0.0 );
1262 0 : bUseProjInfo = TRUE;
1263 : }
1264 : }
1265 :
1266 0 : if (bUseProjInfo) {
1267 0 : CPLFree( poDS->pszProjection );
1268 0 : poDS->pszProjection = NULL;
1269 0 : oPrj.exportToWkt( &(poDS->pszProjection) );
1270 : }
1271 : else {
1272 : CPLError(CE_Warning,CPLE_AppDefined,"Unable to interpret "
1273 : "projection information; check mapProjection info in "
1274 0 : "product.xml!");
1275 : }
1276 : }
1277 :
1278 4 : CPLFree( poDS->pszGCPProjection );
1279 4 : poDS->pszGCPProjection = NULL;
1280 4 : oLL.exportToWkt( &(poDS->pszGCPProjection) );
1281 :
1282 : }
1283 :
1284 : /* -------------------------------------------------------------------- */
1285 : /* Collect GCPs. */
1286 : /* -------------------------------------------------------------------- */
1287 : CPLXMLNode *psGeoGrid =
1288 : CPLGetXMLNode( psImageAttributes,
1289 4 : "geographicInformation.geolocationGrid" );
1290 :
1291 4 : if( psGeoGrid != NULL ) {
1292 : /* count GCPs */
1293 4 : poDS->nGCPCount = 0;
1294 :
1295 20 : for( psNode = psGeoGrid->psChild; psNode != NULL;
1296 : psNode = psNode->psNext )
1297 : {
1298 16 : if( EQUAL(psNode->pszValue,"imageTiePoint") )
1299 16 : poDS->nGCPCount++ ;
1300 : }
1301 :
1302 : poDS->pasGCPList = (GDAL_GCP *)
1303 4 : CPLCalloc(sizeof(GDAL_GCP),poDS->nGCPCount);
1304 :
1305 4 : poDS->nGCPCount = 0;
1306 :
1307 20 : for( psNode = psGeoGrid->psChild; psNode != NULL;
1308 : psNode = psNode->psNext )
1309 : {
1310 : char szID[32];
1311 16 : GDAL_GCP *psGCP = poDS->pasGCPList + poDS->nGCPCount;
1312 :
1313 16 : if( !EQUAL(psNode->pszValue,"imageTiePoint") )
1314 0 : continue;
1315 :
1316 16 : poDS->nGCPCount++ ;
1317 :
1318 16 : sprintf( szID, "%d", poDS->nGCPCount );
1319 16 : psGCP->pszId = CPLStrdup( szID );
1320 16 : psGCP->pszInfo = CPLStrdup("");
1321 : psGCP->dfGCPPixel =
1322 16 : atof(CPLGetXMLValue(psNode,"imageCoordinate.pixel","0"));
1323 : psGCP->dfGCPLine =
1324 16 : atof(CPLGetXMLValue(psNode,"imageCoordinate.line","0"));
1325 : psGCP->dfGCPX =
1326 16 : atof(CPLGetXMLValue(psNode,"geodeticCoordinate.longitude",""));
1327 : psGCP->dfGCPY =
1328 16 : atof(CPLGetXMLValue(psNode,"geodeticCoordinate.latitude",""));
1329 : psGCP->dfGCPZ =
1330 16 : atof(CPLGetXMLValue(psNode,"geodeticCoordinate.height",""));
1331 : }
1332 : }
1333 :
1334 4 : CPLFree( pszPath );
1335 4 : if (pszBeta0LUT) CPLFree(pszBeta0LUT);
1336 4 : if (pszSigma0LUT) CPLFree(pszSigma0LUT);
1337 4 : if (pszGammaLUT) CPLFree(pszGammaLUT);
1338 :
1339 : /* -------------------------------------------------------------------- */
1340 : /* Initialize any PAM information. */
1341 : /* -------------------------------------------------------------------- */
1342 4 : CPLString osDescription;
1343 :
1344 4 : switch (eCalib) {
1345 : case Sigma0:
1346 : osDescription.Printf( "RADARSAT_2_CALIB:SIGMA0:%s",
1347 0 : osMDFilename.c_str() );
1348 0 : break;
1349 : case Beta0:
1350 : osDescription.Printf( "RADARSAT_2_CALIB:BETA0:%s",
1351 2 : osMDFilename.c_str());
1352 2 : break;
1353 : case Gamma:
1354 : osDescription.Printf( "RADARSAT_2_CALIB:GAMMA0:%s",
1355 0 : osMDFilename.c_str() );
1356 0 : break;
1357 : case Uncalib:
1358 : osDescription.Printf( "RADARSAT_2_CALIB:UNCALIB:%s",
1359 0 : osMDFilename.c_str() );
1360 0 : break;
1361 : default:
1362 2 : osDescription = osMDFilename;
1363 : }
1364 :
1365 4 : if( eCalib != None )
1366 : poDS->papszExtraFiles =
1367 2 : CSLAddString( poDS->papszExtraFiles, osMDFilename );
1368 :
1369 : /* -------------------------------------------------------------------- */
1370 : /* Initialize any PAM information. */
1371 : /* -------------------------------------------------------------------- */
1372 4 : poDS->SetDescription( osDescription );
1373 :
1374 4 : poDS->SetPhysicalFilename( osMDFilename );
1375 4 : poDS->SetSubdatasetName( osDescription );
1376 :
1377 4 : poDS->TryLoadXML();
1378 :
1379 : /* -------------------------------------------------------------------- */
1380 : /* Check for overviews. */
1381 : /* -------------------------------------------------------------------- */
1382 4 : poDS->oOvManager.Initialize( poDS, ":::VIRTUAL:::" );
1383 :
1384 4 : return( poDS );
1385 : }
1386 :
1387 : /************************************************************************/
1388 : /* GetGCPCount() */
1389 : /************************************************************************/
1390 :
1391 0 : int RS2Dataset::GetGCPCount()
1392 :
1393 : {
1394 0 : return nGCPCount;
1395 : }
1396 :
1397 : /************************************************************************/
1398 : /* GetGCPProjection() */
1399 : /************************************************************************/
1400 :
1401 0 : const char *RS2Dataset::GetGCPProjection()
1402 :
1403 : {
1404 0 : return pszGCPProjection;
1405 : }
1406 :
1407 : /************************************************************************/
1408 : /* GetGCPs() */
1409 : /************************************************************************/
1410 :
1411 0 : const GDAL_GCP *RS2Dataset::GetGCPs()
1412 :
1413 : {
1414 0 : return pasGCPList;
1415 : }
1416 :
1417 :
1418 : /************************************************************************/
1419 : /* GetProjectionRef() */
1420 : /************************************************************************/
1421 :
1422 0 : const char *RS2Dataset::GetProjectionRef()
1423 :
1424 : {
1425 0 : return( pszProjection );
1426 : }
1427 :
1428 : /************************************************************************/
1429 : /* GetGeoTransform() */
1430 : /************************************************************************/
1431 :
1432 0 : CPLErr RS2Dataset::GetGeoTransform( double * padfTransform )
1433 :
1434 : {
1435 0 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
1436 :
1437 0 : if (bHaveGeoTransform)
1438 0 : return( CE_None );
1439 :
1440 0 : return( CE_Failure );
1441 : }
1442 :
1443 :
1444 : /************************************************************************/
1445 : /* GetMetadata() */
1446 : /************************************************************************/
1447 :
1448 0 : char **RS2Dataset::GetMetadata( const char *pszDomain )
1449 :
1450 : {
1451 0 : if( pszDomain != NULL && EQUALN( pszDomain, "SUBDATASETS", 11 ) &&
1452 : papszSubDatasets != NULL)
1453 0 : return papszSubDatasets;
1454 : else
1455 0 : return GDALDataset::GetMetadata( pszDomain );
1456 : }
1457 :
1458 : /************************************************************************/
1459 : /* GDALRegister_RS2() */
1460 : /************************************************************************/
1461 :
1462 1135 : void GDALRegister_RS2()
1463 :
1464 : {
1465 : GDALDriver *poDriver;
1466 :
1467 1135 : if( GDALGetDriverByName( "RS2" ) == NULL )
1468 : {
1469 1093 : poDriver = new GDALDriver();
1470 :
1471 1093 : poDriver->SetDescription( "RS2" );
1472 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1473 1093 : "RadarSat 2 XML Product" );
1474 1093 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "frmt_rs2.html" );
1475 :
1476 1093 : poDriver->pfnOpen = RS2Dataset::Open;
1477 1093 : poDriver->pfnIdentify = RS2Dataset::Identify;
1478 :
1479 1093 : GetGDALDriverManager()->RegisterDriver( poDriver );
1480 : }
1481 1135 : }
1482 :
|
