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