/****************************************************************************** * $Id: btdataset.cpp,v 1.9 2006/04/11 00:33:15 fwarmerdam Exp $ * * Project: VTP .bt Driver * Purpose: Implementation of VTP .bt elevation format read/write support. * http://www.vterrain.org/Implementation/Formats/BT.html * Author: Frank Warmerdam, warmerdam@pobox.com * ****************************************************************************** * Copyright (c) 2003, Frank Warmerdam * * 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: btdataset.cpp,v $ * Revision 1.9 2006/04/11 00:33:15 fwarmerdam * Added SetUnitType(), contributed by Ray Gardener. * * Revision 1.8 2006/04/09 04:22:28 fwarmerdam * Implement GetUnitType() support when reading, from Ray Gardner. * * Revision 1.7 2005/05/05 13:55:41 fwarmerdam * PAM Enable * * Revision 1.6 2004/05/11 17:27:33 warmerda * Set eAccess properly in Open(). * * Revision 1.5 2004/03/02 15:14:07 warmerda * Fixed CPL_MSB case. * * Revision 1.4 2004/03/01 17:55:10 warmerda * rewrote to *not* use raw layer * * Revision 1.3 2004/01/08 20:13:33 warmerda * Fixed bugs with georeferencing (on read), and with imagery access * (was diagonally flipped). Note, that BT imagery is column organized, * like USGS DEMs, not row oriented as the raw driver is optimized for. * * Revision 1.2 2003/12/11 22:44:11 warmerda * ensure external .prj flag gets set * * Revision 1.1 2003/12/11 06:18:18 warmerda * New * */ #include "rawdataset.h" #include "ogr_spatialref.h" CPL_CVSID("$Id: btdataset.cpp,v 1.9 2006/04/11 00:33:15 fwarmerdam Exp $"); CPL_C_START void GDALRegister_BT(void); CPL_C_END /************************************************************************/ /* ==================================================================== */ /* BTDataset */ /* ==================================================================== */ /************************************************************************/ class BTDataset : public GDALPamDataset { friend class BTRasterBand; FILE *fpImage; // image data file. int bGeoTransformValid; double adfGeoTransform[6]; char *pszProjection; double dfVScale; int nVersionCode; // version times 10. int bHeaderModified; unsigned char abyHeader[256]; float m_fVscale; public: BTDataset(); ~BTDataset(); virtual const char *GetProjectionRef(void); virtual CPLErr SetProjection( const char * ); virtual CPLErr GetGeoTransform( double * ); virtual CPLErr SetGeoTransform( double * ); virtual void FlushCache(); static GDALDataset *Open( GDALOpenInfo * ); static GDALDataset *Create( const char * pszFilename, int nXSize, int nYSize, int nBands, GDALDataType eType, char ** papszOptions ); }; /************************************************************************/ /* ==================================================================== */ /* BTRasterBand */ /* ==================================================================== */ /************************************************************************/ class BTRasterBand : public GDALRasterBand { FILE *fpImage; public: BTRasterBand( GDALDataset * poDS, FILE * fp, GDALDataType eType ); virtual CPLErr IReadBlock( int, int, void * ); virtual CPLErr IWriteBlock( int, int, void * ); virtual const char* GetUnitType(); virtual CPLErr SetUnitType(const char*); }; /************************************************************************/ /* BTRasterBand() */ /************************************************************************/ BTRasterBand::BTRasterBand( GDALDataset *poDS, FILE *fp, GDALDataType eType ) { this->poDS = poDS; this->nBand = 1; this->eDataType = eType; this->fpImage = fp; nBlockXSize = 1; nBlockYSize = poDS->GetRasterYSize(); } /************************************************************************/ /* IReadBlock() */ /************************************************************************/ CPLErr BTRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff, void * pImage ) { int nDataSize = GDALGetDataTypeSize( eDataType ) / 8; int i; CPLAssert( nBlockYOff == 0 ); /* -------------------------------------------------------------------- */ /* Seek to profile. */ /* -------------------------------------------------------------------- */ if( VSIFSeekL( fpImage, 256 + nBlockXOff * nDataSize * nRasterYSize, SEEK_SET ) != 0 ) { CPLError( CE_Failure, CPLE_FileIO, ".bt Seek failed:%s", VSIStrerror( errno ) ); return CE_Failure; } /* -------------------------------------------------------------------- */ /* Read the profile. */ /* -------------------------------------------------------------------- */ if( VSIFReadL( pImage, nDataSize, nRasterYSize, fpImage ) != (size_t) nRasterYSize ) { CPLError( CE_Failure, CPLE_FileIO, ".bt Read failed:%s", VSIStrerror( errno ) ); return CE_Failure; } /* -------------------------------------------------------------------- */ /* Swap on MSB platforms. */ /* -------------------------------------------------------------------- */ #ifdef CPL_MSB GDALSwapWords( pImage, nDataSize, nRasterYSize, nDataSize ); #endif /* -------------------------------------------------------------------- */ /* Vertical flip, since GDAL expects values from top to bottom, */ /* but in .bt they are bottom to top. */ /* -------------------------------------------------------------------- */ for( i = 0; i < nRasterYSize / 2; i++ ) { GByte abyWrk[8]; memcpy( abyWrk, ((GByte *) pImage) + i * nDataSize, nDataSize ); memcpy( ((GByte *) pImage) + i * nDataSize, ((GByte *) pImage) + (nRasterYSize - i - 1) * nDataSize, nDataSize ); memcpy( ((GByte *) pImage) + (nRasterYSize - i - 1) * nDataSize, abyWrk, nDataSize ); } return CE_None; } /************************************************************************/ /* IWriteBlock() */ /************************************************************************/ CPLErr BTRasterBand::IWriteBlock( int nBlockXOff, int nBlockYOff, void * pImage ) { int nDataSize = GDALGetDataTypeSize( eDataType ) / 8; GByte *pabyWrkBlock; int i; CPLAssert( nBlockYOff == 0 ); /* -------------------------------------------------------------------- */ /* Seek to profile. */ /* -------------------------------------------------------------------- */ if( VSIFSeekL( fpImage, 256 + nBlockXOff * nDataSize * nRasterYSize, SEEK_SET ) != 0 ) { CPLError( CE_Failure, CPLE_FileIO, ".bt Seek failed:%s", VSIStrerror( errno ) ); return CE_Failure; } /* -------------------------------------------------------------------- */ /* Allocate working buffer. */ /* -------------------------------------------------------------------- */ pabyWrkBlock = (GByte *) CPLMalloc(nDataSize * nRasterYSize); /* -------------------------------------------------------------------- */ /* Vertical flip data into work buffer, since GDAL expects */ /* values from top to bottom, but in .bt they are bottom to */ /* top. */ /* -------------------------------------------------------------------- */ for( i = 0; i < nRasterYSize; i++ ) { memcpy( pabyWrkBlock + (nRasterYSize - i - 1) * nDataSize, ((GByte *) pImage) + i * nDataSize, nDataSize ); } /* -------------------------------------------------------------------- */ /* Swap on MSB platforms. */ /* -------------------------------------------------------------------- */ #ifdef CPL_MSB GDALSwapWords( pabyWrkBlock, nDataSize, nRasterYSize, nDataSize ); #endif /* -------------------------------------------------------------------- */ /* Read the profile. */ /* -------------------------------------------------------------------- */ if( VSIFWriteL( pabyWrkBlock, nDataSize, nRasterYSize, fpImage ) != (size_t) nRasterYSize ) { CPLFree( pabyWrkBlock ); CPLError( CE_Failure, CPLE_FileIO, ".bt Write failed:%s", VSIStrerror( errno ) ); return CE_Failure; } CPLFree( pabyWrkBlock ); return CE_None; } /************************************************************************/ /* GetUnitType() */ /************************************************************************/ static bool approx_equals(float a, float b) { const float epsilon = (float)1e-5; return (fabs(a-b) <= epsilon); } const char* BTRasterBand::GetUnitType(void) { const BTDataset& ds = *(BTDataset*)poDS; float f = ds.m_fVscale; if(f == 1.0f) return "m"; if(approx_equals(f, 0.3048f)) return "ft"; if(approx_equals(f, 1200.0f/3937.0f)) return "sft"; // todo: the BT spec allows for any value for // ds.m_fVscale, so rigorous adherence would // require testing for all possible units you // may want to support, such as km, yards, miles, etc. // But m/ft/sft seem to be the top three. return ""; } /************************************************************************/ /* SetUnitType() */ /************************************************************************/ CPLErr BTRasterBand::SetUnitType(const char* psz) { BTDataset& ds = *(BTDataset*)poDS; if(EQUAL(psz, "m")) ds.m_fVscale = 1.0f; else if(EQUAL(psz, "ft")) ds.m_fVscale = 0.3048f; else if(EQUAL(psz, "sft")) ds.m_fVscale = 1200.0f / 3937.0f; else return CE_Failure; float fScale = ds.m_fVscale; CPL_LSBPTR32(&fScale); // Update header's elevation scale field. memcpy(ds.abyHeader + 62, &fScale, sizeof(fScale)); ds.bHeaderModified = TRUE; return CE_None; } /************************************************************************/ /* ==================================================================== */ /* BTDataset */ /* ==================================================================== */ /************************************************************************/ /************************************************************************/ /* BTDataset() */ /************************************************************************/ BTDataset::BTDataset() { fpImage = NULL; bGeoTransformValid = FALSE; adfGeoTransform[0] = 0.0; adfGeoTransform[1] = 1.0; adfGeoTransform[2] = 0.0; adfGeoTransform[3] = 0.0; adfGeoTransform[4] = 0.0; adfGeoTransform[5] = 1.0; pszProjection = NULL; bHeaderModified = FALSE; } /************************************************************************/ /* ~BTDataset() */ /************************************************************************/ BTDataset::~BTDataset() { FlushCache(); if( fpImage != NULL ) VSIFCloseL( fpImage ); CPLFree( pszProjection ); } /************************************************************************/ /* FlushCache() */ /* */ /* We override this to include flush out the header block. */ /************************************************************************/ void BTDataset::FlushCache() { GDALDataset::FlushCache(); if( !bHeaderModified ) return; bHeaderModified = FALSE; VSIFSeekL( fpImage, 0, SEEK_SET ); VSIFWriteL( abyHeader, 256, 1, fpImage ); } /************************************************************************/ /* GetGeoTransform() */ /************************************************************************/ CPLErr BTDataset::GetGeoTransform( double * padfTransform ) { memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 ); if( bGeoTransformValid ) return CE_None; else return CE_Failure; } /************************************************************************/ /* SetGeoTransform() */ /************************************************************************/ CPLErr BTDataset::SetGeoTransform( double *padfTransform ) { CPLErr eErr = CE_None; memcpy( adfGeoTransform, padfTransform, sizeof(double) * 6 ); if( adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0 ) { CPLError( CE_Failure, CPLE_AppDefined, ".bt format does not support rotational coefficients in geotransform, ignoring." ); eErr = CE_Failure; } /* -------------------------------------------------------------------- */ /* Compute bounds, and update header info. */ /* -------------------------------------------------------------------- */ double dfLeft, dfRight, dfTop, dfBottom; dfLeft = adfGeoTransform[0]; dfRight = dfLeft + adfGeoTransform[1] * nRasterXSize; dfTop = adfGeoTransform[3]; dfBottom = dfTop + adfGeoTransform[5] * nRasterYSize; memcpy( abyHeader + 28, &dfLeft, 8 ); memcpy( abyHeader + 36, &dfRight, 8 ); memcpy( abyHeader + 44, &dfBottom, 8 ); memcpy( abyHeader + 52, &dfTop, 8 ); CPL_LSBPTR64( abyHeader + 28 ); CPL_LSBPTR64( abyHeader + 36 ); CPL_LSBPTR64( abyHeader + 44 ); CPL_LSBPTR64( abyHeader + 52 ); bHeaderModified = TRUE; return eErr; } /************************************************************************/ /* GetProjectionRef() */ /************************************************************************/ const char *BTDataset::GetProjectionRef() { if( pszProjection == NULL ) return ""; else return pszProjection; } /************************************************************************/ /* SetProjection() */ /************************************************************************/ CPLErr BTDataset::SetProjection( const char *pszNewProjection ) { CPLErr eErr = CE_None; CPLFree( pszProjection ); pszProjection = CPLStrdup( pszNewProjection ); bHeaderModified = TRUE; /* -------------------------------------------------------------------- */ /* Parse projection. */ /* -------------------------------------------------------------------- */ OGRSpatialReference oSRS( pszProjection ); GInt16 nShortTemp; /* -------------------------------------------------------------------- */ /* Linear units. */ /* -------------------------------------------------------------------- */ if( oSRS.IsGeographic() ) nShortTemp = 0; else { double dfLinear = oSRS.GetLinearUnits(); if( ABS(dfLinear - 0.3048) < 0.0000001 ) nShortTemp = 2; else if( ABS(dfLinear - atof(SRS_UL_US_FOOT_CONV)) < 0.00000001 ) nShortTemp = 3; else nShortTemp = 1; } nShortTemp = CPL_LSBWORD16( 1 ); memcpy( abyHeader + 22, &nShortTemp, 2 ); /* -------------------------------------------------------------------- */ /* UTM Zone */ /* -------------------------------------------------------------------- */ int bNorth; nShortTemp = oSRS.GetUTMZone( &bNorth ); if( bNorth ) nShortTemp = -nShortTemp; nShortTemp = CPL_LSBWORD16( nShortTemp ); memcpy( abyHeader + 24, &nShortTemp, 2 ); /* -------------------------------------------------------------------- */ /* Datum */ /* -------------------------------------------------------------------- */ if( oSRS.GetAuthorityName( "GEOGCS|DATUM" ) != NULL && EQUAL(oSRS.GetAuthorityName( "GEOGCS|DATUM" ),"EPSG") ) nShortTemp = atoi(oSRS.GetAuthorityName( "GEOGCS|DATUM" )) + 2000; else nShortTemp = -2; nShortTemp = CPL_LSBWORD16( nShortTemp ); /* datum unknown */ memcpy( abyHeader + 26, &nShortTemp, 2 ); /* -------------------------------------------------------------------- */ /* Write out the projection to a .prj file. */ /* -------------------------------------------------------------------- */ const char *pszPrjFile = CPLResetExtension( GetDescription(), "prj" ); FILE * fp; fp = VSIFOpen( pszPrjFile, "wt" ); if( fp != NULL ) { VSIFPrintf( fp, "%s\n", pszProjection ); VSIFClose( fp ); abyHeader[60] = 1; } else { CPLError( CE_Failure, CPLE_AppDefined, "Unable to write out .prj file." ); eErr = CE_Failure; } return eErr; } /************************************************************************/ /* Open() */ /************************************************************************/ GDALDataset *BTDataset::Open( GDALOpenInfo * poOpenInfo ) { /* -------------------------------------------------------------------- */ /* Verify that this is some form of binterr file. */ /* -------------------------------------------------------------------- */ if( poOpenInfo->nHeaderBytes < 256 || poOpenInfo->fp == NULL ) return NULL; if( strncmp( (const char *) poOpenInfo->pabyHeader, "binterr", 7 ) != 0 ) return NULL; /* -------------------------------------------------------------------- */ /* Create the dataset. */ /* -------------------------------------------------------------------- */ BTDataset *poDS; poDS = new BTDataset(); memcpy( poDS->abyHeader, poOpenInfo->pabyHeader, 256 ); /* -------------------------------------------------------------------- */ /* Get the version. */ /* -------------------------------------------------------------------- */ char szVersion[4]; strncpy( szVersion, (char *) (poDS->abyHeader + 7), 3 ); szVersion[3] = '\0'; poDS->nVersionCode = (int) (atof(szVersion) * 10); /* -------------------------------------------------------------------- */ /* Extract core header information, being careful about the */ /* version. */ /* -------------------------------------------------------------------- */ GInt32 nIntTemp; GInt16 nDataSize; GDALDataType eType; memcpy( &nIntTemp, poDS->abyHeader + 10, 4 ); poDS->nRasterXSize = CPL_LSBWORD32( nIntTemp ); memcpy( &nIntTemp, poDS->abyHeader + 14, 4 ); poDS->nRasterYSize = CPL_LSBWORD32( nIntTemp ); memcpy( &nDataSize, poDS->abyHeader+18, 2 ); nDataSize = CPL_LSBWORD16( nDataSize ); if( poDS->abyHeader[20] != 0 && nDataSize == 4 ) eType = GDT_Float32; else if( poDS->abyHeader[20] == 0 && nDataSize == 4 ) eType = GDT_Int32; else if( poDS->abyHeader[20] == 0 && nDataSize == 2 ) eType = GDT_Int16; else { CPLError( CE_Failure, CPLE_AppDefined, ".bt file data type unknown, got datasize=%d.", nDataSize ); return NULL; } /* rcg, apr 7/06: read offset 62 for vert. units. If zero, assume 1.0 as per spec. */ memcpy( &poDS->m_fVscale, poDS->abyHeader + 62, 4 ); CPL_LSBPTR32(&poDS->m_fVscale); if(poDS->m_fVscale == 0.0f) poDS->m_fVscale = 1.0f; /* -------------------------------------------------------------------- */ /* Try to read a .prj file if it is indicated. */ /* -------------------------------------------------------------------- */ OGRSpatialReference oSRS; if( poDS->nVersionCode >= 12 && poDS->abyHeader[60] != 0 ) { const char *pszPrjFile = CPLResetExtension( poOpenInfo->pszFilename, "prj" ); FILE *fp; fp = VSIFOpen( pszPrjFile, "rt" ); if( fp != NULL ) { char *pszBuffer, *pszBufPtr; int nBufMax = 100000; int nBytes; pszBuffer = (char *) CPLMalloc(nBufMax); nBytes = VSIFRead( pszBuffer, 1, nBufMax-1, fp ); VSIFClose( fp ); pszBuffer[nBytes] = '\0'; pszBufPtr = pszBuffer; if( oSRS.importFromWkt( &pszBufPtr ) != OGRERR_NONE ) { CPLError( CE_Warning, CPLE_AppDefined, "Unable to parse .prj file, coordinate system missing." ); } CPLFree( pszBuffer ); } } /* -------------------------------------------------------------------- */ /* If we didn't find a .prj file, try to use internal info. */ /* -------------------------------------------------------------------- */ if( oSRS.GetRoot() == NULL ) { GInt16 nUTMZone, nDatum, nHUnits; memcpy( &nUTMZone, poDS->abyHeader + 24, 2 ); nUTMZone = CPL_LSBWORD16( nUTMZone ); memcpy( &nDatum, poDS->abyHeader + 26, 2 ); nDatum = CPL_LSBWORD16( nDatum ); memcpy( &nHUnits, poDS->abyHeader + 22, 2 ); nHUnits = CPL_LSBWORD16( nHUnits ); if( nUTMZone != 0 ) oSRS.SetUTM( ABS(nUTMZone), nUTMZone > 0 ); else if( nHUnits != 0 ) oSRS.SetLocalCS( "Unknown" ); if( nHUnits == 1 ) oSRS.SetLinearUnits( SRS_UL_METER, 1.0 ); else if( nHUnits == 2 ) oSRS.SetLinearUnits( SRS_UL_FOOT, atof(SRS_UL_FOOT_CONV) ); else if( nHUnits == 3 ) oSRS.SetLinearUnits( SRS_UL_US_FOOT, atof(SRS_UL_US_FOOT_CONV) ); // Translate some of the more obvious old USGS datum codes if( nDatum == 0 ) nDatum = 6201; else if( nDatum == 1 ) nDatum = 6209; else if( nDatum == 2 ) nDatum = 6210; else if( nDatum == 3 ) nDatum = 6202; else if( nDatum == 4 ) nDatum = 6203; else if( nDatum == 4 ) nDatum = 6203; else if( nDatum == 6 ) nDatum = 6222; else if( nDatum == 7 ) nDatum = 6230; else if( nDatum == 13 ) nDatum = 6267; else if( nDatum == 14 ) nDatum = 6269; else if( nDatum == 17 ) nDatum = 6277; else if( nDatum == 19 ) nDatum = 6284; else if( nDatum == 21 ) nDatum = 6301; else if( nDatum == 22 ) nDatum = 6322; else if( nDatum == 23 ) nDatum = 6326; if( !oSRS.IsLocal() ) { if( nDatum >= 6000 ) { char szName[32]; sprintf( szName, "EPSG:%d", nDatum-2000 ); oSRS.SetWellKnownGeogCS( szName ); } else oSRS.SetWellKnownGeogCS( "WGS84" ); } } /* -------------------------------------------------------------------- */ /* Convert coordinate system back to WKT. */ /* -------------------------------------------------------------------- */ if( oSRS.GetRoot() != NULL ) oSRS.exportToWkt( &poDS->pszProjection ); /* -------------------------------------------------------------------- */ /* Get georeferencing bounds. */ /* -------------------------------------------------------------------- */ if( poDS->nVersionCode >= 11 ) { double dfLeft, dfRight, dfTop, dfBottom; memcpy( &dfLeft, poDS->abyHeader + 28, 8 ); CPL_LSBPTR64( &dfLeft ); memcpy( &dfRight, poDS->abyHeader + 36, 8 ); CPL_LSBPTR64( &dfRight ); memcpy( &dfBottom, poDS->abyHeader + 44, 8 ); CPL_LSBPTR64( &dfBottom ); memcpy( &dfTop, poDS->abyHeader + 52, 8 ); CPL_LSBPTR64( &dfTop ); poDS->adfGeoTransform[0] = dfLeft; poDS->adfGeoTransform[1] = (dfRight - dfLeft) / poDS->nRasterXSize; poDS->adfGeoTransform[2] = 0.0; poDS->adfGeoTransform[3] = dfTop; poDS->adfGeoTransform[4] = 0.0; poDS->adfGeoTransform[5] = (dfBottom - dfTop) / poDS->nRasterYSize; poDS->bGeoTransformValid = TRUE; } /* -------------------------------------------------------------------- */ /* Re-open the file with the desired access. */ /* -------------------------------------------------------------------- */ VSIFClose( poOpenInfo->fp ); poOpenInfo->fp = NULL; if( poOpenInfo->eAccess == GA_Update ) poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "rb+" ); else poDS->fpImage = VSIFOpenL( poOpenInfo->pszFilename, "rb" ); if( poDS->fpImage == NULL ) { CPLError( CE_Failure, CPLE_OpenFailed, "Failed to re-open %s within BT driver.\n", poOpenInfo->pszFilename ); return NULL; } poDS->eAccess = poOpenInfo->eAccess; /* -------------------------------------------------------------------- */ /* Create band information objects */ /* -------------------------------------------------------------------- */ poDS->SetBand( 1, new BTRasterBand( poDS, poDS->fpImage, eType ) ); /* -------------------------------------------------------------------- */ /* Check for overviews. */ /* -------------------------------------------------------------------- */ poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename ); #ifdef notdef poDS->bGeoTransformValid = GDALReadWorldFile( poOpenInfo->pszFilename, ".wld", poDS->adfGeoTransform ); #endif /* -------------------------------------------------------------------- */ /* Initialize any PAM information. */ /* -------------------------------------------------------------------- */ poDS->SetDescription( poOpenInfo->pszFilename ); poDS->TryLoadXML(); return( poDS ); } /************************************************************************/ /* Create() */ /************************************************************************/ GDALDataset *BTDataset::Create( const char * pszFilename, int nXSize, int nYSize, int nBands, GDALDataType eType, char ** papszOptions ) { /* -------------------------------------------------------------------- */ /* Verify input options. */ /* -------------------------------------------------------------------- */ if( eType != GDT_Int16 && eType != GDT_Int32 && eType != GDT_Float32 ) { CPLError( CE_Failure, CPLE_AppDefined, "Attempt to create .bt dataset with an illegal\n" "data type (%s), only Int16, Int32 and Float32 supported.\n", GDALGetDataTypeName(eType) ); return NULL; } if( nBands != 1 ) { CPLError( CE_Failure, CPLE_AppDefined, "Attempt to create .bt dataset with %d bands, only 1 supported", nBands ); return NULL; } /* -------------------------------------------------------------------- */ /* Try to create the file. */ /* -------------------------------------------------------------------- */ FILE *fp; fp = VSIFOpenL( pszFilename, "wb" ); if( fp == NULL ) { CPLError( CE_Failure, CPLE_OpenFailed, "Attempt to create file `%s' failed.\n", pszFilename ); return NULL; } /* -------------------------------------------------------------------- */ /* Setup base header. */ /* -------------------------------------------------------------------- */ unsigned char abyHeader[256]; GInt32 nTemp; GInt16 nShortTemp; memset( abyHeader, 0, 256 ); memcpy( abyHeader, "binterr1.3", 10 ); nTemp = CPL_LSBWORD32( nXSize ); memcpy( abyHeader+10, &nTemp, 4 ); nTemp = CPL_LSBWORD32( nYSize ); memcpy( abyHeader+14, &nTemp, 4 ); nShortTemp = CPL_LSBWORD16( GDALGetDataTypeSize( eType ) / 8 ); memcpy( abyHeader + 18, &nShortTemp, 2 ); if( eType == GDT_Float32 ) abyHeader[20] = 1; else abyHeader[20] = 0; nShortTemp = CPL_LSBWORD16( 1 ); /* meters */ memcpy( abyHeader + 22, &nShortTemp, 2 ); nShortTemp = CPL_LSBWORD16( 0 ); /* not utm */ memcpy( abyHeader + 24, &nShortTemp, 2 ); nShortTemp = CPL_LSBWORD16( -2 ); /* datum unknown */ memcpy( abyHeader + 26, &nShortTemp, 2 ); /* -------------------------------------------------------------------- */ /* Set dummy extents. */ /* -------------------------------------------------------------------- */ double dfLeft=0, dfRight=nXSize, dfTop=nYSize, dfBottom=0; memcpy( abyHeader + 28, &dfLeft, 8 ); memcpy( abyHeader + 36, &dfRight, 8 ); memcpy( abyHeader + 44, &dfBottom, 8 ); memcpy( abyHeader + 52, &dfTop, 8 ); CPL_LSBPTR64( abyHeader + 28 ); CPL_LSBPTR64( abyHeader + 36 ); CPL_LSBPTR64( abyHeader + 44 ); CPL_LSBPTR64( abyHeader + 52 ); /* -------------------------------------------------------------------- */ /* Set dummy scale. */ /* -------------------------------------------------------------------- */ float fScale = 1.0; memcpy( abyHeader + 62, &fScale, 4 ); CPL_LSBPTR32( abyHeader + 62 ); /* -------------------------------------------------------------------- */ /* Write to disk. */ /* -------------------------------------------------------------------- */ VSIFWriteL( (void *) abyHeader, 256, 1, fp ); if( VSIFSeekL( fp, (GDALGetDataTypeSize(eType)/8) * nXSize * nYSize - 1, SEEK_CUR ) != 0 || VSIFWriteL( abyHeader+255, 1, 1, fp ) != 1 ) { CPLError( CE_Failure, CPLE_FileIO, "Failed to extent file to its full size, out of disk space?" ); VSIFCloseL( fp ); VSIUnlink( pszFilename ); return NULL; } VSIFCloseL( fp ); return (GDALDataset *) GDALOpen( pszFilename, GA_Update ); } /************************************************************************/ /* GDALRegister_BT() */ /************************************************************************/ void GDALRegister_BT() { GDALDriver *poDriver; if( GDALGetDriverByName( "BT" ) == NULL ) { poDriver = new GDALDriver(); poDriver->SetDescription( "BT" ); poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "VTP .bt (Binary Terrain) 1.3 Format" ); poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "frmt_various.html#BT" ); poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "bt" ); poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, "Int16 Int32 Float32" ); poDriver->pfnOpen = BTDataset::Open; poDriver->pfnCreate = BTDataset::Create; GetGDALDriverManager()->RegisterDriver( poDriver ); } }