1 /******************************************************************************
2 *
3 * Project: AirSAR Reader
4 * Purpose: Implements read support for AirSAR Polarimetric data.
5 * Author: Frank Warmerdam, warmerdam@pobox.com
6 *
7 ******************************************************************************
8 * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
9 * Copyright (c) 2007-2009, Even Rouault <even dot rouault at spatialys.com>
10 *
11 * Permission is hereby granted, free of charge, to any person obtaining a
12 * copy of this software and associated documentation files (the "Software"),
13 * to deal in the Software without restriction, including without limitation
14 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 * and/or sell copies of the Software, and to permit persons to whom the
16 * Software is furnished to do so, subject to the following conditions:
17 *
18 * The above copyright notice and this permission notice shall be included
19 * in all copies or substantial portions of the Software.
20 *
21 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 * DEALINGS IN THE SOFTWARE.
28 ****************************************************************************/
29
30 #include "cpl_conv.h"
31 #include "cpl_string.h"
32 #include "cpl_vsi.h"
33 #include "gdal_frmts.h"
34 #include "gdal_pam.h"
35
36 CPL_CVSID("$Id: airsardataset.cpp a5d5ed208537a05de4437e97b6a09b7ba44f76c9 2020-03-24 08:27:48 +0100 Kai Pastor $")
37
38 /************************************************************************/
39 /* ==================================================================== */
40 /* AirSARDataset */
41 /* ==================================================================== */
42 /************************************************************************/
43
44 class AirSARRasterBand;
45
46 class AirSARDataset final: public GDALPamDataset
47 {
48 friend class AirSARRasterBand;
49
50 VSILFILE *fp;
51
52 int nLoadedLine;
53 GByte *pabyCompressedLine;
54 double *padfMatrix;
55
56 int nDataStart;
57 int nRecordLength;
58
59 CPLErr LoadLine(int iLine);
60
61 static char **ReadHeader( VSILFILE * fp, int nFileOffset,
62 const char *pszPrefix, int nMaxLines );
63
64 public:
65 AirSARDataset();
66 ~AirSARDataset();
67
68 static GDALDataset *Open( GDALOpenInfo * );
69 };
70
71 /************************************************************************/
72 /* ==================================================================== */
73 /* AirSARRasterBand */
74 /* ==================================================================== */
75 /************************************************************************/
76
77 class AirSARRasterBand final: public GDALPamRasterBand
78 {
79 public:
80 AirSARRasterBand( AirSARDataset *, int );
81 ~AirSARRasterBand() override;
82
83 CPLErr IReadBlock( int, int, void * ) override;
84 };
85
86 /* locations of stokes matrix values within padfMatrix ... same order as they
87 are computed in the document. */
88
89 constexpr int M11 = 0;
90 constexpr int M12 = 1;
91 constexpr int M13 = 2;
92 constexpr int M14 = 3;
93 constexpr int M23 = 4;
94 constexpr int M24 = 5;
95 constexpr int M33 = 6;
96 constexpr int M34 = 7;
97 constexpr int M44 = 8;
98 constexpr int M22 = 9;
99
100 /************************************************************************/
101 /* AirSARRasterBand() */
102 /************************************************************************/
103
AirSARRasterBand(AirSARDataset * poDSIn,int nBandIn)104 AirSARRasterBand::AirSARRasterBand( AirSARDataset *poDSIn,
105 int nBandIn )
106
107 {
108 poDS = poDSIn;
109 nBand = nBandIn;
110
111 nBlockXSize = poDS->GetRasterXSize();
112 nBlockYSize = 1;
113
114 if( this->nBand == 2 || this->nBand == 3 || this->nBand == 5 )
115 eDataType = GDT_CFloat32;
116 else
117 eDataType = GDT_Float32;
118
119 switch( nBand )
120 {
121 case 1:
122 SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_11" );
123 SetDescription( "Covariance_11" );
124 eDataType = GDT_CFloat32;
125 break;
126
127 case 2:
128 SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_12" );
129 SetDescription( "Covariance_12" );
130 eDataType = GDT_CFloat32;
131 break;
132
133 case 3:
134 SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_13" );
135 SetDescription( "Covariance_13" );
136 eDataType = GDT_CFloat32;
137 break;
138
139 case 4:
140 SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_22" );
141 SetDescription( "Covariance_22" );
142 eDataType = GDT_CFloat32;
143 break;
144
145 case 5:
146 SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_23" );
147 SetDescription( "Covariance_23" );
148 eDataType = GDT_CFloat32;
149 break;
150
151 case 6:
152 SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_33" );
153 SetDescription( "Covariance_33" );
154 eDataType = GDT_CFloat32;
155 break;
156 }
157 }
158
159 /************************************************************************/
160 /* ~AirSARRasterBand() */
161 /************************************************************************/
162
~AirSARRasterBand()163 AirSARRasterBand::~AirSARRasterBand()
164
165 {
166 }
167
168 /************************************************************************/
169 /* IReadBlock() */
170 /************************************************************************/
171
IReadBlock(int,int nBlockYOff,void * pImage)172 CPLErr AirSARRasterBand::IReadBlock( int /* nBlockXOff */,
173 int nBlockYOff,
174 void * pImage )
175 {
176 float *pafLine = (float *) pImage;
177 const double SQRT_2 = 1.4142135623730951;
178
179 CPLErr eErr = ((AirSARDataset *)poDS)->LoadLine( nBlockYOff );
180 if( eErr != CE_None )
181 return eErr;
182
183 double *padfMatrix = ((AirSARDataset *) poDS)->padfMatrix;
184
185 if( nBand == 1 ) /* C11 */
186 {
187 for( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
188 {
189 double *m = padfMatrix + 10 * iPixel;
190
191 pafLine[iPixel*2+0] = (float)(m[M11] + m[M22] + 2 * m[M12]);
192 pafLine[iPixel*2+1] = 0.0;
193 }
194 }
195 else if( nBand == 2 ) /* C12 */
196 {
197 for( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
198 {
199 double *m = padfMatrix + 10 * iPixel;
200
201 // real
202 pafLine[iPixel*2 + 0] = (float)(SQRT_2 * (m[M13] + m[M23]));
203
204 // imaginary
205 pafLine[iPixel*2 + 1] = (float)(- SQRT_2 * (m[M24] + m[M14]));
206 }
207 }
208 else if( nBand == 3 ) /* C13 */
209 {
210 for( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
211 {
212 double *m = padfMatrix + 10 * iPixel;
213
214 // real
215 pafLine[iPixel*2 + 0] = (float)(2*m[M33] + m[M22] - m[M11]);
216
217 // imaginary
218 pafLine[iPixel*2 + 1] = (float)(-2 * m[M34]);
219 }
220 }
221 else if( nBand == 4 ) /* C22 */
222 {
223 for( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
224 {
225 double *m = padfMatrix + 10 * iPixel;
226
227 pafLine[iPixel*2+0] = (float)(2 * (m[M11] - m[M22]));
228 pafLine[iPixel*2+1] = 0.0;
229 }
230 }
231 else if( nBand == 5 ) /* C23 */
232 {
233 for( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
234 {
235 double *m = padfMatrix + 10 * iPixel;
236
237 // real
238 pafLine[iPixel*2 + 0] = (float)(SQRT_2 * (m[M13] - m[M23]));
239
240 // imaginary
241 pafLine[iPixel*2 + 1] = (float)(SQRT_2 * (m[M24] - m[M14]));
242 }
243 }
244 else if( nBand == 6 ) /* C33 */
245 {
246 for( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
247 {
248 double *m = padfMatrix + 10 * iPixel;
249
250 pafLine[iPixel*2+0] = (float)(m[M11] + m[M22] - 2 * m[M12]);
251 pafLine[iPixel*2+1] = 0.0;
252 }
253 }
254
255 return CE_None;
256 }
257
258 /************************************************************************/
259 /* ==================================================================== */
260 /* AirSARDataset */
261 /* ==================================================================== */
262 /************************************************************************/
263
264 /************************************************************************/
265 /* AirSARDataset() */
266 /************************************************************************/
267
AirSARDataset()268 AirSARDataset::AirSARDataset() :
269 fp(nullptr),
270 nLoadedLine(-1),
271 pabyCompressedLine(nullptr),
272 padfMatrix(nullptr),
273 nDataStart(0),
274 nRecordLength(0)
275 {}
276
277 /************************************************************************/
278 /* ~AirSARDataset() */
279 /************************************************************************/
280
~AirSARDataset()281 AirSARDataset::~AirSARDataset()
282
283 {
284 FlushCache();
285 CPLFree( pabyCompressedLine );
286 CPLFree( padfMatrix );
287
288 if( fp != nullptr )
289 {
290 VSIFCloseL( fp );
291 fp = nullptr;
292 }
293 }
294
295 /************************************************************************/
296 /* LoadLine() */
297 /************************************************************************/
298
LoadLine(int iLine)299 CPLErr AirSARDataset::LoadLine( int iLine )
300
301 {
302 if( iLine == nLoadedLine )
303 return CE_None;
304
305 /* -------------------------------------------------------------------- */
306 /* allocate working buffers if we don't have them already. */
307 /* -------------------------------------------------------------------- */
308 if( pabyCompressedLine == nullptr )
309 {
310 pabyCompressedLine = (GByte *) VSI_MALLOC2_VERBOSE(nRasterXSize, 10);
311
312 padfMatrix = (double *) VSI_MALLOC2_VERBOSE(10* sizeof(double), nRasterXSize);
313 if (pabyCompressedLine == nullptr ||
314 padfMatrix == nullptr)
315 {
316 CPLFree (pabyCompressedLine);
317 CPLFree (padfMatrix);
318 return CE_Failure;
319 }
320 }
321
322 CPLAssert( nRecordLength == nRasterXSize * 10 );
323
324 /* -------------------------------------------------------------------- */
325 /* Load raw compressed data. */
326 /* -------------------------------------------------------------------- */
327 if( VSIFSeekL( fp, nDataStart + iLine * nRecordLength, SEEK_SET ) != 0
328 || ((int) VSIFReadL( pabyCompressedLine, 10, nRasterXSize, fp ))
329 != nRasterXSize )
330 {
331 CPLError( CE_Failure, CPLE_FileIO,
332 "Error reading %d bytes for line %d at offset %d.\n%s",
333 nRasterXSize * 10, iLine, nDataStart + iLine * nRecordLength,
334 VSIStrerror( errno ) );
335 return CE_Failure;
336 }
337
338 /* -------------------------------------------------------------------- */
339 /* Build stokes matrix */
340 /* -------------------------------------------------------------------- */
341 for( int iPixel = 0; iPixel < nRasterXSize; iPixel++ )
342 {
343 double *M = padfMatrix + 10 * iPixel;
344 signed char *byte = (signed char *) pabyCompressedLine + 10*iPixel - 1;
345 const double gen_fac = 1.0; // should we have a general scale factor?
346
347 M[M11] = (byte[2] / 254.0 + 1.5) * pow(2.0,byte[1]) * gen_fac;
348 M[M12] = byte[3] * M[M11] / 127.0;
349 M[M13] = byte[4] * fabs((double) byte[4]) * M[M11] / (127*127);
350 M[M14] = byte[5] * fabs((double) byte[5]) * M[M11] / (127*127);
351 M[M23] = byte[6] * fabs((double) byte[6]) * M[M11] / (127*127);
352 M[M24] = byte[7] * fabs((double) byte[7]) * M[M11] / (127*127);
353 M[M33] = byte[8] * M[M11] / 127;
354 M[M34] = byte[9] * M[M11] / 127;
355 M[M44] = byte[10] * M[M11] / 127;
356 M[M22] = M[M11] - M[M33] - M[M44];
357 }
358
359 return CE_None;
360 }
361
362 /************************************************************************/
363 /* ReadHeader() */
364 /* */
365 /* Read the AirSAR header. We assume an equal sign separates */
366 /* the keyword name from the value. If not, assume the last */
367 /* "blank delimited" word is the value and everything else is a */
368 /* keyword. */
369 /* */
370 /* The records are 50 characters each. Read till we get an all */
371 /* blank record or some zero bytes. */
372 /************************************************************************/
373
ReadHeader(VSILFILE * fp,int nFileOffset,const char * pszPrefix,int nMaxLines)374 char ** AirSARDataset::ReadHeader( VSILFILE * fp, int nFileOffset,
375 const char *pszPrefix, int nMaxLines )
376
377 {
378 char **papszHeadInfo = nullptr;
379 char szLine[51];
380
381 VSIFSeekL( fp, nFileOffset, SEEK_SET );
382
383 /* ==================================================================== */
384 /* Loop collecting one line at a time. */
385 /* ==================================================================== */
386 for( int iLine = 0; iLine < nMaxLines; iLine++ )
387 {
388 /* -------------------------------------------------------------------- */
389 /* Read a 50 byte header record. */
390 /* -------------------------------------------------------------------- */
391 if( VSIFReadL( szLine, 1, 50, fp ) != 50 )
392 {
393 CPLError( CE_Failure, CPLE_FileIO,
394 "Read error collecting AirSAR header." );
395 CSLDestroy( papszHeadInfo );
396 return nullptr;
397 }
398
399 szLine[50] = '\0';
400
401 /* -------------------------------------------------------------------- */
402 /* Is it all spaces, or does it have a zero byte? */
403 /* -------------------------------------------------------------------- */
404 bool bAllSpaces = true;
405 bool bHasIllegalChars = false;
406
407 for( int i = 0; i < 50; i++ )
408 {
409 if( szLine[i] == '\0' )
410 break;
411
412 if( szLine[i] != ' ' )
413 bAllSpaces = false;
414
415 if( ((unsigned char *) szLine)[i] > 127
416 || ((unsigned char *) szLine)[i] < 10 )
417 bHasIllegalChars = true;
418 }
419
420 if( bAllSpaces || bHasIllegalChars )
421 break;
422
423 /* -------------------------------------------------------------------- */
424 /* Find the pivot between the keyword name and value. */
425 /* -------------------------------------------------------------------- */
426 int iPivot = -1;
427
428 for( int i = 0; i < 50; i++ )
429 {
430 if( szLine[i] == '=' )
431 {
432 iPivot = i;
433 break;
434 }
435 }
436
437 // If no "=" found, split on first double white space
438 if( iPivot == -1 )
439 {
440 for( int i = 48; i >= 0; i-- )
441 {
442 if( szLine[i] == ' ' && szLine[i+1] == ' ' )
443 {
444 iPivot = i;
445 break;
446 }
447 }
448 }
449
450 if( iPivot == -1 ) // Yikes!
451 {
452 CPLDebug( "AIRSAR", "No pivot in line `%s'.",
453 szLine );
454 CPLAssert( iPivot != -1 );
455 break;
456 }
457
458 /* -------------------------------------------------------------------- */
459 /* Trace ahead to the first non-white space value character. */
460 /* -------------------------------------------------------------------- */
461 int iValue = iPivot + 1;
462
463 while( iValue < 50 && szLine[iValue] == ' ' )
464 iValue++;
465
466 /* -------------------------------------------------------------------- */
467 /* Strip any white space off the keyword. */
468 /* -------------------------------------------------------------------- */
469 int iKeyEnd = iPivot - 1;
470
471 while( iKeyEnd > 0 && szLine[iKeyEnd] == ' ' )
472 iKeyEnd--;
473
474 szLine[iKeyEnd+1] = '\0';
475
476 /* -------------------------------------------------------------------- */
477 /* Convert spaces or colons into underscores in the key name. */
478 /* -------------------------------------------------------------------- */
479 for( int i = 0; szLine[i] != '\0'; i++ )
480 {
481 if( szLine[i] == ' ' || szLine[i] == ':' || szLine[i] == ',' )
482 szLine[i] = '_';
483 }
484
485 /* -------------------------------------------------------------------- */
486 /* Prefix key name with provided prefix string. */
487 /* -------------------------------------------------------------------- */
488 char szPrefixedKeyName[55];
489
490 snprintf( szPrefixedKeyName, sizeof(szPrefixedKeyName), "%s_%s", pszPrefix, szLine );
491
492 papszHeadInfo =
493 CSLSetNameValue( papszHeadInfo, szPrefixedKeyName, szLine+iValue );
494 }
495
496 return papszHeadInfo;
497 }
498
499 /************************************************************************/
500 /* Open() */
501 /************************************************************************/
502
Open(GDALOpenInfo * poOpenInfo)503 GDALDataset *AirSARDataset::Open( GDALOpenInfo * poOpenInfo )
504
505 {
506 if( poOpenInfo->fpL == nullptr || poOpenInfo->nHeaderBytes < 800 )
507 return nullptr;
508
509 /* -------------------------------------------------------------------- */
510 /* Check for AirSAR/ keyword. */
511 /* -------------------------------------------------------------------- */
512 if( !STARTS_WITH_CI((char *) poOpenInfo->pabyHeader, "RECORD LENGTH IN BYTES") )
513 return nullptr;
514
515 if( strstr((char *) poOpenInfo->pabyHeader, "COMPRESSED") == nullptr
516 || strstr((char *) poOpenInfo->pabyHeader, "JPL AIRCRAFT") == nullptr )
517 return nullptr;
518
519 /* -------------------------------------------------------------------- */
520 /* Parse the header fields. We turn all the transform the */
521 /* keywords by converting spaces to underscores so they will be */
522 /* "well behaved" as metadata keywords. */
523 /* -------------------------------------------------------------------- */
524 char **papszMD = ReadHeader( poOpenInfo->fpL, 0, "MH", 20 );
525
526 if( papszMD == nullptr )
527 return nullptr;
528
529 /* -------------------------------------------------------------------- */
530 /* Confirm the requested access is supported. */
531 /* -------------------------------------------------------------------- */
532 if( poOpenInfo->eAccess == GA_Update )
533 {
534 CPLError( CE_Failure, CPLE_NotSupported,
535 "The AIRSAR driver does not support update access to existing"
536 " datasets.\n" );
537 CSLDestroy( papszMD );
538 return nullptr;
539 }
540 /* -------------------------------------------------------------------- */
541 /* Create a corresponding GDALDataset. */
542 /* -------------------------------------------------------------------- */
543 AirSARDataset *poDS = new AirSARDataset();
544
545 /* -------------------------------------------------------------------- */
546 /* Extract some key information. */
547 /* -------------------------------------------------------------------- */
548
549 poDS->nRasterXSize =
550 atoi(CSLFetchNameValue(papszMD,"MH_NUMBER_OF_SAMPLES_PER_RECORD"));
551 poDS->nRasterYSize =
552 atoi(CSLFetchNameValue(papszMD,"MH_NUMBER_OF_LINES_IN_IMAGE"));
553
554 poDS->nRecordLength = atoi(
555 CSLFetchNameValue( papszMD, "MH_RECORD_LENGTH_IN_BYTES" ) );
556
557 poDS->nDataStart = atoi(
558 CSLFetchNameValue( papszMD, "MH_BYTE_OFFSET_OF_FIRST_DATA_RECORD" ));
559
560 /* -------------------------------------------------------------------- */
561 /* Adopt the openinfo file pointer. */
562 /* -------------------------------------------------------------------- */
563 poDS->fp = poOpenInfo->fpL;
564 poOpenInfo->fpL = nullptr;
565
566 /* -------------------------------------------------------------------- */
567 /* Read and merge parameter header into metadata. Prefix */
568 /* parameter header values with PH_. */
569 /* -------------------------------------------------------------------- */
570 int nPHOffset = 0;
571
572 if( CSLFetchNameValue( papszMD,
573 "MH_BYTE_OFFSET_OF_PARAMETER_HEADER" ) != nullptr )
574 {
575 nPHOffset = atoi(CSLFetchNameValue(
576 papszMD, "MH_BYTE_OFFSET_OF_PARAMETER_HEADER"));
577 char **papszPHInfo = ReadHeader( poDS->fp, nPHOffset, "PH", 100 );
578
579 papszMD = CSLInsertStrings( papszMD, CSLCount(papszMD), papszPHInfo );
580
581 CSLDestroy( papszPHInfo );
582 }
583
584 /* -------------------------------------------------------------------- */
585 /* Read and merge calibration header into metadata. Prefix */
586 /* parameter header values with CH_. */
587 /* -------------------------------------------------------------------- */
588 if( nPHOffset != 0 )
589 {
590 char **papszCHInfo = ReadHeader( poDS->fp,
591 nPHOffset+poDS->nRecordLength,
592 "CH", 18 );
593
594 papszMD = CSLInsertStrings( papszMD, CSLCount(papszMD), papszCHInfo );
595
596 CSLDestroy( papszCHInfo );
597 }
598
599 /* -------------------------------------------------------------------- */
600 /* Assign metadata to dataset. */
601 /* -------------------------------------------------------------------- */
602 poDS->SetMetadata( papszMD );
603 CSLDestroy( papszMD );
604
605 /* -------------------------------------------------------------------- */
606 /* Create band information objects. */
607 /* -------------------------------------------------------------------- */
608 poDS->SetBand( 1, new AirSARRasterBand( poDS, 1 ));
609 poDS->SetBand( 2, new AirSARRasterBand( poDS, 2 ));
610 poDS->SetBand( 3, new AirSARRasterBand( poDS, 3 ));
611 poDS->SetBand( 4, new AirSARRasterBand( poDS, 4 ));
612 poDS->SetBand( 5, new AirSARRasterBand( poDS, 5 ));
613 poDS->SetBand( 6, new AirSARRasterBand( poDS, 6 ));
614
615 poDS->SetMetadataItem( "MATRIX_REPRESENTATION", "SYMMETRIZED_COVARIANCE" );
616
617 /* -------------------------------------------------------------------- */
618 /* Initialize any PAM information. */
619 /* -------------------------------------------------------------------- */
620 poDS->SetDescription( poOpenInfo->pszFilename );
621 poDS->TryLoadXML();
622
623 poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
624
625 return poDS;
626 }
627
628 /************************************************************************/
629 /* GDALRegister_AirSAR() */
630 /************************************************************************/
631
GDALRegister_AirSAR()632 void GDALRegister_AirSAR()
633
634 {
635 if( GDALGetDriverByName( "AirSAR" ) != nullptr )
636 return;
637
638 GDALDriver *poDriver = new GDALDriver();
639
640 poDriver->SetDescription( "AirSAR" );
641 poDriver->SetMetadataItem( GDAL_DCAP_RASTER, "YES" );
642 poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
643 "AirSAR Polarimetric Image" );
644 poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "drivers/raster/airsar.html" );
645
646 poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
647
648 poDriver->pfnOpen = AirSARDataset::Open;
649
650 GetGDALDriverManager()->RegisterDriver( poDriver );
651 }
652