/******************************************************************************
* $Id: fastdataset.cpp,v 1.19 2005/12/20 20:07:30 dron Exp $
*
* Project: EOSAT FAST Format reader
* Purpose: Reads Landsat FAST-L7A, IRS 1C/1D
* Author: Andrey Kiselev, dron@at1895.spb.edu
*
******************************************************************************
* Copyright (c) 2002, Andrey Kiselev <dron@at1895.spb.edu>
*
* 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: fastdataset.cpp,v $
* Revision 1.19 2005/12/20 20:07:30 dron
* Added support for 16-bit data.
*
* Revision 1.18 2005/09/14 13:18:32 dron
* Avoid warnings.
*
* Revision 1.17 2005/05/05 13:55:42 fwarmerdam
* PAM Enable
*
* Revision 1.16 2004/07/12 18:24:23 gwalter
* Fixed geotransform calculation.
*
* Revision 1.15 2004/04/27 14:25:41 warmerda
* Cast to avoid warning on Solaris.
*
* Revision 1.14 2004/03/16 18:27:39 dron
* Fixes in projection parameters parsing code.
*
* Revision 1.13 2004/02/18 20:22:10 dron
* Create RawRasterBand objects in "large" mode; more datums and ellipsoids.
*
* Revision 1.12 2004/02/17 08:05:45 dron
* Do not calculate projection definition if corner coordinates are not set.
*
* Revision 1.11 2004/02/03 20:38:50 dron
* Fixes in coordinate hadling; recognize more projections.
*
* Revision 1.10 2004/02/01 17:29:34 dron
* Format parsing logic completely rewritten. Start using importFromUSGS().
*
* Revision 1.9 2003/10/17 07:08:21 dron
* Use locale selection option in CPLScanDouble().
*
* Revision 1.8 2003/10/05 15:31:00 dron
* TM projection support implemented.
*
* Revision 1.7 2003/07/08 21:10:19 warmerda
* avoid warnings
*
* Revision 1.6 2003/03/14 17:28:10 dron
* CPLFormCIFilename() used instead of CPLFormFilename() for FAST-L7 datasets.
*
* Revision 1.5 2003/03/05 15:49:59 dron
* Fixed typo when reading SENSOR metadata record.
*
* Revision 1.4 2003/02/18 15:07:49 dron
* IRS-1C/1D support added.
*
* Revision 1.3 2003/02/14 21:05:40 warmerda
* Don't use path for rawdataset.h.
*
* Revision 1.2 2002/12/30 14:55:01 dron
* SetProjCS() removed, added unit setting.
*
* Revision 1.1 2002/10/05 12:35:31 dron
* FAST driver moved to the RAW directory.
*
* Revision 1.5 2002/10/04 16:06:06 dron
* Some redundancy removed.
*
* Revision 1.4 2002/10/04 12:33:02 dron
* Added calibration coefficients extraction.
*
* Revision 1.3 2002/09/04 06:50:37 warmerda
* avoid static driver pointers
*
* Revision 1.2 2002/08/15 09:35:50 dron
* Fixes in georeferencing
*
* Revision 1.1 2002/08/13 16:55:41 dron
* Initial release
*
*
*/
#include "cpl_string.h"
#include "cpl_conv.h"
#include "ogr_spatialref.h"
#include "rawdataset.h"
CPL_CVSID("$Id: fastdataset.cpp,v 1.19 2005/12/20 20:07:30 dron Exp $");
CPL_C_START
void GDALRegister_FAST(void);
CPL_C_END
#define ADM_STD_HEADER_SIZE 4608 // XXX: Format specification says it
#define ADM_HEADER_SIZE 5000 // should be 4608, but some vendors
// ship broken large datasets.
#define ACQUISITION_DATE "ACQUISITION DATE ="
#define ACQUISITION_DATE_SIZE 8
#define SATELLITE_NAME "SATELLITE ="
#define SATELLITE_NAME_SIZE 10
#define SENSOR_NAME "SENSOR ="
#define SENSOR_NAME_SIZE 10
#define FILENAME "FILENAME ="
#define FILENAME_SIZE 29
#define PIXELS "PIXELS PER LINE ="
#define PIXELS_SIZE 5
#define LINES "LINES PER BAND ="
#define LINES_SIZE 5
#define BITS_PER_PIXEL "OUTPUT BITS PER PIXEL ="
#define BITS_PER_PIXEL_SIZE 2
#define PROJECTION_NAME "MAP PROJECTION ="
#define PROJECTION_NAME_SIZE 4
#define ELLIPSOID_NAME "ELLIPSOID ="
#define ELLIPSOID_NAME_SIZE 18
#define DATUM_NAME "DATUM ="
#define DATUM_NAME_SIZE 6
#define ZONE_NUMBER "USGS MAP ZONE ="
#define ZONE_NUMBER_SIZE 6
#define USGS_PARAMETERS "USGS PROJECTION PARAMETERS ="
#define CORNER_UPPER_LEFT "UL ="
#define CORNER_UPPER_RIGHT "UR ="
#define CORNER_LOWER_LEFT "LL ="
#define CORNER_LOWER_RIGHT "LR ="
#define CORNER_VALUE_SIZE 13
#define VALUE_SIZE 24
enum FASTSatellite // Satellites:
{
LANDSAT, // Landsat 7
IRS // IRS 1C/1D
};
/************************************************************************/
/* ==================================================================== */
/* FASTDataset */
/* ==================================================================== */
/************************************************************************/
class FASTDataset : public GDALPamDataset
{
friend class FASTRasterBand;
double adfGeoTransform[6];
char *pszProjection;
FILE *fpHeader;
FILE *fpChannels[6];
const char *pszFilename;
char *pszDirname;
GDALDataType eDataType;
FASTSatellite iSatellite;
public:
FASTDataset();
~FASTDataset();
static GDALDataset *Open( GDALOpenInfo * );
CPLErr GetGeoTransform( double * padfTransform );
const char *GetProjectionRef();
FILE *FOpenChannel( char *pszFilename, int iBand );
};
/************************************************************************/
/* ==================================================================== */
/* FASTRasterBand */
/* ==================================================================== */
/************************************************************************/
class FASTRasterBand : public RawRasterBand
{
friend class FASTDataset;
public:
FASTRasterBand( FASTDataset *, int, FILE *, vsi_l_offset,
int, int, GDALDataType, int );
};
/************************************************************************/
/* FASTRasterBand() */
/************************************************************************/
FASTRasterBand::FASTRasterBand( FASTDataset *poDS, int nBand, FILE * fpRaw,
vsi_l_offset nImgOffset, int nPixelOffset,
int nLineOffset, GDALDataType eDataType,
int bNativeOrder) :
RawRasterBand( poDS, nBand, fpRaw, nImgOffset, nPixelOffset,
nLineOffset, eDataType, bNativeOrder, TRUE)
{
}
/************************************************************************/
/* ==================================================================== */
/* FASTDataset */
/* ==================================================================== */
/************************************************************************/
/************************************************************************/
/* FASTDataset() */
/************************************************************************/
FASTDataset::FASTDataset()
{
fpHeader = NULL;
pszDirname = NULL;
pszProjection = CPLStrdup( "" );
adfGeoTransform[0] = 0.0;
adfGeoTransform[1] = 1.0;
adfGeoTransform[2] = 0.0;
adfGeoTransform[3] = 0.0;
adfGeoTransform[4] = 0.0;
adfGeoTransform[5] = 1.0;
nBands = 0;
}
/************************************************************************/
/* ~FASTDataset() */
/************************************************************************/
FASTDataset::~FASTDataset()
{
int i;
FlushCache();
if ( pszDirname )
CPLFree( pszDirname );
if ( pszProjection )
CPLFree( pszProjection );
for ( i = 0; i < nBands; i++ )
if ( fpChannels[i] )
VSIFCloseL( fpChannels[i] );
if( fpHeader != NULL )
VSIFClose( fpHeader );
}
/************************************************************************/
/* GetGeoTransform() */
/************************************************************************/
CPLErr FASTDataset::GetGeoTransform( double * padfTransform )
{
memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
return CE_None;
}
/************************************************************************/
/* GetProjectionRef() */
/************************************************************************/
const char *FASTDataset::GetProjectionRef()
{
if( pszProjection )
return pszProjection;
else
return "";
}
/************************************************************************/
/* FOpenChannel() */
/************************************************************************/
FILE *FASTDataset::FOpenChannel( char *pszFilename, int iBand )
{
const char *pszChannelFilename = NULL;
char *pszPrefix = CPLStrdup( CPLGetBasename( this->pszFilename ) );
char *pszSuffix = CPLStrdup( CPLGetExtension( this->pszFilename ) );
switch ( iSatellite )
{
case LANDSAT:
if ( pszFilename && !EQUAL( pszFilename, "" ) )
{
pszChannelFilename =
CPLFormCIFilename( pszDirname, pszFilename, NULL );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
}
else
fpChannels[iBand] = NULL;
break;
case IRS:
default:
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "%s.%d", pszPrefix, iBand + 1 ), pszSuffix );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "IMAGERY%d", iBand + 1 ), pszSuffix );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "imagery%d", iBand + 1 ), pszSuffix );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "IMAGERY%d.DAT", iBand + 1 ), NULL );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "imagery%d.dat", iBand + 1 ), NULL );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "IMAGERY%d.dat", iBand + 1 ), NULL );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "imagery%d.DAT", iBand + 1 ), NULL );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "BAND%d", iBand + 1 ), pszSuffix );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "band%d", iBand + 1 ), pszSuffix );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "BAND%d.DAT", iBand + 1 ), NULL );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "band%d.dat", iBand + 1 ), NULL );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "BAND%d.dat", iBand + 1 ), NULL );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
if ( fpChannels[iBand] )
break;
pszChannelFilename = CPLFormFilename( pszDirname,
CPLSPrintf( "band%d.DAT", iBand + 1 ), NULL );
fpChannels[iBand] = VSIFOpenL( pszChannelFilename, "rb" );
break;
}
CPLDebug( "FAST", "Band %d filename=%s", iBand + 1, pszChannelFilename);
CPLFree( pszPrefix );
CPLFree( pszSuffix );
return fpChannels[iBand];
}
/************************************************************************/
/* GetValue() */
/************************************************************************/
static char *GetValue( const char *pszString, const char *pszName,
int iValueSize, int iNormalize )
{
char *pszTemp = strstr( (char *) pszString, pszName );
if ( pszTemp )
{
pszTemp += strlen( pszName );
pszTemp = CPLScanString( pszTemp, iValueSize, TRUE, iNormalize );
}
return pszTemp;
}
/************************************************************************/
/* USGSMnemonicToCode() */
/************************************************************************/
static long USGSMnemonicToCode( const char* pszMnemonic )
{
if ( EQUAL(pszMnemonic, "UTM") )
return 1L;
else if ( EQUAL(pszMnemonic, "LCC") )
return 4L;
else if ( EQUAL(pszMnemonic, "PS") )
return 6L;
else if ( EQUAL(pszMnemonic, "PC") )
return 7L;
else if ( EQUAL(pszMnemonic, "TM") )
return 9L;
else if ( EQUAL(pszMnemonic, "OM") )
return 20L;
else if ( EQUAL(pszMnemonic, "SOM") )
return 22L;
else
return 1L; // UTM by default
}
/************************************************************************/
/* USGSEllipsoidToCode() */
/************************************************************************/
static long USGSEllipsoidToCode( const char* pszMnemonic )
{
if ( EQUAL(pszMnemonic, "CLARKE_1866") )
return 0L;
else if ( EQUAL(pszMnemonic, "CLARKE_1880") )
return 1L;
else if ( EQUAL(pszMnemonic, "BESSEL") )
return 2L;
else if ( EQUAL(pszMnemonic, "INTERNATL_1967") )
return 3L;
else if ( EQUAL(pszMnemonic, "INTERNATL_1909") )
return 4L;
else if ( EQUAL(pszMnemonic, "WGS72") || EQUAL(pszMnemonic, "WGS_72") )
return 5L;
else if ( EQUAL(pszMnemonic, "EVEREST") )
return 6L;
else if ( EQUAL(pszMnemonic, "WGS66") || EQUAL(pszMnemonic, "WGS_66") )
return 7L;
else if ( EQUAL(pszMnemonic, "GRS_80") )
return 8L;
else if ( EQUAL(pszMnemonic, "AIRY") )
return 9L;
else if ( EQUAL(pszMnemonic, "MODIFIED_EVEREST") )
return 10L;
else if ( EQUAL(pszMnemonic, "MODIFIED_AIRY") )
return 11L;
else if ( EQUAL(pszMnemonic, "WGS84") || EQUAL(pszMnemonic, "WGS_84") )
return 12L;
else if ( EQUAL(pszMnemonic, "SOUTHEAST_ASIA") )
return 13L;
else if ( EQUAL(pszMnemonic, "AUSTRALIAN_NATL") )
return 14L;
else if ( EQUAL(pszMnemonic, "KRASSOVSKY") )
return 15L;
else if ( EQUAL(pszMnemonic, "HOUGH") )
return 16L;
else if ( EQUAL(pszMnemonic, "MERCURY_1960") )
return 17L;
else if ( EQUAL(pszMnemonic, "MOD_MERC_1968") )
return 18L;
else if ( EQUAL(pszMnemonic, "6370997_M_SPHERE") )
return 19L;
else
return 0L;
}
/************************************************************************/
/* Open() */
/************************************************************************/
GDALDataset *FASTDataset::Open( GDALOpenInfo * poOpenInfo )
{
int i;
if( poOpenInfo->fp == NULL )
return NULL;
if( !EQUALN((const char *) poOpenInfo->pabyHeader + 52,
"ACQUISITION DATE =", 18) )
return NULL;
/* -------------------------------------------------------------------- */
/* Create a corresponding GDALDataset. */
/* -------------------------------------------------------------------- */
FASTDataset *poDS;
poDS = new FASTDataset();
poDS->fpHeader = poOpenInfo->fp;
poOpenInfo->fp = NULL;
poDS->pszFilename = poOpenInfo->pszFilename;
poDS->pszDirname = CPLStrdup( CPLGetDirname( poOpenInfo->pszFilename ) );
/* -------------------------------------------------------------------- */
/* Read the administrative record. */
/* -------------------------------------------------------------------- */
char *pszTemp;
char *pszHeader = (char *) CPLMalloc( ADM_HEADER_SIZE + 1 );
size_t nBytesRead;
VSIFSeek( poDS->fpHeader, 0, SEEK_SET );
nBytesRead = VSIFRead( pszHeader, 1, ADM_HEADER_SIZE, poDS->fpHeader );
if ( nBytesRead < ADM_STD_HEADER_SIZE )
{
CPLDebug( "FAST", "Header file too short. Reading failed" );
delete poDS;
return NULL;
}
pszHeader[nBytesRead] = '\0';
// Read acquisition date
pszTemp = GetValue( pszHeader, ACQUISITION_DATE,
ACQUISITION_DATE_SIZE, TRUE );
poDS->SetMetadataItem( "ACQUISITION_DATE", pszTemp );
CPLFree( pszTemp );
// Read satellite name (will read the first one only)
pszTemp = GetValue( pszHeader, SATELLITE_NAME, SATELLITE_NAME_SIZE, TRUE );
poDS->SetMetadataItem( "SATELLITE", pszTemp );
if ( EQUALN(pszTemp, "LANDSAT", 7) )
poDS->iSatellite = LANDSAT;
else if ( EQUALN(pszTemp, "IRS", 3) )
poDS->iSatellite = IRS;
else
poDS->iSatellite = IRS;
CPLFree( pszTemp );
// Read sensor name (will read the first one only)
pszTemp = GetValue( pszHeader, SENSOR_NAME, SENSOR_NAME_SIZE, TRUE );
poDS->SetMetadataItem( "SENSOR", pszTemp );
CPLFree( pszTemp );
// Read filenames
pszTemp = pszHeader;
poDS->nBands = 0;
for ( i = 0; i < 6; i++ )
{
char *pszFilename = NULL ;
if ( pszTemp )
pszTemp = strstr( pszTemp, FILENAME );
if ( pszTemp )
{
pszTemp += strlen(FILENAME);
pszFilename = CPLScanString( pszTemp, FILENAME_SIZE, TRUE, FALSE );
}
else
pszTemp = NULL;
if ( poDS->FOpenChannel( pszFilename, poDS->nBands ) )
poDS->nBands++;
if ( pszFilename )
CPLFree( pszFilename );
}
if ( !poDS->nBands )
{
CPLDebug( "FAST", "Failed to find and open band data files." );
delete poDS;
return NULL;
}
// Read number of pixels/lines and bit depth
pszTemp = strstr( pszHeader, PIXELS );
if ( pszTemp )
poDS->nRasterXSize = CPLScanLong( pszTemp + strlen(PIXELS),
PIXELS_SIZE );
else
{
CPLDebug( "FAST", "Failed to find number of pixels in line." );
delete poDS;
return NULL;
}
pszTemp = strstr( pszHeader, LINES );
if ( pszTemp )
poDS->nRasterYSize = CPLScanLong( pszTemp + strlen(LINES), LINES_SIZE );
else
{
CPLDebug( "FAST", "Failed to find number of lines in raster." );
delete poDS;
return NULL;
}
pszTemp = strstr( pszHeader, BITS_PER_PIXEL );
switch( CPLScanLong(pszTemp+strlen(BITS_PER_PIXEL), BITS_PER_PIXEL_SIZE) )
{
case 8:
default:
poDS->eDataType = GDT_Byte;
break;
case 16:
poDS->eDataType = GDT_UInt16;
break;
}
/* -------------------------------------------------------------------- */
/* Read radiometric record. */
/* -------------------------------------------------------------------- */
// Read gains and biases. This is a trick!
pszTemp = strstr( pszHeader, "BIASES" );// It may be "BIASES AND GAINS"
// or "GAINS AND BIASES"
// Now search for the first number occurance after that string
for ( i = 1; i <= poDS->nBands; i++ )
{
char *pszValue = NULL;
pszTemp = strpbrk( pszTemp, "-.0123456789" );
if ( pszTemp )
{
pszValue = CPLScanString( pszTemp, VALUE_SIZE, TRUE, TRUE );
poDS->SetMetadataItem( CPLSPrintf("BIAS%d", i ), pszValue );
}
pszTemp += VALUE_SIZE;
if ( pszValue )
CPLFree( pszValue );
pszTemp = strpbrk( pszTemp, "-.0123456789" );
if ( pszTemp )
{
pszValue = CPLScanString( pszTemp, VALUE_SIZE, TRUE, TRUE );
poDS->SetMetadataItem( CPLSPrintf("GAIN%d", i ), pszValue );
}
pszTemp += VALUE_SIZE;
if ( pszValue )
CPLFree( pszValue );
}
/* -------------------------------------------------------------------- */
/* Read geometric record. */
/* -------------------------------------------------------------------- */
OGRSpatialReference oSRS;
long iProjSys, iZone, iDatum;
// Coordinates of pixel's centers
double dfULX = 0.0, dfULY = 0.0;
double dfURX = 0.0, dfURY = 0.0;
double dfLLX = 0.0, dfLLY = 0.0;
double dfLRX = 0.0, dfLRY = 0.0;
double adfProjParms[15];
// Read projection name
pszTemp = GetValue( pszHeader, PROJECTION_NAME,
PROJECTION_NAME_SIZE, FALSE );
if ( pszTemp && !EQUAL( pszTemp, "" ) )
iProjSys = USGSMnemonicToCode( pszTemp );
else
iProjSys = 1L; // UTM by default
CPLFree( pszTemp );
// Read ellipsoid name
pszTemp = GetValue( pszHeader, ELLIPSOID_NAME, ELLIPSOID_NAME_SIZE, FALSE );
if ( pszTemp && !EQUAL( pszTemp, "" ) )
iDatum = USGSEllipsoidToCode( pszTemp );
else
iDatum = 0L; // Clarke, 1866 (NAD1927) by default
CPLFree( pszTemp );
// Read zone number
pszTemp = GetValue( pszHeader, ZONE_NUMBER, ZONE_NUMBER_SIZE, FALSE );
if ( pszTemp && !EQUAL( pszTemp, "" ) )
iZone = atoi( pszTemp );
else
iZone = 0L;
CPLFree( pszTemp );
// Read 15 USGS projection parameters
for ( i = 0; i < 15; i++ )
adfProjParms[i] = 0.0;
pszTemp = strstr( pszHeader, USGS_PARAMETERS );
if ( pszTemp && !EQUAL( pszTemp, "" ) )
{
pszTemp += strlen( USGS_PARAMETERS );
for ( i = 0; i < 15; i++ )
{
pszTemp = strpbrk( pszTemp, "-.0123456789" );
if ( pszTemp )
{
adfProjParms[i] = CPLScanDouble( pszTemp, VALUE_SIZE, "C" );
#if DEBUG
CPLDebug("FAST", "USGS parameter %2d=%f.", i, adfProjParms[i]);
#endif
}
pszTemp = strpbrk( pszTemp, " \t" );
}
}
// Read corner coordinates
pszTemp = strstr( pszHeader, CORNER_UPPER_LEFT );
if ( pszTemp && !EQUAL( pszTemp, "" ) )
{
pszTemp += strlen( CORNER_UPPER_LEFT ) + 28;
dfULX = CPLScanDouble( pszTemp, CORNER_VALUE_SIZE, "C" );
pszTemp += CORNER_VALUE_SIZE + 1;
dfULY = CPLScanDouble( pszTemp, CORNER_VALUE_SIZE, "C" );
}
pszTemp = strstr( pszHeader, CORNER_UPPER_RIGHT );
if ( pszTemp && !EQUAL( pszTemp, "" ) )
{
pszTemp += strlen( CORNER_UPPER_RIGHT ) + 28;
dfURX = CPLScanDouble( pszTemp, CORNER_VALUE_SIZE, "C" );
pszTemp += CORNER_VALUE_SIZE + 1;
dfURY = CPLScanDouble( pszTemp, CORNER_VALUE_SIZE, "C" );
}
pszTemp = strstr( pszHeader, CORNER_LOWER_LEFT );
if ( pszTemp && !EQUAL( pszTemp, "" ) )
{
pszTemp += strlen( CORNER_LOWER_LEFT ) + 28;
dfLLX = CPLScanDouble( pszTemp, CORNER_VALUE_SIZE, "C" );
pszTemp += CORNER_VALUE_SIZE + 1;
dfLLY = CPLScanDouble( pszTemp, CORNER_VALUE_SIZE, "C" );
}
pszTemp = strstr( pszHeader, CORNER_LOWER_RIGHT );
if ( pszTemp && !EQUAL( pszTemp, "" ) )
{
pszTemp += strlen( CORNER_LOWER_RIGHT ) + 28;
dfLRX = CPLScanDouble( pszTemp, CORNER_VALUE_SIZE, "C" );
pszTemp += CORNER_VALUE_SIZE + 1;
dfLRY = CPLScanDouble( pszTemp, CORNER_VALUE_SIZE, "C" );
}
if ( dfULX != 0.0 && dfULY != 0.0
&& dfURX != 0.0 && dfURY != 0.0
&& dfLLX != 0.0 && dfLLY != 0.0
&& dfLRX != 0.0 && dfLRY != 0.0 )
{
int transform_ok=FALSE;
GDAL_GCP *pasGCPList;
// Strip out zone number from the easting values, if either
if ( dfULX >= 1000000.0 )
dfULX -= (double)iZone * 1000000.0;
if ( dfURX >= 1000000.0 )
dfURX -= (double)iZone * 1000000.0;
if ( dfLLX >= 1000000.0 )
dfLLX -= (double)iZone * 1000000.0;
if ( dfLRX >= 1000000.0 )
dfLRX -= (double)iZone * 1000000.0;
// Create projection definition
oSRS.importFromUSGS( iProjSys, iZone, adfProjParms, iDatum );
oSRS.SetLinearUnits( SRS_UL_METER, 1.0 );
// Read datum name
pszTemp = GetValue( pszHeader, DATUM_NAME, DATUM_NAME_SIZE, FALSE );
if ( EQUAL( pszTemp, "WGS84" ) )
oSRS.SetWellKnownGeogCS( "WGS84" );
else if ( EQUAL( pszTemp, "NAD27" ) )
oSRS.SetWellKnownGeogCS( "NAD27" );
else if ( EQUAL( pszTemp, "NAD83" ) )
oSRS.SetWellKnownGeogCS( "NAD83" );
CPLFree( pszTemp );
if ( poDS->pszProjection )
CPLFree( poDS->pszProjection );
oSRS.exportToWkt( &poDS->pszProjection );
// Generate GCPs
pasGCPList = (GDAL_GCP *) CPLCalloc( sizeof( GDAL_GCP ), 4 );
GDALInitGCPs( 4, pasGCPList );
pasGCPList[0].pszId = "UPPER_LEFT";
pasGCPList[0].dfGCPX = dfULX;
pasGCPList[0].dfGCPY = dfULY;
pasGCPList[0].dfGCPZ = 0.0;
pasGCPList[0].dfGCPPixel = 0.5;
pasGCPList[0].dfGCPLine = 0.5;
pasGCPList[1].pszId = "UPPER_RIGHT";
pasGCPList[1].dfGCPX = dfURX;
pasGCPList[1].dfGCPY = dfURY;
pasGCPList[1].dfGCPZ = 0.0;
pasGCPList[1].dfGCPPixel = poDS->nRasterXSize-0.5;
pasGCPList[1].dfGCPLine = 0.5;
pasGCPList[2].pszId = "LOWER_LEFT";
pasGCPList[2].dfGCPX = dfLLX;
pasGCPList[2].dfGCPY = dfLLY;
pasGCPList[2].dfGCPZ = 0.0;
pasGCPList[2].dfGCPPixel = 0.5;
pasGCPList[2].dfGCPLine = poDS->nRasterYSize-0.5;
pasGCPList[3].pszId = "LOWER_RIGHT";
pasGCPList[3].dfGCPX = dfLRX;
pasGCPList[3].dfGCPY = dfLRY;
pasGCPList[3].dfGCPZ = 0.0;
pasGCPList[3].dfGCPPixel = poDS->nRasterXSize-0.5;
pasGCPList[3].dfGCPLine = poDS->nRasterYSize-0.5;
// Calculate transformation matrix, if accurate
transform_ok = GDALGCPsToGeoTransform(4,pasGCPList,poDS->adfGeoTransform,0);
if (transform_ok == FALSE)
{
poDS->adfGeoTransform[0] = 0.0;
poDS->adfGeoTransform[1] = 1.0;
poDS->adfGeoTransform[2] = 0.0;
poDS->adfGeoTransform[3] = 0.0;
poDS->adfGeoTransform[4] = 0.0;
poDS->adfGeoTransform[5] = 1.0;
if ( poDS->pszProjection )
CPLFree( poDS->pszProjection );
poDS->pszProjection = CPLStrdup("");
}
}
/* -------------------------------------------------------------------- */
/* Create band information objects. */
/* -------------------------------------------------------------------- */
int nPixelOffset = GDALGetDataTypeSize(poDS->eDataType) / 8;
int nLineOffset = poDS->nRasterXSize * nPixelOffset;
for( i = 1; i <= poDS->nBands; i++ )
poDS->SetBand( i, new FASTRasterBand( poDS, i, poDS->fpChannels[i - 1],
0, nPixelOffset, nLineOffset, poDS->eDataType, TRUE));
CPLFree( pszHeader );
/* -------------------------------------------------------------------- */
/* Initialize any PAM information. */
/* -------------------------------------------------------------------- */
poDS->SetDescription( poOpenInfo->pszFilename );
poDS->TryLoadXML();
return( poDS );
}
/************************************************************************/
/* GDALRegister_FAST() */
/************************************************************************/
void GDALRegister_FAST()
{
GDALDriver *poDriver;
if( GDALGetDriverByName( "FAST" ) == NULL )
{
poDriver = new GDALDriver();
poDriver->SetDescription( "FAST" );
poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
"EOSAT FAST Format" );
poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
"frmt_fast.html" );
poDriver->pfnOpen = FASTDataset::Open;
GetGDALDriverManager()->RegisterDriver( poDriver );
}
}
