/****************************************************************************** * $Id: doq1dataset.cpp,v 1.15 2005/07/25 18:07:48 fwarmerdam Exp $ * * Project: USGS DOQ Driver (First Generation Format) * Purpose: Implementation of DOQ1Dataset * Author: Frank Warmerdam, warmerda@home.com * ****************************************************************************** * Copyright (c) 1999, 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: doq1dataset.cpp,v $ * Revision 1.15 2005/07/25 18:07:48 fwarmerdam * Fixed up description handling. * * Revision 1.14 2005/05/05 13:55:41 fwarmerdam * PAM Enable * * Revision 1.13 2002/09/17 13:09:28 warmerda * fix zone in WKT description - bugzilla 199 * * Revision 1.12 2002/09/16 16:17:53 warmerda * XY_ORIGIN is the top left corner of pixel! * * Revision 1.11 2002/09/04 06:50:37 warmerda * avoid static driver pointers * * Revision 1.10 2002/06/25 14:57:43 warmerda * Added overview support. * * Revision 1.9 2002/06/12 21:12:25 warmerda * update to metadata based driver info * * Revision 1.8 2001/07/18 04:51:57 warmerda * added CPL_CVSID * * Revision 1.7 2000/09/25 21:20:13 warmerda * avoid initialization warnings * * Revision 1.6 2000/08/15 19:28:26 warmerda * added help topic * * Revision 1.5 2000/02/28 16:32:20 warmerda * use SetBand method * * Revision 1.4 2000/01/24 03:08:16 shadow * add description reading * * Revision 1.3 2000/01/06 20:58:21 warmerda * Fixed initial poOpenInfo check. * * Revision 1.2 1999/12/30 02:39:43 warmerda * Added projections support. * * Revision 1.1 1999/12/29 20:43:05 warmerda * New */ #include "rawdataset.h" #include "cpl_string.h" CPL_CVSID("$Id: doq1dataset.cpp,v 1.15 2005/07/25 18:07:48 fwarmerdam Exp $"); static double DOQGetField( unsigned char *, int ); static void DOQGetDescription( GDALDataset *, unsigned char * ); CPL_C_START void GDALRegister_DOQ1(void); CPL_C_END #define UTM_FORMAT \ "PROJCS[\"%s / UTM zone %dN\",GEOGCS[%s,PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%d],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],%s]" #define WGS84_DATUM \ "\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]]" #define WGS72_DATUM \ "\"WGS 72\",DATUM[\"WGS_1972\",SPHEROID[\"NWL 10D\",6378135,298.26]]" #define NAD27_DATUM \ "\"NAD27\",DATUM[\"North_American_Datum_1927\",SPHEROID[\"Clarke 1866\",6378206.4,294.978698213901]]" #define NAD83_DATUM \ "\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.257222101]]" /************************************************************************/ /* ==================================================================== */ /* DOQ1Dataset */ /* ==================================================================== */ /************************************************************************/ class DOQ1Dataset : public RawDataset { FILE *fpImage; // image data file. double dfULX, dfULY; double dfXPixelSize, dfYPixelSize; char *pszProjection; public: DOQ1Dataset(); ~DOQ1Dataset(); CPLErr GetGeoTransform( double * padfTransform ); const char *GetProjectionRef( void ); static GDALDataset *Open( GDALOpenInfo * ); }; /************************************************************************/ /* DOQ1Dataset() */ /************************************************************************/ DOQ1Dataset::DOQ1Dataset() { pszProjection = NULL; fpImage = NULL; } /************************************************************************/ /* ~DOQ1Dataset() */ /************************************************************************/ DOQ1Dataset::~DOQ1Dataset() { FlushCache(); CPLFree( pszProjection ); if( fpImage != NULL ) VSIFClose( fpImage ); } /************************************************************************/ /* GetGeoTransform() */ /************************************************************************/ CPLErr DOQ1Dataset::GetGeoTransform( double * padfTransform ) { padfTransform[0] = dfULX; padfTransform[1] = dfXPixelSize; padfTransform[2] = 0.0; padfTransform[3] = dfULY; padfTransform[4] = 0.0; padfTransform[5] = -1 * dfYPixelSize; return( CE_None ); } /************************************************************************/ /* GetProjectionString() */ /************************************************************************/ const char *DOQ1Dataset::GetProjectionRef() { return pszProjection; } /************************************************************************/ /* Open() */ /************************************************************************/ GDALDataset *DOQ1Dataset::Open( GDALOpenInfo * poOpenInfo ) { int nWidth, nHeight, nBandStorage, nBandTypes; /* -------------------------------------------------------------------- */ /* We assume the user is pointing to the binary (ie. .bil) file. */ /* -------------------------------------------------------------------- */ if( poOpenInfo->nHeaderBytes < 212 || poOpenInfo->fp == NULL ) return NULL; /* -------------------------------------------------------------------- */ /* Attempt to extract a few key values from the header. */ /* -------------------------------------------------------------------- */ nWidth = (int) DOQGetField(poOpenInfo->pabyHeader + 150, 6); nHeight = (int) DOQGetField(poOpenInfo->pabyHeader + 144, 6); nBandStorage = (int) DOQGetField(poOpenInfo->pabyHeader + 162, 3); nBandTypes = (int) DOQGetField(poOpenInfo->pabyHeader + 156, 3); /* -------------------------------------------------------------------- */ /* Do these values look coherent for a DOQ file? It would be */ /* nice to do a more comprehensive test than this! */ /* -------------------------------------------------------------------- */ if( nWidth < 500 || nWidth > 25000 || nHeight < 500 || nHeight > 25000 || nBandStorage < 0 || nBandStorage > 4 || nBandTypes < 1 || nBandTypes > 9 ) return NULL; /* -------------------------------------------------------------------- */ /* Check the configuration. We don't currently handle all */ /* variations, only the common ones. */ /* -------------------------------------------------------------------- */ if( nBandTypes > 5 ) { CPLError( CE_Failure, CPLE_OpenFailed, "DOQ Data Type (%d) is not a supported configuration.\n", nBandTypes ); return NULL; } /* -------------------------------------------------------------------- */ /* Create a corresponding GDALDataset. */ /* -------------------------------------------------------------------- */ DOQ1Dataset *poDS; poDS = new DOQ1Dataset(); /* -------------------------------------------------------------------- */ /* Capture some information from the file that is of interest. */ /* -------------------------------------------------------------------- */ poDS->nRasterXSize = nWidth; poDS->nRasterYSize = nHeight; /* -------------------------------------------------------------------- */ /* Assume ownership of the file handled from the GDALOpenInfo. */ /* -------------------------------------------------------------------- */ poDS->fpImage = poOpenInfo->fp; poOpenInfo->fp = NULL; /* -------------------------------------------------------------------- */ /* Compute layout of data. */ /* -------------------------------------------------------------------- */ int nSkipBytes, nBytesPerPixel=0, nBytesPerLine, i; if( nBandTypes < 5 ) nBytesPerPixel = 1; else if( nBandTypes == 5 ) nBytesPerPixel = 3; nBytesPerLine = nBytesPerPixel * nWidth; nSkipBytes = 4 * nBytesPerLine; /* -------------------------------------------------------------------- */ /* Create band information objects. */ /* -------------------------------------------------------------------- */ poDS->nBands = nBytesPerPixel; for( i = 0; i < poDS->nBands; i++ ) { poDS->SetBand( i+1, new RawRasterBand( poDS, i+1, poDS->fpImage, nSkipBytes + i, nBytesPerPixel, nBytesPerLine, GDT_Byte, TRUE ) ); } /* -------------------------------------------------------------------- */ /* Set the description. */ /* -------------------------------------------------------------------- */ DOQGetDescription(poDS, poOpenInfo->pabyHeader); /* -------------------------------------------------------------------- */ /* Establish the projection string. */ /* -------------------------------------------------------------------- */ if( ((int) DOQGetField(poOpenInfo->pabyHeader + 195, 3)) != 1 ) poDS->pszProjection = VSIStrdup(""); else { const char *pszDatumLong, *pszDatumShort; const char *pszUnits; int nZone; nZone = (int) DOQGetField(poOpenInfo->pabyHeader + 198, 6); if( ((int) DOQGetField(poOpenInfo->pabyHeader + 204, 3)) == 1 ) pszUnits = "UNIT[\"US survey foot\",0.304800609601219]"; else pszUnits = "UNIT[\"metre\",1]"; switch( (int) DOQGetField(poOpenInfo->pabyHeader + 167, 2) ) { case 1: pszDatumLong = NAD27_DATUM; pszDatumShort = "NAD 27"; break; case 2: pszDatumLong = WGS72_DATUM; pszDatumShort = "WGS 72"; break; case 3: pszDatumLong = WGS84_DATUM; pszDatumShort = "WGS 84"; break; case 4: pszDatumLong = NAD83_DATUM; pszDatumShort = "NAD 83"; break; default: pszDatumLong = "DATUM[\"unknown\"]"; pszDatumShort = "unknown"; break; } poDS->pszProjection = CPLStrdup(CPLSPrintf( UTM_FORMAT, pszDatumShort, nZone, pszDatumLong, nZone * 6 - 183, pszUnits )); } /* -------------------------------------------------------------------- */ /* Read the georeferencing information. */ /* -------------------------------------------------------------------- */ unsigned char abyRecordData[500]; if( VSIFSeek( poDS->fpImage, nBytesPerLine * 2, SEEK_SET ) != 0 || VSIFRead(abyRecordData,sizeof(abyRecordData),1,poDS->fpImage) != 1 ) { CPLError( CE_Failure, CPLE_FileIO, "Header read error on %s.\n", poOpenInfo->pszFilename ); return NULL; } poDS->dfULX = DOQGetField( abyRecordData + 288, 24 ); poDS->dfULY = DOQGetField( abyRecordData + 312, 24 ); if( VSIFSeek( poDS->fpImage, nBytesPerLine * 3, SEEK_SET ) != 0 || VSIFRead(abyRecordData,sizeof(abyRecordData),1,poDS->fpImage) != 1 ) { CPLError( CE_Failure, CPLE_FileIO, "Header read error on %s.\n", poOpenInfo->pszFilename ); return NULL; } poDS->dfXPixelSize = DOQGetField( abyRecordData + 59, 12 ); poDS->dfYPixelSize = DOQGetField( abyRecordData + 71, 12 ); /* -------------------------------------------------------------------- */ /* Check for overviews. */ /* -------------------------------------------------------------------- */ poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename ); /* -------------------------------------------------------------------- */ /* Initialize any PAM information. */ /* -------------------------------------------------------------------- */ poDS->SetDescription( poOpenInfo->pszFilename ); poDS->TryLoadXML(); return( poDS ); } /************************************************************************/ /* GDALRegister_DOQ1() */ /************************************************************************/ void GDALRegister_DOQ1() { GDALDriver *poDriver; if( GDALGetDriverByName( "DOQ1" ) == NULL ) { poDriver = new GDALDriver(); poDriver->SetDescription( "DOQ1" ); poDriver->SetMetadataItem( GDAL_DMD_LONGNAME, "USGS DOQ (Old Style)" ); poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "frmt_various.html#DOQ1" ); poDriver->pfnOpen = DOQ1Dataset::Open; GetGDALDriverManager()->RegisterDriver( poDriver ); } } /************************************************************************/ /* DOQGetField() */ /************************************************************************/ static double DOQGetField( unsigned char *pabyData, int nBytes ) { char szWork[128]; int i; strncpy( szWork, (const char *) pabyData, nBytes ); szWork[nBytes] = '\0'; for( i = 0; i < nBytes; i++ ) { if( szWork[i] == 'D' || szWork[i] == 'd' ) szWork[i] = 'E'; } return atof(szWork); } /************************************************************************/ /* DOQGetDescription() */ /************************************************************************/ static void DOQGetDescription( GDALDataset *poDS, unsigned char *pabyData ) { char szWork[128]; int i = 0; memset( szWork, ' ', 128 ); strncpy( szWork, "USGS GeoTIFF DOQ 1:12000 Q-Quad of ", 35 ); strncpy( szWork + 35, (const char *) pabyData + 0, 38 ); while ( *(szWork + 72 - i) == ' ' ) { i++; } i--; strncpy( szWork + 73 - i, (const char *) pabyData + 38, 2 ); strncpy( szWork + 76 - i, (const char *) pabyData + 44, 2 ); szWork[77-i] = '\0'; poDS->SetMetadataItem( "DOQ_DESC", szWork ); }