1 : /******************************************************************************
2 : * $Id: cpgdataset.cpp 18073 2009-11-22 01:01:14Z rouault $
3 : *
4 : * Project: Polarimetric Workstation
5 : * Purpose: Convair PolGASP data (.img/.hdr format).
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 "rawdataset.h"
31 : #include "ogr_spatialref.h"
32 : #include "cpl_string.h"
33 :
34 : CPL_CVSID("$Id: cpgdataset.cpp 18073 2009-11-22 01:01:14Z rouault $");
35 :
36 : CPL_C_START
37 : void GDALRegister_CPG(void);
38 : CPL_C_END
39 :
40 :
41 : enum Interleave {BSQ, BIL, BIP};
42 :
43 : /************************************************************************/
44 : /* ==================================================================== */
45 : /* CPGDataset */
46 : /* ==================================================================== */
47 : /************************************************************************/
48 :
49 : class SIRC_QSLCRasterBand;
50 : class CPG_STOKESRasterBand;
51 :
52 : class CPGDataset : public RawDataset
53 : {
54 : friend class SIRC_QSLCRasterBand;
55 : friend class CPG_STOKESRasterBand;
56 :
57 : FILE *afpImage[4];
58 :
59 : int nGCPCount;
60 : GDAL_GCP *pasGCPList;
61 : char *pszGCPProjection;
62 :
63 : double adfGeoTransform[6];
64 : char *pszProjection;
65 :
66 : int nLoadedStokesLine;
67 : float *padfStokesMatrix;
68 :
69 : int nInterleave;
70 : static int AdjustFilename( char **, const char *, const char * );
71 : static int FindType1( const char *pszWorkname );
72 : static int FindType2( const char *pszWorkname );
73 : static int FindType3( const char *pszWorkname );
74 : static GDALDataset *InitializeType1Or2Dataset( const char *pszWorkname );
75 : static GDALDataset *InitializeType3Dataset( const char *pszWorkname );
76 : CPLErr LoadStokesLine( int iLine, int bNativeOrder );
77 :
78 : public:
79 : CPGDataset();
80 : ~CPGDataset();
81 :
82 : virtual int GetGCPCount();
83 : virtual const char *GetGCPProjection();
84 : virtual const GDAL_GCP *GetGCPs();
85 :
86 : virtual const char *GetProjectionRef(void);
87 : virtual CPLErr GetGeoTransform( double * );
88 :
89 : static GDALDataset *Open( GDALOpenInfo * );
90 : };
91 :
92 : /************************************************************************/
93 : /* CPGDataset() */
94 : /************************************************************************/
95 :
96 1 : CPGDataset::CPGDataset()
97 : {
98 : int iBand;
99 :
100 1 : nGCPCount = 0;
101 1 : pasGCPList = NULL;
102 1 : pszProjection = CPLStrdup("");
103 1 : pszGCPProjection = CPLStrdup("");
104 1 : adfGeoTransform[0] = 0.0;
105 1 : adfGeoTransform[1] = 1.0;
106 1 : adfGeoTransform[2] = 0.0;
107 1 : adfGeoTransform[3] = 0.0;
108 1 : adfGeoTransform[4] = 0.0;
109 1 : adfGeoTransform[5] = 1.0;
110 :
111 1 : nLoadedStokesLine = -1;
112 1 : padfStokesMatrix = NULL;
113 :
114 5 : for( iBand = 0; iBand < 4; iBand++ )
115 4 : afpImage[iBand] = NULL;
116 1 : }
117 :
118 : /************************************************************************/
119 : /* ~CPGDataset() */
120 : /************************************************************************/
121 :
122 1 : CPGDataset::~CPGDataset()
123 :
124 : {
125 : int iBand;
126 :
127 1 : FlushCache();
128 :
129 5 : for( iBand = 0; iBand < 4; iBand++ )
130 : {
131 4 : if( afpImage[iBand] != NULL )
132 1 : VSIFClose( afpImage[iBand] );
133 : }
134 :
135 1 : if( nGCPCount > 0 )
136 : {
137 1 : GDALDeinitGCPs( nGCPCount, pasGCPList );
138 1 : CPLFree( pasGCPList );
139 : }
140 :
141 1 : CPLFree( pszProjection );
142 1 : CPLFree( pszGCPProjection );
143 :
144 1 : if (padfStokesMatrix != NULL)
145 0 : CPLFree( padfStokesMatrix );
146 :
147 1 : }
148 :
149 : /************************************************************************/
150 : /* ==================================================================== */
151 : /* SIRC_QSLCPRasterBand */
152 : /* ==================================================================== */
153 : /************************************************************************/
154 :
155 : class SIRC_QSLCRasterBand : public GDALRasterBand
156 4 : {
157 : friend class CPGDataset;
158 :
159 : public:
160 : SIRC_QSLCRasterBand( CPGDataset *, int, GDALDataType );
161 :
162 : virtual CPLErr IReadBlock( int, int, void * );
163 : };
164 :
165 : #define M11 0
166 : #define M12 1
167 : #define M13 2
168 : #define M14 3
169 : #define M21 4
170 : #define M22 5
171 : #define M23 6
172 : #define M24 7
173 : #define M31 8
174 : #define M32 9
175 : #define M33 10
176 : #define M34 11
177 : #define M41 12
178 : #define M42 13
179 : #define M43 14
180 : #define M44 15
181 :
182 : /************************************************************************/
183 : /* ==================================================================== */
184 : /* CPG_STOKESRasterBand */
185 : /* ==================================================================== */
186 : /************************************************************************/
187 :
188 : class CPG_STOKESRasterBand : public GDALRasterBand
189 : {
190 : friend class CPGDataset;
191 :
192 : int nBand;
193 : int bNativeOrder;
194 :
195 : public:
196 : CPG_STOKESRasterBand( GDALDataset *poDS, int nBand,
197 : GDALDataType eType,
198 : int bNativeOrder );
199 : virtual ~CPG_STOKESRasterBand();
200 :
201 : virtual CPLErr IReadBlock( int, int, void * );
202 : };
203 :
204 : /************************************************************************/
205 : /* AdjustFilename() */
206 : /* */
207 : /* Try to find the file with the request polarization and */
208 : /* extention and update the passed filename accordingly. */
209 : /* */
210 : /* Return TRUE if file found otherwise FALSE. */
211 : /************************************************************************/
212 :
213 4 : int CPGDataset::AdjustFilename( char **pszFilename,
214 : const char *pszPolarization,
215 : const char *pszExtension )
216 :
217 : {
218 : VSIStatBuf sStatBuf;
219 : const char *pszNewName;
220 : char *subptr;
221 :
222 : /* eventually we should handle upper/lower case ... */
223 :
224 4 : if ( EQUAL(pszPolarization,"stokes") )
225 : {
226 : pszNewName = CPLResetExtension((const char *) *pszFilename,
227 0 : (const char *) pszExtension);
228 0 : CPLFree(*pszFilename);
229 0 : *pszFilename = CPLStrdup(pszNewName);
230 : }
231 4 : else if (strlen(pszPolarization) == 2)
232 : {
233 1 : subptr = strstr(*pszFilename,"hh");
234 1 : if (subptr == NULL)
235 1 : subptr = strstr(*pszFilename,"hv");
236 1 : if (subptr == NULL)
237 1 : subptr = strstr(*pszFilename,"vv");
238 1 : if (subptr == NULL)
239 1 : subptr = strstr(*pszFilename,"vh");
240 1 : if (subptr == NULL)
241 1 : return FALSE;
242 :
243 0 : strncpy( subptr, pszPolarization, 2);
244 : pszNewName = CPLResetExtension((const char *) *pszFilename,
245 0 : (const char *) pszExtension);
246 0 : CPLFree(*pszFilename);
247 0 : *pszFilename = CPLStrdup(pszNewName);
248 :
249 : }
250 : else
251 : {
252 : pszNewName = CPLResetExtension((const char *) *pszFilename,
253 3 : (const char *) pszExtension);
254 3 : CPLFree(*pszFilename);
255 3 : *pszFilename = CPLStrdup(pszNewName);
256 : }
257 3 : return VSIStat( *pszFilename, &sStatBuf ) == 0;
258 : }
259 :
260 : /************************************************************************/
261 : /* Search for the various types of Convair filesets */
262 : /* Return TRUE for a match, FALSE for no match */
263 : /************************************************************************/
264 11982 : int CPGDataset::FindType1( const char *pszFilename )
265 : {
266 : int nNameLen;
267 :
268 11982 : nNameLen = strlen(pszFilename);
269 :
270 11982 : if ((strstr(pszFilename,"sso") == NULL) &&
271 : (strstr(pszFilename,"polgasp") == NULL))
272 11982 : return FALSE;
273 :
274 0 : if (( strlen(pszFilename) < 5) ||
275 : (!EQUAL(pszFilename+nNameLen-4,".hdr")
276 : && !EQUAL(pszFilename+nNameLen-4,".img")))
277 0 : return FALSE;
278 :
279 : /* Expect all bands and headers to be present */
280 0 : char* pszTemp = CPLStrdup(pszFilename);
281 :
282 : int bNotFound = !AdjustFilename( &pszTemp, "hh", "img" )
283 : || !AdjustFilename( &pszTemp, "hh", "hdr" )
284 : || !AdjustFilename( &pszTemp, "hv", "img" )
285 : || !AdjustFilename( &pszTemp, "hv", "hdr" )
286 : || !AdjustFilename( &pszTemp, "vh", "img" )
287 : || !AdjustFilename( &pszTemp, "vh", "hdr" )
288 : || !AdjustFilename( &pszTemp, "vv", "img" )
289 0 : || !AdjustFilename( &pszTemp, "vv", "hdr" );
290 :
291 0 : CPLFree(pszTemp);
292 :
293 0 : if (bNotFound)
294 0 : return FALSE;
295 :
296 0 : return TRUE;
297 : }
298 :
299 11982 : int CPGDataset::FindType2( const char *pszFilename )
300 : {
301 : int nNameLen;
302 :
303 11982 : nNameLen = strlen( pszFilename );
304 :
305 11982 : if (( strlen(pszFilename) < 9) ||
306 : (!EQUAL(pszFilename+nNameLen-8,"SIRC.hdr")
307 : && !EQUAL(pszFilename+nNameLen-8,"SIRC.img")))
308 11981 : return FALSE;
309 :
310 1 : char* pszTemp = CPLStrdup(pszFilename);
311 : int bNotFound = !AdjustFilename( &pszTemp, "", "img" )
312 1 : || !AdjustFilename( &pszTemp, "", "hdr" );
313 1 : CPLFree(pszTemp);
314 :
315 1 : if (bNotFound)
316 0 : return FALSE;
317 :
318 1 : return TRUE;
319 : }
320 :
321 0 : int CPGDataset::FindType3( const char *pszFilename )
322 : {
323 : int nNameLen;
324 :
325 0 : nNameLen = strlen( pszFilename );
326 :
327 0 : if ((strstr(pszFilename,"sso") == NULL) &&
328 : (strstr(pszFilename,"polgasp") == NULL))
329 0 : return FALSE;
330 :
331 0 : if (( strlen(pszFilename) < 9) ||
332 : (!EQUAL(pszFilename+nNameLen-4,".img")
333 : && !EQUAL(pszFilename+nNameLen-8,".img_def")))
334 0 : return FALSE;
335 :
336 0 : char* pszTemp = CPLStrdup(pszFilename);
337 : int bNotFound = !AdjustFilename( &pszTemp, "stokes", "img" )
338 0 : || !AdjustFilename( &pszTemp, "stokes", "img_def" );
339 0 : CPLFree(pszTemp);
340 :
341 0 : if (bNotFound)
342 0 : return FALSE;
343 :
344 0 : return TRUE;
345 : }
346 :
347 : /************************************************************************/
348 : /* LoadStokesLine() */
349 : /************************************************************************/
350 :
351 0 : CPLErr CPGDataset::LoadStokesLine( int iLine, int bNativeOrder )
352 :
353 : {
354 : int offset, nBytesToRead, band_index;
355 0 : int nDataSize = GDALGetDataTypeSize(GDT_Float32)/8;
356 :
357 0 : if( iLine == nLoadedStokesLine )
358 0 : return CE_None;
359 :
360 : /* -------------------------------------------------------------------- */
361 : /* allocate working buffers if we don't have them already. */
362 : /* -------------------------------------------------------------------- */
363 0 : if( padfStokesMatrix == NULL )
364 : {
365 0 : padfStokesMatrix = (float *) CPLMalloc(sizeof(float) * nRasterXSize*16);
366 : }
367 :
368 : /* -------------------------------------------------------------------- */
369 : /* Load all the pixel data associated with this scanline. */
370 : /* Retains same interleaving as original dataset. */
371 : /* -------------------------------------------------------------------- */
372 0 : if ( nInterleave == BIP )
373 : {
374 0 : offset = nRasterXSize*iLine*nDataSize*16;
375 0 : nBytesToRead = nDataSize*nRasterXSize*16;
376 0 : if (( VSIFSeek( afpImage[0], offset, SEEK_SET ) != 0 ) ||
377 : (int) VSIFRead( ( GByte *) padfStokesMatrix, 1, nBytesToRead,
378 : afpImage[0] ) != nBytesToRead )
379 : {
380 : CPLError( CE_Failure, CPLE_FileIO,
381 : "Error reading %d bytes of Stokes Convair at offset %d.\n"
382 : "Reading file %s failed.",
383 0 : nBytesToRead, offset, GetDescription() );
384 0 : CPLFree( padfStokesMatrix );
385 0 : padfStokesMatrix = NULL;
386 0 : nLoadedStokesLine = -1;
387 0 : return CE_Failure;
388 : }
389 : }
390 0 : else if ( nInterleave == BIL )
391 : {
392 0 : for ( band_index = 0; band_index < 16; band_index++)
393 : {
394 : offset = nDataSize * (nRasterXSize*iLine +
395 0 : nRasterXSize*band_index);
396 0 : nBytesToRead = nDataSize*nRasterXSize;
397 0 : if (( VSIFSeek( afpImage[0], offset, SEEK_SET ) != 0 ) ||
398 : (int) VSIFRead(
399 : ( GByte *) padfStokesMatrix + nBytesToRead*band_index,
400 : 1, nBytesToRead,
401 : afpImage[0] ) != nBytesToRead )
402 : {
403 : CPLError( CE_Failure, CPLE_FileIO,
404 : "Error reading %d bytes of Stokes Convair at offset %d.\n"
405 : "Reading file %s failed.",
406 0 : nBytesToRead, offset, GetDescription() );
407 0 : CPLFree( padfStokesMatrix );
408 0 : padfStokesMatrix = NULL;
409 0 : nLoadedStokesLine = -1;
410 0 : return CE_Failure;
411 :
412 : }
413 : }
414 : }
415 : else
416 : {
417 0 : for ( band_index = 0; band_index < 16; band_index++)
418 : {
419 : offset = nDataSize * (nRasterXSize*iLine +
420 0 : nRasterXSize*nRasterYSize*band_index);
421 0 : nBytesToRead = nDataSize*nRasterXSize;
422 0 : if (( VSIFSeek( afpImage[0], offset, SEEK_SET ) != 0 ) ||
423 : (int) VSIFRead(
424 : ( GByte *) padfStokesMatrix + nBytesToRead*band_index,
425 : 1, nBytesToRead,
426 : afpImage[0] ) != nBytesToRead )
427 : {
428 : CPLError( CE_Failure, CPLE_FileIO,
429 : "Error reading %d bytes of Stokes Convair at offset %d.\n"
430 : "Reading file %s failed.",
431 0 : nBytesToRead, offset, GetDescription() );
432 0 : CPLFree( padfStokesMatrix );
433 0 : padfStokesMatrix = NULL;
434 0 : nLoadedStokesLine = -1;
435 0 : return CE_Failure;
436 :
437 : }
438 : }
439 : }
440 :
441 0 : if (!bNativeOrder)
442 0 : GDALSwapWords( padfStokesMatrix,nDataSize,nRasterXSize*16, nDataSize );
443 :
444 0 : nLoadedStokesLine = iLine;
445 :
446 0 : return CE_None;
447 : }
448 :
449 : /************************************************************************/
450 : /* Parse header information and initialize dataset for the */
451 : /* appropriate Convair dataset style. */
452 : /* Returns dataset if successful; NULL if there was a problem. */
453 : /************************************************************************/
454 :
455 1 : GDALDataset* CPGDataset::InitializeType1Or2Dataset( const char *pszFilename )
456 : {
457 :
458 : /* -------------------------------------------------------------------- */
459 : /* Read the .hdr file (the hh one for the .sso and polgasp cases) */
460 : /* and parse it. */
461 : /* -------------------------------------------------------------------- */
462 : char **papszHdrLines;
463 : int iLine;
464 1 : int nLines = 0, nSamples = 0;
465 1 : int nError = 0;
466 1 : int nNameLen = 0;
467 :
468 : /* Parameters required for pseudo-geocoding. GCPs map */
469 : /* slant range to ground range at 16 points. */
470 1 : int iGeoParamsFound = 0, itransposed = 0;
471 1 : double dfaltitude = 0.0, dfnear_srd = 0.0;
472 1 : double dfsample_size = 0.0, dfsample_size_az = 0.0;
473 :
474 : /* Parameters in geogratis geocoded images */
475 1 : int iUTMParamsFound = 0, iUTMZone=0, iCorner=0;
476 1 : double dfnorth = 0.0, dfeast = 0.0;
477 :
478 1 : char* pszWorkname = CPLStrdup(pszFilename);
479 1 : AdjustFilename( &pszWorkname, "hh", "hdr" );
480 1 : papszHdrLines = CSLLoad( pszWorkname );
481 :
482 8 : for( iLine = 0; papszHdrLines && papszHdrLines[iLine] != NULL; iLine++ )
483 : {
484 7 : char **papszTokens = CSLTokenizeString( papszHdrLines[iLine] );
485 :
486 : /* Note: some cv580 file seem to have comments with #, hence the >=
487 : * instead of = for token checking, and the equalN for the corner.
488 : */
489 :
490 7 : if( CSLCount( papszTokens ) < 2 )
491 : {
492 : /* ignore */;
493 : }
494 7 : else if ( ( CSLCount( papszTokens ) >= 3 ) &&
495 0 : EQUAL(papszTokens[0],"reference") &&
496 0 : EQUAL(papszTokens[1],"north") )
497 : {
498 0 : dfnorth = atof(papszTokens[2]);
499 0 : iUTMParamsFound++;
500 : }
501 7 : else if ( ( CSLCount( papszTokens ) >= 3 ) &&
502 0 : EQUAL(papszTokens[0],"reference") &&
503 0 : EQUAL(papszTokens[1],"east") )
504 : {
505 0 : dfeast = atof(papszTokens[2]);
506 0 : iUTMParamsFound++;
507 : }
508 7 : else if ( ( CSLCount( papszTokens ) >= 5 ) &&
509 0 : EQUAL(papszTokens[0],"reference") &&
510 0 : EQUAL(papszTokens[1],"projection") &&
511 0 : EQUAL(papszTokens[2],"UTM") &&
512 0 : EQUAL(papszTokens[3],"zone") )
513 : {
514 0 : iUTMZone = atoi(papszTokens[4]);
515 0 : iUTMParamsFound++;
516 : }
517 7 : else if ( ( CSLCount( papszTokens ) >= 3 ) &&
518 0 : EQUAL(papszTokens[0],"reference") &&
519 0 : EQUAL(papszTokens[1],"corner") &&
520 0 : EQUALN(papszTokens[2],"Upper_Left",10) )
521 : {
522 0 : iCorner = 0;
523 0 : iUTMParamsFound++;
524 : }
525 7 : else if( EQUAL(papszTokens[0],"number_lines") )
526 1 : nLines = atoi(papszTokens[1]);
527 :
528 6 : else if( EQUAL(papszTokens[0],"number_samples") )
529 1 : nSamples = atoi(papszTokens[1]);
530 :
531 20 : else if( (EQUAL(papszTokens[0],"header_offset")
532 0 : && atoi(papszTokens[1]) != 0)
533 5 : || (EQUAL(papszTokens[0],"number_channels")
534 0 : && (atoi(papszTokens[1]) != 1)
535 0 : && (atoi(papszTokens[1]) != 10))
536 5 : || (EQUAL(papszTokens[0],"datatype")
537 0 : && atoi(papszTokens[1]) != 1)
538 5 : || (EQUAL(papszTokens[0],"number_format")
539 0 : && !EQUAL(papszTokens[1],"float32")
540 0 : && !EQUAL(papszTokens[1],"int8")))
541 : {
542 : CPLError( CE_Failure, CPLE_AppDefined,
543 : "Keyword %s has value %s which does not match CPG driver expectation.",
544 0 : papszTokens[0], papszTokens[1] );
545 0 : nError = 1;
546 : }
547 5 : else if( EQUAL(papszTokens[0],"altitude") )
548 : {
549 1 : dfaltitude = atof(papszTokens[1]);
550 1 : iGeoParamsFound++;
551 : }
552 4 : else if( EQUAL(papszTokens[0],"near_srd") )
553 : {
554 1 : dfnear_srd = atof(papszTokens[1]);
555 1 : iGeoParamsFound++;
556 : }
557 :
558 3 : else if( EQUAL(papszTokens[0],"sample_size") )
559 : {
560 1 : dfsample_size = atof(papszTokens[1]);
561 1 : iGeoParamsFound++;
562 1 : iUTMParamsFound++;
563 : }
564 2 : else if( EQUAL(papszTokens[0],"sample_size_az") )
565 : {
566 1 : dfsample_size_az = atof(papszTokens[1]);
567 1 : iGeoParamsFound++;
568 1 : iUTMParamsFound++;
569 : }
570 1 : else if( EQUAL(papszTokens[0],"transposed") )
571 : {
572 1 : itransposed = atoi(papszTokens[1]);
573 1 : iGeoParamsFound++;
574 1 : iUTMParamsFound++;
575 : }
576 :
577 :
578 :
579 7 : CSLDestroy( papszTokens );
580 : }
581 1 : CSLDestroy( papszHdrLines );
582 : /* -------------------------------------------------------------------- */
583 : /* Check for successful completion. */
584 : /* -------------------------------------------------------------------- */
585 1 : if( nError )
586 : {
587 0 : CPLFree(pszWorkname);
588 0 : return NULL;
589 : }
590 :
591 1 : if( nLines <= 0 || nSamples <= 0 )
592 : {
593 : CPLError( CE_Failure, CPLE_AppDefined,
594 : "Did not find valid number_lines or number_samples keywords in %s.",
595 0 : pszWorkname );
596 0 : CPLFree(pszWorkname);
597 0 : return NULL;
598 : }
599 :
600 : /* -------------------------------------------------------------------- */
601 : /* Initialize dataset. */
602 : /* -------------------------------------------------------------------- */
603 1 : int iBand=0;
604 : CPGDataset *poDS;
605 :
606 1 : poDS = new CPGDataset();
607 :
608 1 : poDS->nRasterXSize = nSamples;
609 1 : poDS->nRasterYSize = nLines;
610 :
611 : /* -------------------------------------------------------------------- */
612 : /* Open the four bands. */
613 : /* -------------------------------------------------------------------- */
614 : static const char *apszPolarizations[4] = { "hh", "hv", "vv", "vh" };
615 :
616 1 : nNameLen = strlen(pszWorkname);
617 :
618 2 : if ( EQUAL(pszWorkname+nNameLen-7,"IRC.hdr") ||
619 : EQUAL(pszWorkname+nNameLen-7,"IRC.img") )
620 : {
621 :
622 1 : AdjustFilename( &pszWorkname, "" , "img" );
623 1 : poDS->afpImage[0] = VSIFOpen( pszWorkname, "rb" );
624 1 : if( poDS->afpImage[0] == NULL )
625 : {
626 : CPLError( CE_Failure, CPLE_OpenFailed,
627 : "Failed to open .img file: %s",
628 0 : pszWorkname );
629 0 : CPLFree(pszWorkname);
630 0 : delete poDS;
631 0 : return NULL;
632 : }
633 10 : for( iBand = 0; iBand < 4; iBand++ )
634 : {
635 : SIRC_QSLCRasterBand *poBand;
636 :
637 4 : poBand = new SIRC_QSLCRasterBand( poDS, iBand+1, GDT_CFloat32 );
638 4 : poDS->SetBand( iBand+1, poBand );
639 : poBand->SetMetadataItem( "POLARIMETRIC_INTERP",
640 4 : apszPolarizations[iBand] );
641 : }
642 : }
643 : else
644 : {
645 0 : for( iBand = 0; iBand < 4; iBand++ )
646 : {
647 : RawRasterBand *poBand;
648 :
649 0 : AdjustFilename( &pszWorkname, apszPolarizations[iBand], "img" );
650 :
651 0 : poDS->afpImage[iBand] = VSIFOpen( pszWorkname, "rb" );
652 0 : if( poDS->afpImage[iBand] == NULL )
653 : {
654 : CPLError( CE_Failure, CPLE_OpenFailed,
655 : "Failed to open .img file: %s",
656 0 : pszWorkname );
657 0 : CPLFree(pszWorkname);
658 0 : delete poDS;
659 0 : return NULL;
660 : }
661 :
662 : poBand =
663 0 : new RawRasterBand( poDS, iBand+1, poDS->afpImage[iBand],
664 : 0, 8, 8*nSamples,
665 0 : GDT_CFloat32, !CPL_IS_LSB, FALSE );
666 0 : poDS->SetBand( iBand+1, poBand );
667 :
668 : poBand->SetMetadataItem( "POLARIMETRIC_INTERP",
669 0 : apszPolarizations[iBand] );
670 : }
671 : }
672 :
673 : /* Set an appropriate matrix representation metadata item for the set */
674 1 : if ( poDS->GetRasterCount() == 4 ) {
675 1 : poDS->SetMetadataItem( "MATRIX_REPRESENTATION", "SCATTERING" );
676 : }
677 :
678 : /* ------------------------------------------------------------------------- */
679 : /* Add georeferencing or pseudo-geocoding, if enough information found. */
680 : /* ------------------------------------------------------------------------- */
681 1 : if (iUTMParamsFound == 7)
682 : {
683 0 : OGRSpatialReference oUTM;
684 : double dfnorth_center;
685 :
686 0 : poDS->adfGeoTransform[1] = 0.0;
687 0 : poDS->adfGeoTransform[2] = 0.0;
688 0 : poDS->adfGeoTransform[4] = 0.0;
689 0 : poDS->adfGeoTransform[5] = 0.0;
690 :
691 0 : if (itransposed == 1)
692 : {
693 : printf("Warning- did not have a convair SIRC-style test dataset\n"
694 0 : "with transposed=1 for testing. Georefencing may be wrong.\n");
695 0 : dfnorth_center = dfnorth - nSamples*dfsample_size/2.0;
696 0 : poDS->adfGeoTransform[0] = dfeast;
697 0 : poDS->adfGeoTransform[2] = dfsample_size_az;
698 0 : poDS->adfGeoTransform[3] = dfnorth;
699 0 : poDS->adfGeoTransform[4] = -1*dfsample_size;
700 : }
701 : else
702 : {
703 0 : dfnorth_center = dfnorth - nLines*dfsample_size/2.0;
704 0 : poDS->adfGeoTransform[0] = dfeast;
705 0 : poDS->adfGeoTransform[1] = dfsample_size_az;
706 0 : poDS->adfGeoTransform[3] = dfnorth;
707 0 : poDS->adfGeoTransform[5] = -1*dfsample_size;
708 : }
709 0 : if (dfnorth_center < 0)
710 0 : oUTM.SetUTM(iUTMZone, 0);
711 : else
712 0 : oUTM.SetUTM(iUTMZone, 1);
713 :
714 : /* Assuming WGS84 */
715 0 : oUTM.SetWellKnownGeogCS( "WGS84" );
716 0 : CPLFree( poDS->pszProjection );
717 0 : poDS->pszProjection = NULL;
718 0 : oUTM.exportToWkt( &(poDS->pszProjection) );
719 :
720 :
721 :
722 : }
723 1 : else if (iGeoParamsFound == 5)
724 : {
725 : int ngcp;
726 : double dfgcpLine, dfgcpPixel, dfgcpX, dfgcpY, dftemp;
727 :
728 1 : poDS->nGCPCount = 16;
729 1 : poDS->pasGCPList = (GDAL_GCP *) CPLCalloc(sizeof(GDAL_GCP),poDS->nGCPCount);
730 1 : GDALInitGCPs(poDS->nGCPCount, poDS->pasGCPList);
731 :
732 17 : for( ngcp = 0; ngcp < 16; ngcp ++ )
733 : {
734 : char szID[32];
735 :
736 16 : sprintf(szID,"%d",ngcp+1);
737 16 : if (itransposed == 1)
738 : {
739 0 : if (ngcp < 4)
740 0 : dfgcpPixel = 0.0;
741 0 : else if (ngcp < 8)
742 0 : dfgcpPixel = nSamples/3.0;
743 0 : else if (ngcp < 12)
744 0 : dfgcpPixel = 2.0*nSamples/3.0;
745 : else
746 0 : dfgcpPixel = nSamples;
747 :
748 0 : dfgcpLine = nLines*( ngcp % 4 )/3.0;
749 :
750 0 : dftemp = dfnear_srd + (dfsample_size*dfgcpLine);
751 : /* -1 so that 0,0 maps to largest Y */
752 0 : dfgcpY = -1*sqrt( dftemp*dftemp - dfaltitude*dfaltitude );
753 0 : dfgcpX = dfgcpPixel*dfsample_size_az;
754 :
755 : }
756 : else
757 : {
758 16 : if (ngcp < 4)
759 4 : dfgcpLine = 0.0;
760 12 : else if (ngcp < 8)
761 4 : dfgcpLine = nLines/3.0;
762 8 : else if (ngcp < 12)
763 4 : dfgcpLine = 2.0*nLines/3.0;
764 : else
765 4 : dfgcpLine = nLines;
766 :
767 16 : dfgcpPixel = nSamples*( ngcp % 4 )/3.0;
768 :
769 16 : dftemp = dfnear_srd + (dfsample_size*dfgcpPixel);
770 16 : dfgcpX = sqrt( dftemp*dftemp - dfaltitude*dfaltitude );
771 16 : dfgcpY = (nLines - dfgcpLine)*dfsample_size_az;
772 :
773 : }
774 16 : poDS->pasGCPList[ngcp].dfGCPX = dfgcpX;
775 16 : poDS->pasGCPList[ngcp].dfGCPY = dfgcpY;
776 16 : poDS->pasGCPList[ngcp].dfGCPZ = 0.0;
777 :
778 16 : poDS->pasGCPList[ngcp].dfGCPPixel = dfgcpPixel;
779 16 : poDS->pasGCPList[ngcp].dfGCPLine = dfgcpLine;
780 :
781 16 : CPLFree(poDS->pasGCPList[ngcp].pszId);
782 16 : poDS->pasGCPList[ngcp].pszId = CPLStrdup( szID );
783 :
784 : }
785 :
786 1 : CPLFree(poDS->pszGCPProjection);
787 1 : poDS->pszGCPProjection = CPLStrdup("LOCAL_CS[\"Ground range view / unreferenced meters\",UNIT[\"Meter\",1.0]]");
788 :
789 : }
790 :
791 1 : CPLFree(pszWorkname);
792 :
793 1 : return poDS;
794 : }
795 :
796 0 : GDALDataset *CPGDataset::InitializeType3Dataset( const char *pszFilename )
797 : {
798 :
799 : char **papszHdrLines;
800 0 : int iLine, iBytesPerPixel = 0, iInterleave=-1;
801 0 : int nLines = 0, nSamples = 0, nBands = 0;
802 0 : int nError = 0;
803 :
804 : /* Parameters in geogratis geocoded images */
805 0 : int iUTMParamsFound = 0, iUTMZone=0;
806 0 : double dfnorth = 0.0, dfeast = 0.0, dfOffsetX = 0.0, dfOffsetY = 0.0;
807 0 : double dfxsize = 0.0, dfysize = 0.0;
808 :
809 0 : char* pszWorkname = CPLStrdup(pszFilename);
810 0 : AdjustFilename( &pszWorkname, "stokes", "img_def" );
811 0 : papszHdrLines = CSLLoad( pszWorkname );
812 :
813 0 : for( iLine = 0; papszHdrLines && papszHdrLines[iLine] != NULL; iLine++ )
814 : {
815 0 : char **papszTokens = CSLTokenizeString2( papszHdrLines[iLine],
816 : " \t",
817 0 : CSLT_HONOURSTRINGS & CSLT_ALLOWEMPTYTOKENS );
818 :
819 : /* Note: some cv580 file seem to have comments with #, hence the >=
820 : * instead of = for token checking, and the equalN for the corner.
821 : */
822 :
823 0 : if ( ( CSLCount( papszTokens ) >= 3 ) &&
824 0 : EQUAL(papszTokens[0],"data") &&
825 0 : EQUAL(papszTokens[1],"organization:"))
826 : {
827 :
828 0 : if( EQUALN(papszTokens[2], "BSQ", 3) )
829 0 : iInterleave = BSQ;
830 0 : else if( EQUALN(papszTokens[2], "BIL", 3) )
831 0 : iInterleave = BIL;
832 0 : else if( EQUALN(papszTokens[2], "BIP", 3) )
833 0 : iInterleave = BIP;
834 : else
835 : {
836 : CPLError( CE_Failure, CPLE_AppDefined,
837 : "The interleaving type of the file (%s) is not supported.",
838 0 : papszTokens[2] );
839 0 : nError = 1;
840 : }
841 :
842 : }
843 0 : else if ( ( CSLCount( papszTokens ) >= 3 ) &&
844 0 : EQUAL(papszTokens[0],"data") &&
845 0 : EQUAL(papszTokens[1],"state:") )
846 : {
847 :
848 0 : if( !EQUALN(papszTokens[2], "RAW", 3) &&
849 0 : !EQUALN(papszTokens[2], "GEO", 3) )
850 : {
851 : CPLError( CE_Failure, CPLE_AppDefined,
852 : "The data state of the file (%s) is not supported.\n. Only RAW and GEO are currently recognized.",
853 0 : papszTokens[2] );
854 0 : nError = 1;
855 : }
856 :
857 :
858 : }
859 0 : else if ( ( CSLCount( papszTokens ) >= 4 ) &&
860 0 : EQUAL(papszTokens[0],"data") &&
861 0 : EQUAL(papszTokens[1],"origin") &&
862 0 : EQUAL(papszTokens[2],"point:") )
863 : {
864 0 : if (!EQUALN(papszTokens[3], "Upper_Left", 10))
865 : {
866 : CPLError( CE_Failure, CPLE_AppDefined,
867 : "Unexpected value (%s) for data origin point- expect Upper_Left.",
868 0 : papszTokens[3] );
869 0 : nError = 1;
870 : }
871 0 : iUTMParamsFound++;
872 : }
873 0 : else if ( ( CSLCount( papszTokens ) >= 5 ) &&
874 0 : EQUAL(papszTokens[0],"map") &&
875 0 : EQUAL(papszTokens[1],"projection:") &&
876 0 : EQUAL(papszTokens[2],"UTM") &&
877 0 : EQUAL(papszTokens[3],"zone") )
878 : {
879 0 : iUTMZone = atoi(papszTokens[4]);
880 0 : iUTMParamsFound++;
881 : }
882 0 : else if ( ( CSLCount( papszTokens ) >= 4 ) &&
883 0 : EQUAL(papszTokens[0],"project") &&
884 0 : EQUAL(papszTokens[1],"origin:") )
885 : {
886 0 : dfeast = atof(papszTokens[2]);
887 0 : dfnorth = atof(papszTokens[3]);
888 0 : iUTMParamsFound+=2;
889 : }
890 0 : else if ( ( CSLCount( papszTokens ) >= 4 ) &&
891 0 : EQUAL(papszTokens[0],"file") &&
892 0 : EQUAL(papszTokens[1],"start:"))
893 : {
894 0 : dfOffsetX = atof(papszTokens[2]);
895 0 : dfOffsetY = atof(papszTokens[3]);
896 0 : iUTMParamsFound+=2;
897 : }
898 0 : else if ( ( CSLCount( papszTokens ) >= 6 ) &&
899 0 : EQUAL(papszTokens[0],"pixel") &&
900 0 : EQUAL(papszTokens[1],"size") &&
901 0 : EQUAL(papszTokens[2],"on") &&
902 0 : EQUAL(papszTokens[3],"ground:"))
903 : {
904 0 : dfxsize = atof(papszTokens[4]);
905 0 : dfysize = atof(papszTokens[5]);
906 0 : iUTMParamsFound+=2;
907 :
908 : }
909 0 : else if ( ( CSLCount( papszTokens ) >= 4 ) &&
910 0 : EQUAL(papszTokens[0],"number") &&
911 0 : EQUAL(papszTokens[1],"of") &&
912 0 : EQUAL(papszTokens[2],"pixels:"))
913 : {
914 0 : nSamples = atoi(papszTokens[3]);
915 : }
916 0 : else if ( ( CSLCount( papszTokens ) >= 4 ) &&
917 0 : EQUAL(papszTokens[0],"number") &&
918 0 : EQUAL(papszTokens[1],"of") &&
919 0 : EQUAL(papszTokens[2],"lines:"))
920 : {
921 0 : nLines = atoi(papszTokens[3]);
922 : }
923 0 : else if ( ( CSLCount( papszTokens ) >= 4 ) &&
924 0 : EQUAL(papszTokens[0],"number") &&
925 0 : EQUAL(papszTokens[1],"of") &&
926 0 : EQUAL(papszTokens[2],"bands:"))
927 : {
928 0 : nBands = atoi(papszTokens[3]);
929 0 : if ( nBands != 16)
930 : {
931 : CPLError( CE_Failure, CPLE_AppDefined,
932 : "Number of bands has a value %s which does not match CPG driver\nexpectation (expect a value of 16).",
933 0 : papszTokens[3] );
934 0 : nError = 1;
935 : }
936 : }
937 0 : else if ( ( CSLCount( papszTokens ) >= 4 ) &&
938 0 : EQUAL(papszTokens[0],"bytes") &&
939 0 : EQUAL(papszTokens[1],"per") &&
940 0 : EQUAL(papszTokens[2],"pixel:"))
941 : {
942 0 : iBytesPerPixel = atoi(papszTokens[3]);
943 0 : if (iBytesPerPixel != 4)
944 : {
945 : CPLError( CE_Failure, CPLE_AppDefined,
946 : "Bytes per pixel has a value %s which does not match CPG driver\nexpectation (expect a value of 4).",
947 0 : papszTokens[1] );
948 0 : nError = 1;
949 : }
950 : }
951 0 : CSLDestroy( papszTokens );
952 : }
953 :
954 0 : CSLDestroy( papszHdrLines );
955 :
956 : /* -------------------------------------------------------------------- */
957 : /* Check for successful completion. */
958 : /* -------------------------------------------------------------------- */
959 0 : if( nError )
960 : {
961 0 : CPLFree(pszWorkname);
962 0 : return NULL;
963 : }
964 :
965 0 : if (!GDALCheckDatasetDimensions(nSamples, nLines) ||
966 : !GDALCheckBandCount(nBands, FALSE) || iInterleave == -1 ||
967 : iBytesPerPixel == 0 )
968 : {
969 : CPLError( CE_Failure, CPLE_AppDefined,
970 : "%s is missing a required parameter (number of pixels, number of lines,\nnumber of bands, bytes per pixel, or data organization).",
971 0 : pszWorkname );
972 0 : CPLFree(pszWorkname);
973 0 : return NULL;
974 : }
975 :
976 : /* -------------------------------------------------------------------- */
977 : /* Initialize dataset. */
978 : /* -------------------------------------------------------------------- */
979 0 : int iBand=0;
980 : CPGDataset *poDS;
981 :
982 0 : poDS = new CPGDataset();
983 :
984 0 : poDS->nRasterXSize = nSamples;
985 0 : poDS->nRasterYSize = nLines;
986 :
987 0 : if( iInterleave == BSQ )
988 0 : poDS->nInterleave = BSQ;
989 0 : else if( iInterleave == BIL )
990 0 : poDS->nInterleave = BIL;
991 : else
992 0 : poDS->nInterleave = BIP;
993 :
994 : /* -------------------------------------------------------------------- */
995 : /* Open the 16 bands. */
996 : /* -------------------------------------------------------------------- */
997 :
998 0 : AdjustFilename( &pszWorkname, "stokes" , "img" );
999 0 : poDS->afpImage[0] = VSIFOpen( pszWorkname, "rb" );
1000 0 : if( poDS->afpImage[0] == NULL )
1001 : {
1002 : CPLError( CE_Failure, CPLE_OpenFailed,
1003 : "Failed to open .img file: %s",
1004 0 : pszWorkname );
1005 0 : CPLFree(pszWorkname);
1006 0 : delete poDS;
1007 0 : return NULL;
1008 : }
1009 0 : for( iBand = 0; iBand < 16; iBand++ )
1010 : {
1011 : CPG_STOKESRasterBand *poBand;
1012 :
1013 : poBand = new CPG_STOKESRasterBand( poDS, iBand+1, GDT_CFloat32,
1014 0 : !CPL_IS_LSB );
1015 0 : poDS->SetBand( iBand+1, poBand );
1016 : }
1017 :
1018 : /* -------------------------------------------------------------------- */
1019 : /* Set appropriate MATRIX_REPRESENTATION. */
1020 : /* -------------------------------------------------------------------- */
1021 0 : if ( poDS->GetRasterCount() == 6 ) {
1022 : poDS->SetMetadataItem( "MATRIX_REPRESENTATION",
1023 0 : "COVARIANCE" );
1024 : }
1025 :
1026 :
1027 : /* ------------------------------------------------------------------------- */
1028 : /* Add georeferencing, if enough information found. */
1029 : /* ------------------------------------------------------------------------- */
1030 0 : if (iUTMParamsFound == 8)
1031 : {
1032 0 : OGRSpatialReference oUTM;
1033 : double dfnorth_center;
1034 :
1035 :
1036 0 : dfnorth_center = dfnorth - nLines*dfysize/2.0;
1037 0 : poDS->adfGeoTransform[0] = dfeast + dfOffsetX;
1038 0 : poDS->adfGeoTransform[1] = dfxsize;
1039 0 : poDS->adfGeoTransform[2] = 0.0;
1040 0 : poDS->adfGeoTransform[3] = dfnorth + dfOffsetY;
1041 0 : poDS->adfGeoTransform[4] = 0.0;
1042 0 : poDS->adfGeoTransform[5] = -1*dfysize;
1043 :
1044 0 : if (dfnorth_center < 0)
1045 0 : oUTM.SetUTM(iUTMZone, 0);
1046 : else
1047 0 : oUTM.SetUTM(iUTMZone, 1);
1048 :
1049 : /* Assuming WGS84 */
1050 0 : oUTM.SetWellKnownGeogCS( "WGS84" );
1051 0 : CPLFree( poDS->pszProjection );
1052 0 : poDS->pszProjection = NULL;
1053 0 : oUTM.exportToWkt( &(poDS->pszProjection) );
1054 : }
1055 :
1056 0 : return poDS;
1057 : }
1058 :
1059 : /************************************************************************/
1060 : /* Open() */
1061 : /************************************************************************/
1062 :
1063 11982 : GDALDataset *CPGDataset::Open( GDALOpenInfo * poOpenInfo )
1064 :
1065 : {
1066 : /* -------------------------------------------------------------------- */
1067 : /* Is this a PolGASP fileset? We expect fileset to follow */
1068 : /* one of these patterns: */
1069 : /* 1) <stuff>hh<stuff2>.img, <stuff>hh<stuff2>.hdr, */
1070 : /* <stuff>hv<stuff2>.img, <stuff>hv<stuff2>.hdr, */
1071 : /* <stuff>vh<stuff2>.img, <stuff>vh<stuff2>.hdr, */
1072 : /* <stuff>vv<stuff2>.img, <stuff>vv<stuff2>.hdr, */
1073 : /* where <stuff> or <stuff2> should contain the */
1074 : /* substring "sso" or "polgasp" */
1075 : /* 2) <stuff>SIRC.hdr and <stuff>SIRC.img */
1076 : /* 3) <stuff>.img and <stuff>.img_def */
1077 : /* where <stuff> should contain the */
1078 : /* substring "sso" or "polgasp" */
1079 : /* -------------------------------------------------------------------- */
1080 11982 : int nNameLen = strlen(poOpenInfo->pszFilename);
1081 11982 : int CPGType = 0;
1082 :
1083 11982 : if ( FindType1( poOpenInfo->pszFilename ))
1084 0 : CPGType = 1;
1085 11982 : else if ( FindType2( poOpenInfo->pszFilename ))
1086 1 : CPGType = 2;
1087 :
1088 : /* Stokes matrix convair data: not quite working yet- something
1089 : * is wrong in the interpretation of the matrix elements in terms
1090 : * of hh, hv, vv, vh. Data will load if the next two lines are
1091 : * uncommented, but values will be incorrect. Expect C11 = hh*conj(hh),
1092 : * C12 = hh*conj(hv), etc. Used geogratis data in both scattering
1093 : * matrix and stokes format for comparison.
1094 : */
1095 : //else if ( FindType3( poOpenInfo->pszFilename ))
1096 : // CPGType = 3;
1097 :
1098 : /* Set working name back to original */
1099 :
1100 11982 : if ( CPGType == 0 )
1101 : {
1102 11981 : nNameLen = strlen(poOpenInfo->pszFilename);
1103 11981 : if ( (nNameLen > 8) &&
1104 : ( ( strstr(poOpenInfo->pszFilename,"sso") != NULL ) ||
1105 : ( strstr(poOpenInfo->pszFilename,"polgasp") != NULL ) ) &&
1106 : ( EQUAL(poOpenInfo->pszFilename+nNameLen-4,"img") ||
1107 : EQUAL(poOpenInfo->pszFilename+nNameLen-4,"hdr") ||
1108 : EQUAL(poOpenInfo->pszFilename+nNameLen-7,"img_def") ) )
1109 : {
1110 : CPLError( CE_Failure, CPLE_OpenFailed,
1111 : "Apparent attempt to open Convair PolGASP data failed as\n"
1112 : "one or more of the required files is missing (eight files\n"
1113 0 : "are expected for scattering matrix format, two for Stokes)." );
1114 : }
1115 11981 : else if ( (nNameLen > 8) &&
1116 : ( strstr(poOpenInfo->pszFilename,"SIRC") != NULL ) &&
1117 : ( EQUAL(poOpenInfo->pszFilename+nNameLen-4,"img") ||
1118 : EQUAL(poOpenInfo->pszFilename+nNameLen-4,"hdr")))
1119 : {
1120 : CPLError( CE_Failure, CPLE_OpenFailed,
1121 : "Apparent attempt to open SIRC Convair PolGASP data failed \n"
1122 0 : "as one of the expected files is missing (hdr or img)!" );
1123 : }
1124 11981 : return NULL;
1125 : }
1126 :
1127 : /* -------------------------------------------------------------------- */
1128 : /* Confirm the requested access is supported. */
1129 : /* -------------------------------------------------------------------- */
1130 1 : if( poOpenInfo->eAccess == GA_Update )
1131 : {
1132 : CPLError( CE_Failure, CPLE_NotSupported,
1133 : "The CPG driver does not support update access to existing"
1134 0 : " datasets.\n" );
1135 0 : return NULL;
1136 : }
1137 :
1138 : /* Read the header info and create the dataset */
1139 : CPGDataset *poDS;
1140 :
1141 1 : if ( CPGType < 3 )
1142 1 : poDS = (CPGDataset *) InitializeType1Or2Dataset( poOpenInfo->pszFilename );
1143 : else
1144 0 : poDS = (CPGDataset *) InitializeType3Dataset( poOpenInfo->pszFilename );
1145 :
1146 1 : if (poDS == NULL)
1147 0 : return NULL;
1148 : /* -------------------------------------------------------------------- */
1149 : /* Check for overviews. */
1150 : /* -------------------------------------------------------------------- */
1151 : // Need to think about this.
1152 : // poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
1153 :
1154 : /* -------------------------------------------------------------------- */
1155 : /* Initialize any PAM information. */
1156 : /* -------------------------------------------------------------------- */
1157 1 : poDS->SetDescription( poOpenInfo->pszFilename );
1158 1 : poDS->TryLoadXML();
1159 :
1160 1 : return( poDS );
1161 : }
1162 :
1163 : /************************************************************************/
1164 : /* GetGCPCount() */
1165 : /************************************************************************/
1166 :
1167 0 : int CPGDataset::GetGCPCount()
1168 :
1169 : {
1170 0 : return nGCPCount;
1171 : }
1172 :
1173 : /************************************************************************/
1174 : /* GetGCPProjection() */
1175 : /************************************************************************/
1176 :
1177 0 : const char *CPGDataset::GetGCPProjection()
1178 :
1179 : {
1180 0 : return pszGCPProjection;
1181 : }
1182 :
1183 : /************************************************************************/
1184 : /* GetGCPs() */
1185 : /************************************************************************/
1186 :
1187 0 : const GDAL_GCP *CPGDataset::GetGCPs()
1188 :
1189 : {
1190 0 : return pasGCPList;
1191 : }
1192 :
1193 : /************************************************************************/
1194 : /* GetProjectionRef() */
1195 : /************************************************************************/
1196 :
1197 0 : const char *CPGDataset::GetProjectionRef()
1198 :
1199 : {
1200 0 : return( pszProjection );
1201 : }
1202 :
1203 : /************************************************************************/
1204 : /* GetGeoTransform() */
1205 : /************************************************************************/
1206 :
1207 0 : CPLErr CPGDataset::GetGeoTransform( double * padfTransform )
1208 :
1209 : {
1210 0 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
1211 0 : return( CE_None );
1212 : }
1213 :
1214 : /************************************************************************/
1215 : /* SIRC_QSLCRasterBand() */
1216 : /************************************************************************/
1217 :
1218 4 : SIRC_QSLCRasterBand::SIRC_QSLCRasterBand( CPGDataset *poGDS, int nBand,
1219 4 : GDALDataType eType )
1220 :
1221 : {
1222 4 : this->poDS = poGDS;
1223 4 : this->nBand = nBand;
1224 :
1225 4 : eDataType = eType;
1226 :
1227 4 : nBlockXSize = poGDS->nRasterXSize;
1228 4 : nBlockYSize = 1;
1229 :
1230 4 : if( nBand == 1 )
1231 1 : SetMetadataItem( "POLARIMETRIC_INTERP", "HH" );
1232 3 : else if( nBand == 2 )
1233 1 : SetMetadataItem( "POLARIMETRIC_INTERP", "HV" );
1234 2 : else if( nBand == 3 )
1235 1 : SetMetadataItem( "POLARIMETRIC_INTERP", "VH" );
1236 1 : else if( nBand == 4 )
1237 1 : SetMetadataItem( "POLARIMETRIC_INTERP", "VV" );
1238 4 : }
1239 :
1240 : /************************************************************************/
1241 : /* IReadBlock() */
1242 : /************************************************************************/
1243 :
1244 : /* From: http://southport.jpl.nasa.gov/software/dcomp/dcomp.html
1245 :
1246 : ysca = sqrt{ [ (Byte(2) / 254 ) + 1.5] 2Byte(1) }
1247 :
1248 : Re(SHH) = byte(3) ysca/127
1249 :
1250 : Im(SHH) = byte(4) ysca/127
1251 :
1252 : Re(SHV) = byte(5) ysca/127
1253 :
1254 : Im(SHV) = byte(6) ysca/127
1255 :
1256 : Re(SVH) = byte(7) ysca/127
1257 :
1258 : Im(SVH) = byte(8) ysca/127
1259 :
1260 : Re(SVV) = byte(9) ysca/127
1261 :
1262 : Im(SVV) = byte(10) ysca/127
1263 :
1264 : */
1265 :
1266 1 : CPLErr SIRC_QSLCRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
1267 : void * pImage )
1268 :
1269 : {
1270 1 : int offset, nBytesPerSample=10;
1271 : GByte *pabyRecord;
1272 1 : CPGDataset *poGDS = (CPGDataset *) poDS;
1273 : static float afPowTable[256];
1274 : static int bPowTableInitialized = FALSE;
1275 :
1276 1 : offset = nBlockXSize* nBlockYOff*nBytesPerSample;
1277 :
1278 : /* -------------------------------------------------------------------- */
1279 : /* Load all the pixel data associated with this scanline. */
1280 : /* -------------------------------------------------------------------- */
1281 1 : int nBytesToRead = nBytesPerSample * nBlockXSize;
1282 :
1283 1 : pabyRecord = (GByte *) CPLMalloc( nBytesToRead );
1284 :
1285 1 : if( VSIFSeek( poGDS->afpImage[0], offset, SEEK_SET ) != 0
1286 : || (int) VSIFRead( pabyRecord, 1, nBytesToRead,
1287 : poGDS->afpImage[0] ) != nBytesToRead )
1288 : {
1289 : CPLError( CE_Failure, CPLE_FileIO,
1290 : "Error reading %d bytes of SIRC Convair at offset %d.\n"
1291 : "Reading file %s failed.",
1292 0 : nBytesToRead, offset, poGDS->GetDescription() );
1293 0 : CPLFree( pabyRecord );
1294 0 : return CE_Failure;
1295 : }
1296 :
1297 : /* -------------------------------------------------------------------- */
1298 : /* Initialize our power table if this is our first time through. */
1299 : /* -------------------------------------------------------------------- */
1300 1 : if( !bPowTableInitialized )
1301 : {
1302 : int i;
1303 :
1304 1 : bPowTableInitialized = TRUE;
1305 :
1306 257 : for( i = 0; i < 256; i++ )
1307 : {
1308 256 : afPowTable[i] = (float) pow( 2.0, i-128 );
1309 : }
1310 : }
1311 :
1312 : /* -------------------------------------------------------------------- */
1313 : /* Copy the desired band out based on the size of the type, and */
1314 : /* the interleaving mode. */
1315 : /* -------------------------------------------------------------------- */
1316 : int iX;
1317 :
1318 2 : for( iX = 0; iX < nBlockXSize; iX++ )
1319 : {
1320 1 : unsigned char *pabyGroup = pabyRecord + iX * nBytesPerSample;
1321 1 : signed char *Byte = (signed char*)pabyGroup-1; /* A ones based alias */
1322 : double dfReSHH, dfImSHH, dfReSHV, dfImSHV,
1323 : dfReSVH, dfImSVH, dfReSVV, dfImSVV, dfScale;
1324 :
1325 1 : dfScale = sqrt( (Byte[2] / 254 + 1.5) * afPowTable[Byte[1] + 128] );
1326 :
1327 1 : if( nBand == 1 )
1328 : {
1329 1 : dfReSHH = Byte[3] * dfScale / 127.0;
1330 1 : dfImSHH = Byte[4] * dfScale / 127.0;
1331 :
1332 1 : ((float *) pImage)[iX*2 ] = (float) dfReSHH;
1333 1 : ((float *) pImage)[iX*2+1] = (float) dfImSHH;
1334 : }
1335 0 : else if( nBand == 2 )
1336 : {
1337 0 : dfReSHV = Byte[5] * dfScale / 127.0;
1338 0 : dfImSHV = Byte[6] * dfScale / 127.0;
1339 :
1340 0 : ((float *) pImage)[iX*2 ] = (float) dfReSHV;
1341 0 : ((float *) pImage)[iX*2+1] = (float) dfImSHV;
1342 : }
1343 0 : else if( nBand == 3 )
1344 : {
1345 0 : dfReSVH = Byte[7] * dfScale / 127.0;
1346 0 : dfImSVH = Byte[8] * dfScale / 127.0;
1347 :
1348 0 : ((float *) pImage)[iX*2 ] = (float) dfReSVH;
1349 0 : ((float *) pImage)[iX*2+1] = (float) dfImSVH;
1350 : }
1351 0 : else if( nBand == 4 )
1352 : {
1353 0 : dfReSVV = Byte[9] * dfScale / 127.0;
1354 0 : dfImSVV = Byte[10]* dfScale / 127.0;
1355 :
1356 0 : ((float *) pImage)[iX*2 ] = (float) dfReSVV;
1357 0 : ((float *) pImage)[iX*2+1] = (float) dfImSVV;
1358 : }
1359 : }
1360 :
1361 1 : CPLFree( pabyRecord );
1362 :
1363 1 : return CE_None;
1364 : }
1365 :
1366 : /************************************************************************/
1367 : /* CPG_STOKESRasterBand() */
1368 : /************************************************************************/
1369 :
1370 0 : CPG_STOKESRasterBand::CPG_STOKESRasterBand( GDALDataset *poDS, int nBand,
1371 : GDALDataType eType,
1372 0 : int bNativeOrder )
1373 :
1374 : {
1375 : static const char *apszPolarizations[16] = { "Covariance_11",
1376 : "Covariance_12",
1377 : "Covariance_13",
1378 : "Covariance_14",
1379 : "Covariance_21",
1380 : "Covariance_22",
1381 : "Covariance_23",
1382 : "Covariance_24",
1383 : "Covariance_31",
1384 : "Covariance_32",
1385 : "Covariance_33",
1386 : "Covariance_34",
1387 : "Covariance_41",
1388 : "Covariance_42",
1389 : "Covariance_43",
1390 : "Covariance_44" };
1391 :
1392 0 : this->poDS = poDS;
1393 0 : this->nBand = nBand;
1394 0 : this->eDataType = eType;
1395 0 : this->bNativeOrder = bNativeOrder;
1396 :
1397 0 : nBlockXSize = poDS->GetRasterXSize();
1398 0 : nBlockYSize = 1;
1399 :
1400 0 : SetMetadataItem( "POLARIMETRIC_INTERP",apszPolarizations[nBand-1] );
1401 0 : SetDescription( apszPolarizations[nBand-1] );
1402 0 : }
1403 :
1404 : /************************************************************************/
1405 : /* ~CPG_STOKESRasterBand() */
1406 : /************************************************************************/
1407 :
1408 0 : CPG_STOKESRasterBand::~CPG_STOKESRasterBand()
1409 :
1410 : {
1411 0 : }
1412 :
1413 : /************************************************************************/
1414 : /* IReadBlock() */
1415 : /************************************************************************/
1416 :
1417 : /* Convert from Stokes to Covariance representation */
1418 :
1419 0 : CPLErr CPG_STOKESRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
1420 : void * pImage )
1421 :
1422 : {
1423 : int iPixel;
1424 : int m11, m12, m13, m14, m21, m22, m23, m24, step;
1425 : int m31, m32, m33, m34, m41, m42, m43, m44;
1426 0 : CPGDataset *poGDS = (CPGDataset *) poDS;
1427 : float *M;
1428 : float *pafLine;
1429 : CPLErr eErr;
1430 :
1431 0 : CPLAssert( nBlockXOff == 0 );
1432 :
1433 0 : eErr = poGDS->LoadStokesLine(nBlockYOff, bNativeOrder);
1434 0 : if( eErr != CE_None )
1435 0 : return eErr;
1436 :
1437 0 : M = poGDS->padfStokesMatrix;
1438 0 : pafLine = ( float * ) pImage;
1439 :
1440 0 : if ( poGDS->nInterleave == BIP)
1441 : {
1442 0 : step = 16;
1443 0 : m11 = M11;
1444 0 : m12 = M12;
1445 0 : m13 = M13;
1446 0 : m14 = M14;
1447 0 : m21 = M21;
1448 0 : m22 = M22;
1449 0 : m23 = M23;
1450 0 : m24 = M24;
1451 0 : m31 = M31;
1452 0 : m32 = M32;
1453 0 : m33 = M33;
1454 0 : m34 = M34;
1455 0 : m41 = M41;
1456 0 : m42 = M42;
1457 0 : m43 = M43;
1458 0 : m44 = M44;
1459 : }
1460 : else
1461 : {
1462 0 : step = 1;
1463 0 : m11=0;
1464 0 : m12=nRasterXSize;
1465 0 : m13=nRasterXSize*2;
1466 0 : m14=nRasterXSize*3;
1467 0 : m21=nRasterXSize*4;
1468 0 : m22=nRasterXSize*5;
1469 0 : m23=nRasterXSize*6;
1470 0 : m24=nRasterXSize*7;
1471 0 : m31=nRasterXSize*8;
1472 0 : m32=nRasterXSize*9;
1473 0 : m33=nRasterXSize*10;
1474 0 : m34=nRasterXSize*11;
1475 0 : m41=nRasterXSize*12;
1476 0 : m42=nRasterXSize*13;
1477 0 : m43=nRasterXSize*14;
1478 0 : m44=nRasterXSize*15;
1479 : }
1480 0 : if ( nBand == 1 ) /* C11 */
1481 : {
1482 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1483 : {
1484 0 : pafLine[iPixel*2+0] = M[m11]-M[m22]-M[m33]+M[m44];
1485 0 : pafLine[iPixel*2+1] = 0.0;
1486 0 : m11 += step;
1487 0 : m22 += step;
1488 0 : m33 += step;
1489 0 : m44 += step;
1490 : }
1491 : }
1492 0 : else if ( nBand == 2 ) /* C12 */
1493 : {
1494 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1495 : {
1496 0 : pafLine[iPixel*2+0] = M[m13]-M[m23];
1497 0 : pafLine[iPixel*2+1] = M[m14]-M[m24];
1498 0 : m13 += step;
1499 0 : m23 += step;
1500 0 : m14 += step;
1501 0 : m24 += step;
1502 : }
1503 : }
1504 0 : else if ( nBand == 3 ) /* C13 */
1505 : {
1506 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1507 : {
1508 0 : pafLine[iPixel*2+0] = M[m33]-M[m44];
1509 0 : pafLine[iPixel*2+1] = M[m43]+M[m34];
1510 0 : m33 += step;
1511 0 : m44 += step;
1512 0 : m43 += step;
1513 0 : m34 += step;
1514 : }
1515 : }
1516 0 : else if ( nBand == 4 ) /* C14 */
1517 : {
1518 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1519 : {
1520 0 : pafLine[iPixel*2+0] = M[m31]-M[m32];
1521 0 : pafLine[iPixel*2+1] = M[m41]-M[m42];
1522 0 : m31 += step;
1523 0 : m32 += step;
1524 0 : m41 += step;
1525 0 : m42 += step;
1526 : }
1527 : }
1528 0 : else if ( nBand == 5 ) /* C21 */
1529 : {
1530 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1531 : {
1532 0 : pafLine[iPixel*2+0] = M[m13]-M[m23];
1533 0 : pafLine[iPixel*2+1] = M[m24]-M[m14];
1534 0 : m13 += step;
1535 0 : m23 += step;
1536 0 : m14 += step;
1537 0 : m24 += step;
1538 : }
1539 : }
1540 0 : else if ( nBand == 6 ) /* C22 */
1541 : {
1542 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1543 : {
1544 0 : pafLine[iPixel*2+0] = M[m11]+M[m22]-M[m33]-M[m44];
1545 0 : pafLine[iPixel*2+1] = 0.0;
1546 0 : m11 += step;
1547 0 : m22 += step;
1548 0 : m33 += step;
1549 0 : m44 += step;
1550 : }
1551 : }
1552 0 : else if ( nBand == 7 ) /* C23 */
1553 : {
1554 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1555 : {
1556 0 : pafLine[iPixel*2+0] = M[m31]+M[m32];
1557 0 : pafLine[iPixel*2+1] = M[m41]+M[m42];
1558 0 : m31 += step;
1559 0 : m32 += step;
1560 0 : m41 += step;
1561 0 : m42 += step;
1562 : }
1563 : }
1564 0 : else if ( nBand == 8 ) /* C24 */
1565 : {
1566 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1567 : {
1568 0 : pafLine[iPixel*2+0] = M[m33]+M[m44];
1569 0 : pafLine[iPixel*2+1] = M[m43]-M[m34];
1570 0 : m33 += step;
1571 0 : m44 += step;
1572 0 : m43 += step;
1573 0 : m34 += step;
1574 : }
1575 : }
1576 0 : else if ( nBand == 9 ) /* C31 */
1577 : {
1578 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1579 : {
1580 0 : pafLine[iPixel*2+0] = M[m33]-M[m44];
1581 0 : pafLine[iPixel*2+1] = -1*M[m43]-M[m34];
1582 0 : m33 += step;
1583 0 : m44 += step;
1584 0 : m43 += step;
1585 0 : m34 += step;
1586 : }
1587 : }
1588 0 : else if ( nBand == 10 ) /* C32 */
1589 : {
1590 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1591 : {
1592 0 : pafLine[iPixel*2+0] = M[m31]+M[m32];
1593 0 : pafLine[iPixel*2+1] = -1*M[m41]-M[m42];
1594 0 : m31 += step;
1595 0 : m32 += step;
1596 0 : m41 += step;
1597 0 : m42 += step;
1598 : }
1599 : }
1600 0 : else if ( nBand == 11 ) /* C33 */
1601 : {
1602 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1603 : {
1604 0 : pafLine[iPixel*2+0] = M[m11]+M[m22]+M[m33]+M[m44];
1605 0 : pafLine[iPixel*2+1] = 0.0;
1606 0 : m11 += step;
1607 0 : m22 += step;
1608 0 : m33 += step;
1609 0 : m44 += step;
1610 : }
1611 :
1612 : }
1613 0 : else if ( nBand == 12 ) /* C34 */
1614 : {
1615 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1616 : {
1617 0 : pafLine[iPixel*2+0] = M[m13]-M[m23];
1618 0 : pafLine[iPixel*2+1] = -1*M[m14]-M[m24];
1619 0 : m13 += step;
1620 0 : m23 += step;
1621 0 : m14 += step;
1622 0 : m24 += step;
1623 : }
1624 : }
1625 0 : else if ( nBand == 13 ) /* C41 */
1626 : {
1627 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1628 : {
1629 0 : pafLine[iPixel*2+0] = M[m31]-M[m32];
1630 0 : pafLine[iPixel*2+1] = M[m42]-M[m41];
1631 0 : m31 += step;
1632 0 : m32 += step;
1633 0 : m41 += step;
1634 0 : m42 += step;
1635 : }
1636 : }
1637 0 : else if ( nBand == 14 ) /* C42 */
1638 : {
1639 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1640 : {
1641 0 : pafLine[iPixel*2+0] = M[m33]+M[m44];
1642 0 : pafLine[iPixel*2+1] = M[m34]-M[m43];
1643 0 : m33 += step;
1644 0 : m44 += step;
1645 0 : m43 += step;
1646 0 : m34 += step;
1647 : }
1648 : }
1649 0 : else if ( nBand == 15 ) /* C43 */
1650 : {
1651 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1652 : {
1653 0 : pafLine[iPixel*2+0] = M[m13]-M[m23];
1654 0 : pafLine[iPixel*2+1] = M[m14]+M[m24];
1655 0 : m13 += step;
1656 0 : m23 += step;
1657 0 : m14 += step;
1658 0 : m24 += step;
1659 : }
1660 : }
1661 : else /* C44 */
1662 : {
1663 0 : for ( iPixel = 0; iPixel < nRasterXSize; iPixel++ )
1664 : {
1665 0 : pafLine[iPixel*2+0] = M[m11]-M[m22]+M[m33]-M[m44];
1666 0 : pafLine[iPixel*2+1] = 0.0;
1667 0 : m11 += step;
1668 0 : m22 += step;
1669 0 : m33 += step;
1670 0 : m44 += step;
1671 : }
1672 : }
1673 :
1674 0 : return CE_None;
1675 : }
1676 :
1677 : /************************************************************************/
1678 : /* GDALRegister_CPG() */
1679 : /************************************************************************/
1680 :
1681 582 : void GDALRegister_CPG()
1682 :
1683 : {
1684 : GDALDriver *poDriver;
1685 :
1686 582 : if( GDALGetDriverByName( "CPG" ) == NULL )
1687 : {
1688 561 : poDriver = new GDALDriver();
1689 :
1690 561 : poDriver->SetDescription( "CPG" );
1691 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1692 561 : "Convair PolGASP" );
1693 :
1694 561 : poDriver->pfnOpen = CPGDataset::Open;
1695 :
1696 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
1697 : }
1698 582 : }
1699 :
|
