diff --git a/Utilities/GDAL/frmts/aaigrid/GNUmakefile b/Utilities/GDAL/frmts/aaigrid/GNUmakefile new file mode 100644 index 0000000000000000000000000000000000000000..567ca4b8a746d33471a0428a2fc7ea8b37df06c8 --- /dev/null +++ b/Utilities/GDAL/frmts/aaigrid/GNUmakefile @@ -0,0 +1,15 @@ + +include ../../GDALmake.opt + +OBJ = aaigriddataset.o + +CPPFLAGS := $(GDAL_INCLUDE) $(CPPFLAGS) $(XTRA_OPT) + +default: $(OBJ) + +clean: + rm -f *.o + +all: $(OBJ) + +install-obj: $(O_OBJ) diff --git a/Utilities/GDAL/frmts/aaigrid/aaigriddataset.cpp b/Utilities/GDAL/frmts/aaigrid/aaigriddataset.cpp new file mode 100644 index 0000000000000000000000000000000000000000..fe76f6b5272591b96eebea861b65ca3c2240e51b --- /dev/null +++ b/Utilities/GDAL/frmts/aaigrid/aaigriddataset.cpp @@ -0,0 +1,936 @@ +/****************************************************************************** + * $Id: aaigriddataset.cpp,v 1.34 2006/03/03 04:01:52 fwarmerdam Exp $ + * + * Project: GDAL + * Purpose: Implements Arc/Info ASCII Grid Format. + * Author: Frank Warmerdam, warmerdam@pobox.com + * + ****************************************************************************** + * Copyright (c) 2001, Frank Warmerdam (warmerdam@pobox.com) + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ***************************************************************************** + * + * $Log: aaigriddataset.cpp,v $ + * Revision 1.34 2006/03/03 04:01:52 fwarmerdam + * Added support for DX and DY instead of CELLSIZE as apparently is produced + * by Golden Surfer ascii grids. + * http://bugs.mapwindow.org/show_bug.cgi?id=88 + * + * Revision 1.33 2005/12/01 04:50:12 fwarmerdam + * Fixed dependency on poOpenInfo->fp. + * + * Revision 1.32 2005/12/01 04:26:15 fwarmerdam + * Use large file API. + * + * Revision 1.31 2005/08/23 16:32:23 fwarmerdam + * Fixed to support tabs for white space too. + * + * Revision 1.30 2005/05/05 14:01:36 fwarmerdam + * PAM Enable + * + * Revision 1.29 2005/04/15 18:34:35 fwarmerdam + * Added AREA_OR_POINT metadata. + * + * Revision 1.28 2005/02/08 18:15:34 fwarmerdam + * Removed tail recursion in IReadBlock() to avoid stack overflows. + * + * Revision 1.27 2005/02/08 17:06:27 fwarmerdam + * Fixed up Delete() method to avoid error if there is no .prj. + * + * Revision 1.26 2003/12/02 18:01:09 warmerda + * Rewrote line reading function to avoid calls to CPLReadLine() since + * some files have all the data for the whole file in one long line. + * + * Revision 1.25 2003/12/02 17:06:11 warmerda + * Write out integers as integers. + * + * Revision 1.24 2003/07/08 21:12:07 warmerda + * avoid warnings + * + * Revision 1.23 2003/05/27 17:34:22 warmerda + * fixed problem with scanlines split over multiple text lines + * + * Revision 1.22 2003/05/02 16:06:36 dron + * Memory leak fixed. + * + * Revision 1.21 2003/04/02 19:05:03 dron + * Fixes for large file support in Windows environment. + * + * Revision 1.20 2003/03/27 19:52:58 dron + * Implemented Delete() method and large file support. + * + * Revision 1.19 2003/02/06 04:55:35 warmerda + * use default georef info if bounds missing + * + * Revision 1.18 2002/11/23 18:54:17 warmerda + * added CREATIONDATATYPES metadata for drivers + * + * Revision 1.17 2002/10/02 13:10:16 warmerda + * Fixed bug in setting of Y offset derived from yllcenter, was off 1 pixel. + * As per GRASS RT bug https://intevation.de/rt/webrt?serial_num=1332. + * + * Revision 1.16 2002/09/04 06:50:36 warmerda + * avoid static driver pointers + * + * Revision 1.15 2002/07/09 21:06:54 warmerda + * fixed free of SrcDS pszProjection in CreateCopy + * + * Revision 1.14 2002/06/15 00:07:23 aubin + * mods to enable 64bit file i/o + * + * Revision 1.13 2002/06/14 12:28:41 dron + * No data handling at writing. + * + * Revision 1.12 2002/06/13 13:48:28 dron + * Added writing of .prj files for Arc/Info 8 + * + * Revision 1.11 2002/06/12 21:12:24 warmerda + * update to metadata based driver info + * + * Revision 1.10 2002/06/11 13:13:02 dron + * Write support implemented with CreateCopy() function + * + * Revision 1.9 2002/06/04 17:37:23 dron + * Header records may follow in any order now. + * + * Revision 1.8 2001/11/20 14:52:26 warmerda + * Fix from Markus for center of pixel positioning. + * + * Revision 1.7 2001/11/11 23:50:59 warmerda + * added required class keyword to friend declarations + * + * Revision 1.6 2001/11/06 14:34:22 warmerda + * Fixed bug in YLLCENTER handling. Added case for alternate line ordering. + * + * Revision 1.5 2001/07/18 04:51:56 warmerda + * added CPL_CVSID + * + * Revision 1.4 2001/03/24 20:58:30 warmerda + * Fixed typo. + * + * Revision 1.3 2001/03/13 19:39:03 warmerda + * Added help link + * + * Revision 1.2 2001/03/12 15:35:05 warmerda + * Fixed prj file naming. + * + * Revision 1.1 2001/03/12 15:15:30 warmerda + * New + * + */ + +#include "gdal_pam.h" +#include <ctype.h> +#include "cpl_string.h" +#include "ogr_spatialref.h" + +CPL_CVSID("$Id: aaigriddataset.cpp,v 1.34 2006/03/03 04:01:52 fwarmerdam Exp $"); + +CPL_C_START +void GDALRegister_AAIGrid(void); +CPL_C_END + +/************************************************************************/ +/* ==================================================================== */ +/* AAIGDataset */ +/* ==================================================================== */ +/************************************************************************/ + +class AAIGRasterBand; + +class CPL_DLL AAIGDataset : public GDALPamDataset +{ + friend class AAIGRasterBand; + + FILE *fp; + + double adfGeoTransform[6]; + char **papszPrj; + char *pszProjection; + + int bNoDataSet; + double dfNoDataValue; + + unsigned char achReadBuf[256]; + GUIntBig nBufferOffset; + int nOffsetInBuffer; + + char Getc(); + GUIntBig Tell(); + int Seek( GUIntBig nOffset ); + + public: + AAIGDataset(); + ~AAIGDataset(); + + static GDALDataset *Open( GDALOpenInfo * ); + static CPLErr Delete( const char *pszFilename ); + static CPLErr Remove( const char *pszFilename, int bRepError ); + + virtual CPLErr GetGeoTransform( double * ); + virtual const char *GetProjectionRef(void); +}; + +/************************************************************************/ +/* ==================================================================== */ +/* AAIGRasterBand */ +/* ==================================================================== */ +/************************************************************************/ + +class AAIGRasterBand : public GDALPamRasterBand +{ + friend class AAIGDataset; + + GUIntBig *panLineOffset; + + public: + + AAIGRasterBand( AAIGDataset *, int, GDALDataType ); + virtual ~AAIGRasterBand(); + + virtual double GetNoDataValue( int * ); + virtual CPLErr SetNoDataValue( double ); + virtual CPLErr IReadBlock( int, int, void * ); +}; + +/************************************************************************/ +/* AAIGRasterBand() */ +/************************************************************************/ + +AAIGRasterBand::AAIGRasterBand( AAIGDataset *poDS, int nDataStart, + GDALDataType eTypeIn ) + +{ + this->poDS = poDS; + + nBand = 1; + eDataType = eTypeIn; + + nBlockXSize = poDS->nRasterXSize; + nBlockYSize = 1; + + panLineOffset = + (GUIntBig *) CPLCalloc( poDS->nRasterYSize, sizeof(GUIntBig) ); + panLineOffset[0] = nDataStart; +} + +/************************************************************************/ +/* ~AAIGRasterBand() */ +/************************************************************************/ + +AAIGRasterBand::~AAIGRasterBand() + +{ + CPLFree( panLineOffset ); +} + +/************************************************************************/ +/* IReadBlock() */ +/************************************************************************/ + +CPLErr AAIGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff, + void * pImage ) + +{ + AAIGDataset *poODS = (AAIGDataset *) poDS; + int iPixel; + + if( nBlockYOff < 0 || nBlockYOff > poODS->nRasterYSize - 1 + || nBlockXOff != 0 ) + return CE_Failure; + + if( panLineOffset[nBlockYOff] == 0 ) + { + int iPrevLine; + + for( iPrevLine = 1; iPrevLine <= nBlockYOff; iPrevLine++ ) + if( panLineOffset[iPrevLine] == 0 ) + IReadBlock( nBlockXOff, iPrevLine-1, NULL ); + } + + if( panLineOffset[nBlockYOff] == 0 ) + return CE_Failure; + + + if( poODS->Seek( panLineOffset[nBlockYOff] ) != 0 ) + { + CPLError( CE_Failure, CPLE_FileIO, + "Can't seek to offset %ld in input file to read data.", + panLineOffset[nBlockYOff] ); + return CE_Failure; + } + + for( iPixel = 0; iPixel < poODS->nRasterXSize; ) + { + char szToken[500]; + char chNext; + int iTokenChar = 0; + + /* suck up any pre-white space. */ + do { + chNext = poODS->Getc(); + } while( isspace( chNext ) ); + + while( !isspace(chNext) ) + { + if( iTokenChar == sizeof(szToken)-2 ) + { + CPLError( CE_Failure, CPLE_FileIO, + "Token too long at scanline %d.", + nBlockYOff ); + return CE_Failure; + } + + szToken[iTokenChar++] = chNext; + chNext = poODS->Getc(); + } + + if( chNext == '\0' ) + { + CPLError( CE_Failure, CPLE_FileIO, + "File short, can't read line %d.", + nBlockYOff ); + return CE_Failure; + } + + szToken[iTokenChar] = '\0'; + + if( pImage != NULL ) + { + if( eDataType == GDT_Float32 ) + ((float *) pImage)[iPixel] = (float) atof(szToken); + else + ((GInt16 *) pImage)[iPixel] = (GInt16) atoi(szToken); + } + + iPixel++; + } + + if( nBlockYOff < poODS->nRasterYSize - 1 ) + panLineOffset[nBlockYOff + 1] = poODS->Tell(); + + return CE_None; +} + +/************************************************************************/ +/* GetNoDataValue() */ +/************************************************************************/ + +double AAIGRasterBand::GetNoDataValue( int * pbSuccess ) + +{ + AAIGDataset *poODS = (AAIGDataset *) poDS; + + if( pbSuccess ) + *pbSuccess = poODS->bNoDataSet; + + return poODS->dfNoDataValue; +} + +/************************************************************************/ +/* SetNoDataValue() */ +/************************************************************************/ + +CPLErr AAIGRasterBand::SetNoDataValue( double dfNoData ) + +{ + AAIGDataset *poODS = (AAIGDataset *) poDS; + + poODS->bNoDataSet = TRUE; + poODS->dfNoDataValue = dfNoData; + + return CE_None; +} + +/************************************************************************/ +/* ==================================================================== */ +/* AAIGDataset */ +/* ==================================================================== */ +/************************************************************************/ + + +/************************************************************************/ +/* AAIGDataset() */ +/************************************************************************/ + +AAIGDataset::AAIGDataset() + +{ + papszPrj = NULL; + pszProjection = CPLStrdup(""); + fp = NULL; + bNoDataSet = FALSE; + dfNoDataValue = -9999.0; + + adfGeoTransform[0] = 0.0; + adfGeoTransform[1] = 1.0; + adfGeoTransform[2] = 0.0; + adfGeoTransform[3] = 0.0; + adfGeoTransform[4] = 0.0; + adfGeoTransform[5] = 1.0; + + nOffsetInBuffer = 256; + nBufferOffset = 0; +} + +/************************************************************************/ +/* ~AAIGDataset() */ +/************************************************************************/ + +AAIGDataset::~AAIGDataset() + +{ + FlushCache(); + + if( fp != NULL ) + VSIFCloseL( fp ); + + CPLFree( pszProjection ); + CSLDestroy( papszPrj ); +} + +/************************************************************************/ +/* Tell() */ +/************************************************************************/ + +GUIntBig AAIGDataset::Tell() + +{ + return nBufferOffset + nOffsetInBuffer; +} + +/************************************************************************/ +/* Seek() */ +/************************************************************************/ + +int AAIGDataset::Seek( GUIntBig nNewOffset ) + +{ + nOffsetInBuffer = sizeof(achReadBuf); + return VSIFSeekL( fp, nNewOffset, SEEK_SET ); +} + +/************************************************************************/ +/* Getc() */ +/* */ +/* Read a single character from the input file (efficiently we */ +/* hope). */ +/************************************************************************/ + +char AAIGDataset::Getc() + +{ + if( nOffsetInBuffer < (int) sizeof(achReadBuf) ) + return achReadBuf[nOffsetInBuffer++]; + + nBufferOffset = VSIFTellL( fp ); + if( VSIFReadL( achReadBuf, 1, sizeof(achReadBuf), fp ) < 1 ) + return EOF; + + nOffsetInBuffer = 0; + + return achReadBuf[nOffsetInBuffer++]; +} + +/************************************************************************/ +/* Open() */ +/************************************************************************/ + +GDALDataset *AAIGDataset::Open( GDALOpenInfo * poOpenInfo ) + +{ + int i, j; + GDALDataType eDataType; + char **papszTokens; + +/* -------------------------------------------------------------------- */ +/* Does this look like an AI grid file? */ +/* -------------------------------------------------------------------- */ + if( poOpenInfo->nHeaderBytes < 100 + || !( EQUALN((const char *) poOpenInfo->pabyHeader,"ncols",5) || + EQUALN((const char *) poOpenInfo->pabyHeader,"nrows",5) || + EQUALN((const char *) poOpenInfo->pabyHeader,"xllcorner",9)|| + EQUALN((const char *) poOpenInfo->pabyHeader,"yllcorner",9)|| + EQUALN((const char *) poOpenInfo->pabyHeader,"xllcenter",9)|| + EQUALN((const char *) poOpenInfo->pabyHeader,"yllcenter",9)|| + EQUALN((const char *) poOpenInfo->pabyHeader,"dx",2)|| + EQUALN((const char *) poOpenInfo->pabyHeader,"dy",2)|| + EQUALN((const char *) poOpenInfo->pabyHeader,"cellsize",8)) ) + return NULL; + + papszTokens = + CSLTokenizeString2( (const char *) poOpenInfo->pabyHeader, + " \n\r\t", 0 ); + +/* -------------------------------------------------------------------- */ +/* Create a corresponding GDALDataset. */ +/* -------------------------------------------------------------------- */ + AAIGDataset *poDS; + + poDS = new AAIGDataset(); + +/* -------------------------------------------------------------------- */ +/* Parse the header. */ +/* -------------------------------------------------------------------- */ + double dfCellDX, dfCellDY; + + if ( (i = CSLFindString( papszTokens, "ncols" )) < 0 ) + { + CSLDestroy( papszTokens ); + return NULL; + } + poDS->nRasterXSize = atoi(papszTokens[i + 1]); + if ( (i = CSLFindString( papszTokens, "nrows" )) < 0 ) + { + CSLDestroy( papszTokens ); + return NULL; + } + poDS->nRasterYSize = atoi(papszTokens[i + 1]); + if ( (i = CSLFindString( papszTokens, "cellsize" )) < 0 ) + { + int iDX, iDY; + if( (iDX = CSLFindString(papszTokens,"dx")) < 0 + || (iDY = CSLFindString(papszTokens,"dy")) < 0 ) + { + CSLDestroy( papszTokens ); + return NULL; + } + + dfCellDX = atof( papszTokens[iDX+1] ); + dfCellDY = atof( papszTokens[iDY+1] ); + } + else + dfCellDX = dfCellDY = atof( papszTokens[i + 1] ); + + if ((i = CSLFindString( papszTokens, "xllcorner" )) >= 0 && + (j = CSLFindString( papszTokens, "yllcorner" )) >= 0 ) + { + poDS->adfGeoTransform[0] = atof( papszTokens[i + 1] ); + poDS->adfGeoTransform[1] = dfCellDX; + poDS->adfGeoTransform[2] = 0.0; + poDS->adfGeoTransform[3] = atof( papszTokens[j + 1] ) + + poDS->nRasterYSize * dfCellDY; + poDS->adfGeoTransform[4] = 0.0; + poDS->adfGeoTransform[5] = - dfCellDY; + } + else if ((i = CSLFindString( papszTokens, "xllcenter" )) >= 0 && + (j = CSLFindString( papszTokens, "yllcenter" )) >= 0 ) + { + poDS->SetMetadataItem( GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT ); + + poDS->adfGeoTransform[0] = atof(papszTokens[i + 1]) - 0.5 * dfCellDX; + poDS->adfGeoTransform[1] = dfCellDX; + poDS->adfGeoTransform[2] = 0.0; + poDS->adfGeoTransform[3] = atof( papszTokens[j + 1] ) + - 0.5 * dfCellDY + + poDS->nRasterYSize * dfCellDY; + poDS->adfGeoTransform[4] = 0.0; + poDS->adfGeoTransform[5] = - dfCellDY; + } + else + { + poDS->adfGeoTransform[0] = 0.0; + poDS->adfGeoTransform[1] = dfCellDX; + poDS->adfGeoTransform[2] = 0.0; + poDS->adfGeoTransform[3] = 0.0; + poDS->adfGeoTransform[4] = 0.0; + poDS->adfGeoTransform[5] = - dfCellDY; + } + + if( (i = CSLFindString( papszTokens, "NODATA_value" )) >= 0 ) + { + poDS->bNoDataSet = TRUE; + poDS->dfNoDataValue = atof(papszTokens[i + 1]); + } + + CSLDestroy( papszTokens ); + +/* -------------------------------------------------------------------- */ +/* Open file with large file API. */ +/* -------------------------------------------------------------------- */ + poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "r" ); + if( poDS->fp == NULL ) + { + CPLError( CE_Failure, CPLE_OpenFailed, + "VSIFOpenL(%s) failed unexpectedly.", + poOpenInfo->pszFilename ); + return NULL; + } + +/* -------------------------------------------------------------------- */ +/* Find the start of real data. */ +/* -------------------------------------------------------------------- */ + int nStartOfData; + + for( i = 2; TRUE ; i++ ) + { + if( poOpenInfo->pabyHeader[i] == '\0' ) + { + CPLError( CE_Failure, CPLE_AppDefined, + "Couldn't find data values in ASCII Grid file.\n" ); + return NULL; + } + + if( poOpenInfo->pabyHeader[i-1] == '\n' + || poOpenInfo->pabyHeader[i-2] == '\n' ) + { + if( !isalpha(poOpenInfo->pabyHeader[i]) ) + { + nStartOfData = i; + if( strstr((const char *)poOpenInfo->pabyHeader+i,".") != NULL) + eDataType = GDT_Float32; + else + eDataType = GDT_Int16; + + break; + } + } + } + +/* -------------------------------------------------------------------- */ +/* Create band information objects. */ +/* -------------------------------------------------------------------- */ + poDS->SetBand( 1, new AAIGRasterBand( poDS, nStartOfData, eDataType ) ); + +/* -------------------------------------------------------------------- */ +/* Try to read projection file. */ +/* -------------------------------------------------------------------- */ + char *pszDirname, *pszBasename; + const char *pszPrjFilename; + VSIStatBufL sStatBuf; + + pszDirname = CPLStrdup(CPLGetPath(poOpenInfo->pszFilename)); + pszBasename = CPLStrdup(CPLGetBasename(poOpenInfo->pszFilename)); + + pszPrjFilename = CPLFormFilename( pszDirname, pszBasename, "prj" ); + if( VSIStatL( pszPrjFilename, &sStatBuf ) == 0 ) + { + OGRSpatialReference oSRS; + + poDS->papszPrj = CSLLoad( pszPrjFilename ); + + if( oSRS.importFromESRI( poDS->papszPrj ) == OGRERR_NONE ) + { + CPLFree( poDS->pszProjection ); + oSRS.exportToWkt( &(poDS->pszProjection) ); + } + } + + CPLFree( pszDirname ); + CPLFree( pszBasename ); + +/* -------------------------------------------------------------------- */ +/* Initialize any PAM information. */ +/* -------------------------------------------------------------------- */ + poDS->SetDescription( poOpenInfo->pszFilename ); + poDS->TryLoadXML(); + + return( poDS ); +} + +/************************************************************************/ +/* GetGeoTransform() */ +/************************************************************************/ + +CPLErr AAIGDataset::GetGeoTransform( double * padfTransform ) + +{ + memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 ); + return( CE_None ); +} + +/************************************************************************/ +/* GetProjectionRef() */ +/************************************************************************/ + +const char *AAIGDataset::GetProjectionRef() + +{ + return pszProjection; +} + +/************************************************************************/ +/* AAIGCreateCopy() */ +/************************************************************************/ + +static GDALDataset * +AAIGCreateCopy( const char * pszFilename, GDALDataset *poSrcDS, + int bStrict, char ** papszOptions, + GDALProgressFunc pfnProgress, void * pProgressData ) + +{ + int nBands = poSrcDS->GetRasterCount(); + int nXSize = poSrcDS->GetRasterXSize(); + int nYSize = poSrcDS->GetRasterYSize(); + +/* -------------------------------------------------------------------- */ +/* Some rudimentary checks */ +/* -------------------------------------------------------------------- */ + if( nBands != 1 ) + { + CPLError( CE_Failure, CPLE_NotSupported, + "AAIG driver doesn't support %d bands. Must be 1 band.\n", + nBands ); + + return NULL; + } + + if( !pfnProgress( 0.0, NULL, pProgressData ) ) + return NULL; + +/* -------------------------------------------------------------------- */ +/* Create the dataset. */ +/* -------------------------------------------------------------------- */ + FILE *fpImage; + + fpImage = VSIFOpenL( pszFilename, "wt" ); + if( fpImage == NULL ) + { + CPLError( CE_Failure, CPLE_OpenFailed, + "Unable to create file %s.\n", + pszFilename ); + return NULL; + } + +/* -------------------------------------------------------------------- */ +/* Write ASCII Grid file header */ +/* -------------------------------------------------------------------- */ + double adfGeoTransform[6]; + char szHeader[2000]; + + poSrcDS->GetGeoTransform( adfGeoTransform ); + + if( ABS(adfGeoTransform[1]+adfGeoTransform[5]) < 0.0000001 ) + sprintf( szHeader, + "ncols %d\n" + "nrows %d\n" + "xllcorner %.12f\n" + "yllcorner %.12f\n" + "cellsize %.12f\n", + nXSize, nYSize, + adfGeoTransform[0], + adfGeoTransform[3]- nYSize * adfGeoTransform[1], + adfGeoTransform[1] ); + else + sprintf( szHeader, + "ncols %d\n" + "nrows %d\n" + "xllcorner %.12f\n" + "yllcorner %.12f\n" + "dx %.12f\n" + "dy %.12f\n", + nXSize, nYSize, + adfGeoTransform[0], + adfGeoTransform[3]- nYSize * adfGeoTransform[1], + adfGeoTransform[1], + fabs(adfGeoTransform[5]) ); + +/* -------------------------------------------------------------------- */ +/* Handle nodata (optionally). */ +/* -------------------------------------------------------------------- */ + GDALRasterBand * poBand = poSrcDS->GetRasterBand( 1 ); + double dfNoData; + int bSuccess; + + // Write `nodata' value to header if it is exists in source dataset + dfNoData = poBand->GetNoDataValue( &bSuccess ); + if ( bSuccess ) + sprintf( szHeader+strlen(szHeader), "NODATA_value %6.20g\n", + dfNoData ); + + VSIFWriteL( szHeader, 1, strlen(szHeader), fpImage ); + +/* -------------------------------------------------------------------- */ +/* Loop over image, copying image data. */ +/* -------------------------------------------------------------------- */ + double *padfScanline; + int iLine, iPixel; + CPLErr eErr = CE_None; + + // Write scanlines to output file + padfScanline = (double *) CPLMalloc( nXSize * + GDALGetDataTypeSize(GDT_CFloat64) / 8 ); + for( iLine = 0; eErr == CE_None && iLine < nYSize; iLine++ ) + { + eErr = poBand->RasterIO( GF_Read, 0, iLine, nXSize, 1, + padfScanline, nXSize, 1, GDT_CFloat64, + sizeof(double), sizeof(double) * nXSize ); + + if( poBand->GetRasterDataType() == GDT_Byte + || poBand->GetRasterDataType() == GDT_Int16 + || poBand->GetRasterDataType() == GDT_UInt16 + || poBand->GetRasterDataType() == GDT_Int32 ) + { + for ( iPixel = 0; iPixel < nXSize; iPixel++ ) + { + sprintf( szHeader, " %d", (int) padfScanline[iPixel] ); + if( VSIFWriteL( szHeader, strlen(szHeader), 1, fpImage ) != 1 ) + { + eErr = CE_Failure; + CPLError( CE_Failure, CPLE_AppDefined, + "Write failed, disk full?\n" ); + break; + } + } + } + else + { + for ( iPixel = 0; iPixel < nXSize; iPixel++ ) + { + sprintf( szHeader, " %6.20g", padfScanline[iPixel] ); + if( VSIFWriteL( szHeader, strlen(szHeader), 1, fpImage ) != 1 ) + { + eErr = CE_Failure; + CPLError( CE_Failure, CPLE_AppDefined, + "Write failed, disk full?\n" ); + break; + } + } + } + VSIFWriteL( (void *) "\n", 1, 1, fpImage ); + + if( eErr == CE_None && + !pfnProgress((iLine + 1) / ((double) nYSize), NULL, pProgressData) ) + { + eErr = CE_Failure; + CPLError( CE_Failure, CPLE_UserInterrupt, + "User terminated CreateCopy()" ); + } + } + + CPLFree( padfScanline ); + VSIFCloseL( fpImage ); + +/* -------------------------------------------------------------------- */ +/* Try to write projection file. */ +/* -------------------------------------------------------------------- */ + const char *pszOriginalProjection; + + pszOriginalProjection = (char *)poSrcDS->GetProjectionRef(); + if( !EQUAL( pszOriginalProjection, "" ) ) + { + char *pszDirname, *pszBasename; + const char *pszPrjFilename; + char *pszESRIProjection = NULL; + FILE *fp; + OGRSpatialReference oSRS; + + pszDirname = CPLStrdup( CPLGetPath(pszFilename) ); + pszBasename = CPLStrdup( CPLGetBasename(pszFilename) ); + + pszPrjFilename = CPLFormFilename( pszDirname, pszBasename, "prj" ); + fp = VSIFOpenL( pszPrjFilename, "wt" ); + + oSRS.importFromWkt( (char **) &pszOriginalProjection ); + oSRS.morphToESRI(); + oSRS.exportToWkt( &pszESRIProjection ); + VSIFWriteL( pszESRIProjection, 1, strlen(pszESRIProjection), fp ); + + VSIFCloseL( fp ); + CPLFree( pszDirname ); + CPLFree( pszBasename ); + CPLFree( pszESRIProjection ); + } + +/* -------------------------------------------------------------------- */ +/* Re-open dataset, and copy any auxilary pam information. */ +/* -------------------------------------------------------------------- */ + GDALPamDataset *poDS = (GDALPamDataset *) + GDALOpen( pszFilename, GA_ReadOnly ); + + if( poDS ) + poDS->CloneInfo( poSrcDS, GCIF_PAM_DEFAULT ); + + return poDS; +} + +/************************************************************************/ +/* Remove() */ +/* Called from the Delete() */ +/************************************************************************/ + +CPLErr AAIGDataset::Remove( const char * pszFilename, int bRepError ) + +{ + VSIStatBufL sStat; + + if( VSIStatL( pszFilename, &sStat ) == 0 && VSI_ISREG( sStat.st_mode ) ) + { + if( VSIUnlink( pszFilename ) == 0 ) + return CE_None; + else + { + CPLError( CE_Failure, CPLE_AppDefined, + "Attempt to unlink %s failed.\n", pszFilename ); + return CE_Failure; + } + } + else if( bRepError ) + { + + CPLError( CE_Failure, CPLE_AppDefined, + "Unable to delete %s, not a file.\n", pszFilename ); + return CE_Failure; + } + + return CE_None; +} + +/************************************************************************/ +/* Delete() */ +/************************************************************************/ + +CPLErr AAIGDataset::Delete( const char *pszFilename ) + +{ + Remove( CPLResetExtension( pszFilename, "prj" ), FALSE ); + return Remove( pszFilename, TRUE ); +} + +/************************************************************************/ +/* GDALRegister_AAIGrid() */ +/************************************************************************/ + +void GDALRegister_AAIGrid() + +{ + GDALDriver *poDriver; + + if( GDALGetDriverByName( "AAIGrid" ) == NULL ) + { + poDriver = new GDALDriver(); + + poDriver->SetDescription( "AAIGrid" ); + poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, + "Arc/Info ASCII Grid" ); + poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, + "frmt_various.html#AAIGrid" ); + poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "asc" ); + poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, + "Byte UInt16 Int16 Float32" ); + + poDriver->pfnOpen = AAIGDataset::Open; + poDriver->pfnCreateCopy = AAIGCreateCopy; + poDriver->pfnDelete = AAIGDataset::Delete; + + GetGDALDriverManager()->RegisterDriver( poDriver ); + } +} + diff --git a/Utilities/GDAL/frmts/aaigrid/makefile.vc b/Utilities/GDAL/frmts/aaigrid/makefile.vc new file mode 100644 index 0000000000000000000000000000000000000000..f16a3085739746ff473ac6f9da58f8607883b420 --- /dev/null +++ b/Utilities/GDAL/frmts/aaigrid/makefile.vc @@ -0,0 +1,13 @@ + +OBJ = aaigriddataset.obj + +GDAL_ROOT = ..\.. + +!INCLUDE $(GDAL_ROOT)\nmake.opt + +default: $(OBJ) + copy *.obj ..\o + +clean: + -del *.obj + diff --git a/Utilities/GDAL/frmts/airsar/GNUmakefile b/Utilities/GDAL/frmts/airsar/GNUmakefile new file mode 100644 index 0000000000000000000000000000000000000000..d0d08c9b51a6e26dfc8d7c1c96a0efc8c981b694 --- /dev/null +++ b/Utilities/GDAL/frmts/airsar/GNUmakefile @@ -0,0 +1,13 @@ + +include ../../GDALmake.opt + +OBJ = airsardataset.o + +CPPFLAGS := $(GDAL_INCLUDE) $(CPPFLAGS) + +default: $(OBJ) + +clean: + rm -f *.o $(O_OBJ) + +install-obj: $(O_OBJ) diff --git a/Utilities/GDAL/frmts/airsar/airsardataset.cpp b/Utilities/GDAL/frmts/airsar/airsardataset.cpp new file mode 100644 index 0000000000000000000000000000000000000000..9b0e7cc60f9223e9fafea7c0d97e4eac74d81a03 --- /dev/null +++ b/Utilities/GDAL/frmts/airsar/airsardataset.cpp @@ -0,0 +1,666 @@ +/****************************************************************************** + * $Id: airsardataset.cpp,v 1.6 2005/05/05 15:52:48 fwarmerdam Exp $ + * + * Project: AirSAR Reader + * Purpose: Implements read support for AirSAR Polarimetric data. + * Author: Frank Warmerdam, warmerdam@pobox.com + * + ****************************************************************************** + * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com> + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************** + * + * $Log: airsardataset.cpp,v $ + * Revision 1.6 2005/05/05 15:52:48 fwarmerdam + * PAM Enabled + * + * Revision 1.5 2004/11/11 00:16:01 gwalter + * Polarmetric->Polarimetric. + * + * Revision 1.4 2004/10/12 15:42:25 fwarmerdam + * Change all bands to be complex so that overviews will build properly + * + * Revision 1.3 2004/10/08 12:22:09 fwarmerdam + * Cast bytes passed to fabs(): Martin Daly + * + * Revision 1.2 2004/10/05 18:25:45 fwarmerdam + * Added help pointer. + * + * Revision 1.1 2004/10/05 18:12:26 fwarmerdam + * New + * + */ + +#include "gdal_pam.h" +#include "cpl_string.h" +#include "cpl_conv.h" +#include "cpl_vsi.h" + +CPL_CVSID("$Id: airsardataset.cpp,v 1.6 2005/05/05 15:52:48 fwarmerdam Exp $"); + +CPL_C_START +void GDALRegister_AirSAR(void); +CPL_C_END + +/************************************************************************/ +/* ==================================================================== */ +/* AirSARDataset */ +/* ==================================================================== */ +/************************************************************************/ + +class AirSARRasterBand; + +class AirSARDataset : public GDALPamDataset +{ + friend class AirSARRasterBand; + + FILE *fp; + + int nLoadedLine; + GByte *pabyCompressedLine; + double *padfMatrix; + + int nDataStart; + int nRecordLength; + + CPLErr LoadLine(int iLine); + + static char **ReadHeader( FILE * fp, int nFileOffset, + const char *pszPrefix, int nMaxLines ); + + public: + AirSARDataset(); + ~AirSARDataset(); + + static GDALDataset *Open( GDALOpenInfo * ); +}; + +/************************************************************************/ +/* ==================================================================== */ +/* AirSARRasterBand */ +/* ==================================================================== */ +/************************************************************************/ + +class AirSARRasterBand : public GDALPamRasterBand +{ + public: + AirSARRasterBand( AirSARDataset *, int ); + virtual ~AirSARRasterBand(); + + virtual CPLErr IReadBlock( int, int, void * ); +}; + +/* locations of stokes matrix values within padfMatrix ... same order as they + are computed in the document. */ + +#define M11 0 +#define M12 1 +#define M13 2 +#define M14 3 +#define M23 4 +#define M24 5 +#define M33 6 +#define M34 7 +#define M44 8 +#define M22 9 + +/************************************************************************/ +/* AirSARRasterBand() */ +/************************************************************************/ + +AirSARRasterBand::AirSARRasterBand( AirSARDataset *poDS, + int nBand ) + +{ + this->poDS = poDS; + this->nBand = nBand; + + nBlockXSize = poDS->GetRasterXSize(); + nBlockYSize = 1; + + if( this->nBand == 2 || this->nBand == 3 || this->nBand == 5 ) + eDataType = GDT_CFloat32; + else + eDataType = GDT_Float32; + + switch( nBand ) + { + case 1: + SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_11" ); + SetDescription( "Covariance_11" ); + eDataType = GDT_CFloat32; + break; + + case 2: + SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_12" ); + SetDescription( "Covariance_12" ); + eDataType = GDT_CFloat32; + break; + + case 3: + SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_13" ); + SetDescription( "Covariance_13" ); + eDataType = GDT_CFloat32; + break; + + case 4: + SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_22" ); + SetDescription( "Covariance_22" ); + eDataType = GDT_CFloat32; + break; + + case 5: + SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_23" ); + SetDescription( "Covariance_23" ); + eDataType = GDT_CFloat32; + break; + + case 6: + SetMetadataItem( "POLARIMETRIC_INTERP", "Covariance_33" ); + SetDescription( "Covariance_33" ); + eDataType = GDT_CFloat32; + break; + } +} + +/************************************************************************/ +/* ~AirSARRasterBand() */ +/************************************************************************/ + +AirSARRasterBand::~AirSARRasterBand() + +{ +} + +/************************************************************************/ +/* IReadBlock() */ +/************************************************************************/ + +CPLErr AirSARRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff, + void * pImage ) + +{ + CPLErr eErr; + float *pafLine = (float *) pImage; + int iPixel; + double *padfMatrix; + + eErr = ((AirSARDataset *)poDS)->LoadLine( nBlockYOff ); + if( eErr != CE_None ) + return eErr; + + padfMatrix = ((AirSARDataset *) poDS)->padfMatrix; + +#define SQRT_2 1.4142135623730951 + + if( nBand == 1 ) /* C11 */ + { + for( iPixel = 0; iPixel < nRasterXSize; iPixel++ ) + { + double *m = padfMatrix + 10 * iPixel; + + pafLine[iPixel*2+0] = m[M11] + m[M22] + 2 * m[M12]; + pafLine[iPixel*2+1] = 0.0; + } + } + else if( nBand == 2 ) /* C12 */ + { + for( iPixel = 0; iPixel < nRasterXSize; iPixel++ ) + { + double *m = padfMatrix + 10 * iPixel; + + // real + pafLine[iPixel*2 + 0] = SQRT_2 * (m[M13] + m[M23]); + + // imaginary + pafLine[iPixel*2 + 1] = - SQRT_2 * (m[M24] + m[M14]); + } + } + else if( nBand == 3 ) /* C13 */ + { + for( iPixel = 0; iPixel < nRasterXSize; iPixel++ ) + { + double *m = padfMatrix + 10 * iPixel; + + // real + pafLine[iPixel*2 + 0] = 2*m[M33] + m[M22] - m[M11]; + + // imaginary + pafLine[iPixel*2 + 1] = -2 * m[M34]; + } + } + else if( nBand == 4 ) /* C22 */ + { + for( iPixel = 0; iPixel < nRasterXSize; iPixel++ ) + { + double *m = padfMatrix + 10 * iPixel; + + pafLine[iPixel*2+0] = 2 * (m[M11] - m[M22]); + pafLine[iPixel*2+1] = 0.0; + } + } + else if( nBand == 5 ) /* C23 */ + { + for( iPixel = 0; iPixel < nRasterXSize; iPixel++ ) + { + double *m = padfMatrix + 10 * iPixel; + + // real + pafLine[iPixel*2 + 0] = SQRT_2 * (m[M13] - m[M23]); + + // imaginary + pafLine[iPixel*2 + 1] = SQRT_2 * (m[M23] - m[M14]); + } + } + else if( nBand == 6 ) /* C33 */ + { + for( iPixel = 0; iPixel < nRasterXSize; iPixel++ ) + { + double *m = padfMatrix + 10 * iPixel; + + pafLine[iPixel*2+0] = m[M11] + m[M22] - 2 * m[M12]; + pafLine[iPixel*2+1] = 0.0; + } + } + + return CE_None; +} + +/************************************************************************/ +/* ==================================================================== */ +/* AirSARDataset */ +/* ==================================================================== */ +/************************************************************************/ + +/************************************************************************/ +/* AirSARDataset() */ +/************************************************************************/ + +AirSARDataset::AirSARDataset() + +{ + fp = NULL; + + nLoadedLine = -1; + pabyCompressedLine = NULL; + padfMatrix = NULL; +} + +/************************************************************************/ +/* ~AirSARDataset() */ +/************************************************************************/ + +AirSARDataset::~AirSARDataset() + +{ + FlushCache(); + if( pabyCompressedLine != NULL ) + { + CPLFree( pabyCompressedLine ); + CPLFree( padfMatrix ); + } + + if( fp != NULL ) + { + VSIFClose( fp ); + fp = NULL; + } +} + +/************************************************************************/ +/* LoadLine() */ +/************************************************************************/ + +CPLErr AirSARDataset::LoadLine( int iLine ) + +{ + if( iLine == nLoadedLine ) + return CE_None; + +/* -------------------------------------------------------------------- */ +/* allocate working buffers if we don't have them already. */ +/* -------------------------------------------------------------------- */ + if( pabyCompressedLine == NULL ) + { + pabyCompressedLine = (GByte *) CPLMalloc(nRasterXSize * 10); + + padfMatrix = (double *) CPLMalloc(sizeof(double) * nRasterXSize*10); + } + + CPLAssert( nRecordLength == nRasterXSize * 10 ); + +/* -------------------------------------------------------------------- */ +/* Load raw compressed data. */ +/* -------------------------------------------------------------------- */ + if( VSIFSeek( fp, nDataStart + iLine * nRecordLength, SEEK_SET ) != 0 + || ((int) VSIFRead( pabyCompressedLine, 10, nRasterXSize, fp )) + != nRasterXSize ) + { + CPLError( CE_Failure, CPLE_FileIO, + "Error reading %d bytes for line %d at offset %d.\n%s", + nRasterXSize * 10, iLine, nDataStart + iLine * nRecordLength, + VSIStrerror( errno ) ); + return CE_Failure; + } + +/* -------------------------------------------------------------------- */ +/* Build stokes matrix */ +/* -------------------------------------------------------------------- */ + for( int iPixel = 0; iPixel < nRasterXSize; iPixel++ ) + { + double *M = padfMatrix + 10 * iPixel; + signed char *byte = (signed char *) pabyCompressedLine + 10*iPixel - 1; + double gen_fac = 1.0; // should we have a general scale factor? + + M[M11] = (byte[2] / 254.0 + 1.5) * pow(2.0,byte[1]) * gen_fac; + M[M12] = byte[3] * M[M11] / 127.0; + M[M13] = byte[4] * fabs((double) byte[4]) * M[M11] / (127*127); + M[M14] = byte[5] * fabs((double) byte[5]) * M[M11] / (127*127); + M[M23] = byte[6] * fabs((double) byte[6]) * M[M11] / (127*127); + M[M24] = byte[7] * fabs((double) byte[7]) * M[M11] / (127*127); + M[M33] = byte[8] * M[M11] / 127; + M[M34] = byte[9] * M[M11] / 127; + M[M44] = byte[10] * M[M11] / 127; + M[M22] = M[M11] - M[M33] - M[M44]; + } + + return CE_None; +} + +/************************************************************************/ +/* ReadHeader() */ +/* */ +/* Read the AirSAR header. We assume an equal sign seperates */ +/* the keyword name from the value. If not, assume the last */ +/* "blank delimited" word is the value and everything else is a */ +/* keyword. */ +/* */ +/* The records are 50 characters each. Read till we get an all */ +/* blank record or some zero bytes. */ +/************************************************************************/ + +char ** AirSARDataset::ReadHeader( FILE * fp, int nFileOffset, + const char *pszPrefix, int nMaxLines ) + +{ + char **papszHeadInfo = NULL; + char szLine[51]; + int iLine; + + VSIFSeek( fp, nFileOffset, SEEK_SET ); + +/* ==================================================================== */ +/* Loop collecting one line at a time. */ +/* ==================================================================== */ + for( iLine = 0; iLine < nMaxLines; iLine++ ) + { +/* -------------------------------------------------------------------- */ +/* Read a 50 byte header record. */ +/* -------------------------------------------------------------------- */ + if( VSIFRead( szLine, 1, 50, fp ) != 50 ) + { + CPLError( CE_Failure, CPLE_FileIO, + "Read error collecting AirSAR header." ); + return NULL; + } + + szLine[50] = '\0'; + +/* -------------------------------------------------------------------- */ +/* Is it all spaces, or does it have a zero byte? */ +/* -------------------------------------------------------------------- */ + int bAllSpaces = TRUE; + int bHasIllegalChars = FALSE; + int i; + + for( i = 0; i < 50; i++ ) + { + if( szLine[i] == '\0' ) + break; + + if( szLine[i] != ' ' ) + bAllSpaces = FALSE; + + if( ((unsigned char *) szLine)[i] > 127 + || ((unsigned char *) szLine)[i] < 10 ) + bHasIllegalChars = TRUE; + } + + if( bAllSpaces || bHasIllegalChars ) + break; + +/* -------------------------------------------------------------------- */ +/* Find the pivot between the keyword name and value. */ +/* -------------------------------------------------------------------- */ + int iPivot = -1; + + for( i = 0; i < 50; i++ ) + { + if( szLine[i] == '=' ) + { + iPivot = i; + break; + } + } + + // If no "=" found, split on first double white space + if( iPivot == -1 ) + { + for( i = 48; i >= 0; i-- ) + { + if( szLine[i] == ' ' && szLine[i+1] == ' ' ) + { + iPivot = i; + break; + } + } + } + + if( iPivot == -1 ) // Yikes! + { + CPLDebug( "AIRSAR", "No pivot in line `%s'.", + szLine ); + CPLAssert( iPivot != -1 ); + break; + } + +/* -------------------------------------------------------------------- */ +/* Trace ahead to the first non-white space value character. */ +/* -------------------------------------------------------------------- */ + int iValue = iPivot + 1; + + while( iValue < 50 && szLine[iValue] == ' ' ) + iValue++; + +/* -------------------------------------------------------------------- */ +/* Strip any white space off the keyword. */ +/* -------------------------------------------------------------------- */ + int iKeyEnd = iPivot - 1; + + while( iKeyEnd > 0 && szLine[iKeyEnd] == ' ' ) + iKeyEnd--; + + szLine[iKeyEnd+1] = '\0'; + +/* -------------------------------------------------------------------- */ +/* Convert spaces or colons into underscores in the key name. */ +/* -------------------------------------------------------------------- */ + for( i = 0; szLine[i] != '\0'; i++ ) + { + if( szLine[i] == ' ' || szLine[i] == ':' || szLine[i] == ',' ) + szLine[i] = '_'; + } + +/* -------------------------------------------------------------------- */ +/* Prefix key name with provided prefix string. */ +/* -------------------------------------------------------------------- */ + char szPrefixedKeyName[55]; + + sprintf( szPrefixedKeyName, "%s_%s", pszPrefix, szLine ); + + papszHeadInfo = + CSLSetNameValue( papszHeadInfo, szPrefixedKeyName, szLine+iValue ); + + } + + return papszHeadInfo; +} + + +/************************************************************************/ +/* Open() */ +/************************************************************************/ + +GDALDataset *AirSARDataset::Open( GDALOpenInfo * poOpenInfo ) + +{ + if( poOpenInfo->fp == NULL || poOpenInfo->nHeaderBytes < 800 ) + return NULL; + +/* -------------------------------------------------------------------- */ +/* Check for AirSAR/ keyword. */ +/* -------------------------------------------------------------------- */ + if( !EQUALN((char *) poOpenInfo->pabyHeader, "RECORD LENGTH IN BYTES",22) ) + return NULL; + + if( strstr((char *) poOpenInfo->pabyHeader, "COMPRESSED") == NULL + || strstr((char *) poOpenInfo->pabyHeader, "JPL AIRCRAFT") == NULL ) + return NULL; + +/* -------------------------------------------------------------------- */ +/* Parse the header fields. We turn all the transform the */ +/* keywords by converting spaces to underscores so they will be */ +/* "well behaved" as metadata keywords. */ +/* -------------------------------------------------------------------- */ + char **papszMD = ReadHeader( poOpenInfo->fp, 0, "MH", 20 ); + + if( papszMD == NULL ) + return NULL; + +/* -------------------------------------------------------------------- */ +/* Create a corresponding GDALDataset. */ +/* -------------------------------------------------------------------- */ + AirSARDataset *poDS; + + poDS = new AirSARDataset(); + +/* -------------------------------------------------------------------- */ +/* Extract some key information. */ +/* -------------------------------------------------------------------- */ + + poDS->nRasterXSize = + atoi(CSLFetchNameValue(papszMD,"MH_NUMBER_OF_SAMPLES_PER_RECORD")); + poDS->nRasterYSize = + atoi(CSLFetchNameValue(papszMD,"MH_NUMBER_OF_LINES_IN_IMAGE")); + + poDS->nRecordLength = atoi( + CSLFetchNameValue( papszMD, "MH_RECORD_LENGTH_IN_BYTES" ) ); + + poDS->nDataStart = atoi( + CSLFetchNameValue( papszMD, "MH_BYTE_OFFSET_OF_FIRST_DATA_RECORD" )); + +/* -------------------------------------------------------------------- */ +/* Adopt the openinfo file pointer. */ +/* -------------------------------------------------------------------- */ + poDS->fp = poOpenInfo->fp; + poOpenInfo->fp = NULL; + +/* -------------------------------------------------------------------- */ +/* Read and merge parameter header into metadata. Prefix */ +/* parameter header values with PH_. */ +/* -------------------------------------------------------------------- */ + int nPHOffset = 0; + + if( CSLFetchNameValue( papszMD, + "MH_BYTE_OFFSET_OF_PARAMETER_HEADER" ) != NULL ) + { + nPHOffset = atoi(CSLFetchNameValue( + papszMD, "MH_BYTE_OFFSET_OF_PARAMETER_HEADER")); + char **papszPHInfo = ReadHeader( poDS->fp, nPHOffset, "PH", 100 ); + + papszMD = CSLInsertStrings( papszMD, CSLCount(papszMD), papszPHInfo ); + + CSLDestroy( papszPHInfo ); + } + +/* -------------------------------------------------------------------- */ +/* Read and merge calibration header into metadata. Prefix */ +/* parameter header values with CH_. */ +/* -------------------------------------------------------------------- */ + if( nPHOffset != 0 ) + { + char **papszCHInfo = ReadHeader( poDS->fp, + nPHOffset+poDS->nRecordLength, + "CH", 18 ); + + papszMD = CSLInsertStrings( papszMD, CSLCount(papszMD), papszCHInfo ); + + CSLDestroy( papszCHInfo ); + } + +/* -------------------------------------------------------------------- */ +/* Assign metadata to dataset. */ +/* -------------------------------------------------------------------- */ + poDS->SetMetadata( papszMD ); + CSLDestroy( papszMD ); + +/* -------------------------------------------------------------------- */ +/* Create band information objects. */ +/* -------------------------------------------------------------------- */ + poDS->SetBand( 1, new AirSARRasterBand( poDS, 1 )); + poDS->SetBand( 2, new AirSARRasterBand( poDS, 2 )); + poDS->SetBand( 3, new AirSARRasterBand( poDS, 3 )); + poDS->SetBand( 4, new AirSARRasterBand( poDS, 4 )); + poDS->SetBand( 5, new AirSARRasterBand( poDS, 5 )); + poDS->SetBand( 6, new AirSARRasterBand( poDS, 6 )); + + poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename ); + +/* -------------------------------------------------------------------- */ +/* Initialize any PAM information. */ +/* -------------------------------------------------------------------- */ + poDS->SetDescription( poOpenInfo->pszFilename ); + poDS->TryLoadXML(); + + return( poDS ); +} + +/************************************************************************/ +/* GDALRegister_AirSAR() */ +/************************************************************************/ + +void GDALRegister_AirSAR() + +{ + GDALDriver *poDriver; + + if( GDALGetDriverByName( "AirSAR" ) == NULL ) + { + poDriver = new GDALDriver(); + + poDriver->SetDescription( "AirSAR" ); + poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, + "AirSAR Polarimetric Image" ); + poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "frmt_airsar.html" ); + poDriver->pfnOpen = AirSARDataset::Open; + + GetGDALDriverManager()->RegisterDriver( poDriver ); + } +} diff --git a/Utilities/GDAL/frmts/airsar/frmt_airsar.html b/Utilities/GDAL/frmts/airsar/frmt_airsar.html new file mode 100644 index 0000000000000000000000000000000000000000..b363cb81a311ff8a5b3699a2c4476753ac75c829 --- /dev/null +++ b/Utilities/GDAL/frmts/airsar/frmt_airsar.html @@ -0,0 +1,57 @@ +<html> +<head> +<title>AIRSAR -- AIRSAR Polarimetric Format</title> +</head> + +<body bgcolor="#ffffff"> + +<h1>AIRSAR -- AIRSAR Polarimetric Format</h1> + +Most variants of the AIRSAR Polarimetric Format produced by the AIRSAR +Integrated Processor are supported for reading by GDAL. AIRSAR products +normally include various associated data files, but only the imagery +data themselves is supported. Normally these are named <i>mission</i>_l.dat +(L-Band) or <i>mission</i>_c.dat (C-Band).<p> + +AIRSAR format contains a polarimetric image in compressed stokes matrix form. +Internally GDAL decompresses the data into a stokes matrix, and then converts +that form into a covariance matrix. The returned six bands are the six values +needed to define the 3x3 Hermitian covariance matrix. The convention used +to represent the covariance matrix in terms of the scattering matrix +elements HH, HV (=VH), and VV is indicated below. Note that the +non-diagonal elements of the matrix are complex values, while the diagonal +values are real (though represented as complex bands).<p> + +<ul> +<li> Band 1: Covariance_11 (Float32) = HH*conj(HH) +<li> Band 2: Covariance_12 (CFloat32) = sqrt(2)*HH*conj(HV) +<li> Band 3: Covariance_13 (CFloat32) = HH*conj(VV) +<li> Band 4: Covariance_22 (Float32) = 2*HV*conj(HV) +<li> Band 5: Covariance_23 (CFloat32) = sqrt(2)*HV*conj(VV) +<li> Band 6: Covariance_33 (Float32) = VV*conj(VV) +</ul> + + +The identities of the bands are also reflected in metadata and in the band +descriptions. <p> + +The AIRSAR product format includes (potentially) several headers of +information. This information is captured and represented as metadata on the +file as a whole. Information items from the main header are prefixed with +"MH_", items from the parameter header are prefixed with "PH_" and +information from the calibration header are prefixed with "CH_". +The metadata item names are derived automatically from the names of the +fields within the header itself. <p> + +No effort is made to read files associated with the AIRSAR product such as +<i>mission</i>_l.mocomp, <i>mission</i>_meta.airsar or +<i>mission</i>_meta.podaac.<p> + +See Also:<p> + +<ul> +<li> <a href="http://airsar.jpl.nasa.gov/documents/dataformat.htm">AIRSAR Data Format</a> +</ul> + +</body> +</html> diff --git a/Utilities/GDAL/frmts/airsar/makefile.vc b/Utilities/GDAL/frmts/airsar/makefile.vc new file mode 100644 index 0000000000000000000000000000000000000000..ad95a78c90dfefd4ba303a513ce4dc6ab9d4fbcc --- /dev/null +++ b/Utilities/GDAL/frmts/airsar/makefile.vc @@ -0,0 +1,13 @@ + +OBJ = airsardataset.obj + +GDAL_ROOT = ..\.. + +!INCLUDE $(GDAL_ROOT)\nmake.opt + +default: $(OBJ) + copy *.obj ..\o + +clean: + -del *.obj +