1 : /******************************************************************************
2 : * $Id: doq1dataset.cpp 21717 2011-02-13 20:16:30Z rouault $
3 : *
4 : * Project: USGS DOQ Driver (First Generation Format)
5 : * Purpose: Implementation of DOQ1Dataset
6 : * Author: Frank Warmerdam, warmerda@home.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 1999, Frank Warmerdam
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "rawdataset.h"
31 : #include "cpl_string.h"
32 :
33 : CPL_CVSID("$Id: doq1dataset.cpp 21717 2011-02-13 20:16:30Z rouault $");
34 :
35 : static double DOQGetField( unsigned char *, int );
36 : static void DOQGetDescription( GDALDataset *, unsigned char * );
37 :
38 : CPL_C_START
39 : void GDALRegister_DOQ1(void);
40 : CPL_C_END
41 :
42 : #define UTM_FORMAT \
43 : "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]"
44 :
45 : #define WGS84_DATUM \
46 : "\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]]"
47 :
48 : #define WGS72_DATUM \
49 : "\"WGS 72\",DATUM[\"WGS_1972\",SPHEROID[\"NWL 10D\",6378135,298.26]]"
50 :
51 : #define NAD27_DATUM \
52 : "\"NAD27\",DATUM[\"North_American_Datum_1927\",SPHEROID[\"Clarke 1866\",6378206.4,294.978698213901]]"
53 :
54 : #define NAD83_DATUM \
55 : "\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.257222101]]"
56 :
57 : /************************************************************************/
58 : /* ==================================================================== */
59 : /* DOQ1Dataset */
60 : /* ==================================================================== */
61 : /************************************************************************/
62 :
63 : class DOQ1Dataset : public RawDataset
64 : {
65 : VSILFILE *fpImage; // image data file.
66 :
67 : double dfULX, dfULY;
68 : double dfXPixelSize, dfYPixelSize;
69 :
70 : char *pszProjection;
71 :
72 : public:
73 : DOQ1Dataset();
74 : ~DOQ1Dataset();
75 :
76 : CPLErr GetGeoTransform( double * padfTransform );
77 : const char *GetProjectionRef( void );
78 :
79 : static GDALDataset *Open( GDALOpenInfo * );
80 : };
81 :
82 : /************************************************************************/
83 : /* DOQ1Dataset() */
84 : /************************************************************************/
85 :
86 1 : DOQ1Dataset::DOQ1Dataset()
87 : {
88 1 : pszProjection = NULL;
89 1 : fpImage = NULL;
90 1 : }
91 :
92 : /************************************************************************/
93 : /* ~DOQ1Dataset() */
94 : /************************************************************************/
95 :
96 1 : DOQ1Dataset::~DOQ1Dataset()
97 :
98 : {
99 1 : FlushCache();
100 :
101 1 : CPLFree( pszProjection );
102 1 : if( fpImage != NULL )
103 1 : VSIFCloseL( fpImage );
104 1 : }
105 :
106 : /************************************************************************/
107 : /* GetGeoTransform() */
108 : /************************************************************************/
109 :
110 0 : CPLErr DOQ1Dataset::GetGeoTransform( double * padfTransform )
111 :
112 : {
113 0 : padfTransform[0] = dfULX;
114 0 : padfTransform[1] = dfXPixelSize;
115 0 : padfTransform[2] = 0.0;
116 0 : padfTransform[3] = dfULY;
117 0 : padfTransform[4] = 0.0;
118 0 : padfTransform[5] = -1 * dfYPixelSize;
119 :
120 0 : return( CE_None );
121 : }
122 :
123 : /************************************************************************/
124 : /* GetProjectionString() */
125 : /************************************************************************/
126 :
127 0 : const char *DOQ1Dataset::GetProjectionRef()
128 :
129 : {
130 0 : return pszProjection;
131 : }
132 :
133 : /************************************************************************/
134 : /* Open() */
135 : /************************************************************************/
136 :
137 12340 : GDALDataset *DOQ1Dataset::Open( GDALOpenInfo * poOpenInfo )
138 :
139 : {
140 : int nWidth, nHeight, nBandStorage, nBandTypes;
141 :
142 : /* -------------------------------------------------------------------- */
143 : /* We assume the user is pointing to the binary (ie. .bil) file. */
144 : /* -------------------------------------------------------------------- */
145 12340 : if( poOpenInfo->nHeaderBytes < 212 )
146 11652 : return NULL;
147 :
148 : /* -------------------------------------------------------------------- */
149 : /* Attempt to extract a few key values from the header. */
150 : /* -------------------------------------------------------------------- */
151 688 : nWidth = (int) DOQGetField(poOpenInfo->pabyHeader + 150, 6);
152 688 : nHeight = (int) DOQGetField(poOpenInfo->pabyHeader + 144, 6);
153 688 : nBandStorage = (int) DOQGetField(poOpenInfo->pabyHeader + 162, 3);
154 688 : nBandTypes = (int) DOQGetField(poOpenInfo->pabyHeader + 156, 3);
155 :
156 : /* -------------------------------------------------------------------- */
157 : /* Do these values look coherent for a DOQ file? It would be */
158 : /* nice to do a more comprehensive test than this! */
159 : /* -------------------------------------------------------------------- */
160 688 : if( nWidth < 500 || nWidth > 25000
161 : || nHeight < 500 || nHeight > 25000
162 : || nBandStorage < 0 || nBandStorage > 4
163 : || nBandTypes < 1 || nBandTypes > 9 )
164 687 : return NULL;
165 :
166 : /* -------------------------------------------------------------------- */
167 : /* Check the configuration. We don't currently handle all */
168 : /* variations, only the common ones. */
169 : /* -------------------------------------------------------------------- */
170 1 : if( nBandTypes > 5 )
171 : {
172 : CPLError( CE_Failure, CPLE_OpenFailed,
173 : "DOQ Data Type (%d) is not a supported configuration.\n",
174 0 : nBandTypes );
175 0 : return NULL;
176 : }
177 :
178 : /* -------------------------------------------------------------------- */
179 : /* Confirm the requested access is supported. */
180 : /* -------------------------------------------------------------------- */
181 1 : if( poOpenInfo->eAccess == GA_Update )
182 : {
183 : CPLError( CE_Failure, CPLE_NotSupported,
184 : "The DOQ1 driver does not support update access to existing"
185 0 : " datasets.\n" );
186 0 : return NULL;
187 : }
188 :
189 : /* -------------------------------------------------------------------- */
190 : /* Create a corresponding GDALDataset. */
191 : /* -------------------------------------------------------------------- */
192 : DOQ1Dataset *poDS;
193 :
194 1 : poDS = new DOQ1Dataset();
195 :
196 : /* -------------------------------------------------------------------- */
197 : /* Capture some information from the file that is of interest. */
198 : /* -------------------------------------------------------------------- */
199 1 : poDS->nRasterXSize = nWidth;
200 1 : poDS->nRasterYSize = nHeight;
201 :
202 1 : poDS->fpImage = VSIFOpenL(poOpenInfo->pszFilename, "rb");
203 1 : if (poDS->fpImage == NULL)
204 : {
205 0 : delete poDS;
206 0 : return NULL;
207 : }
208 :
209 : /* -------------------------------------------------------------------- */
210 : /* Compute layout of data. */
211 : /* -------------------------------------------------------------------- */
212 1 : int nSkipBytes, nBytesPerPixel=0, nBytesPerLine, i;
213 :
214 1 : if( nBandTypes < 5 )
215 1 : nBytesPerPixel = 1;
216 0 : else if( nBandTypes == 5 )
217 0 : nBytesPerPixel = 3;
218 :
219 1 : nBytesPerLine = nBytesPerPixel * nWidth;
220 1 : nSkipBytes = 4 * nBytesPerLine;
221 :
222 : /* -------------------------------------------------------------------- */
223 : /* Create band information objects. */
224 : /* -------------------------------------------------------------------- */
225 1 : poDS->nBands = nBytesPerPixel;
226 4 : for( i = 0; i < poDS->nBands; i++ )
227 : {
228 : poDS->SetBand( i+1,
229 : new RawRasterBand( poDS, i+1, poDS->fpImage,
230 : nSkipBytes + i, nBytesPerPixel, nBytesPerLine,
231 1 : GDT_Byte, TRUE, TRUE ) );
232 : }
233 :
234 : /* -------------------------------------------------------------------- */
235 : /* Set the description. */
236 : /* -------------------------------------------------------------------- */
237 1 : DOQGetDescription(poDS, poOpenInfo->pabyHeader);
238 :
239 : /* -------------------------------------------------------------------- */
240 : /* Establish the projection string. */
241 : /* -------------------------------------------------------------------- */
242 1 : if( ((int) DOQGetField(poOpenInfo->pabyHeader + 195, 3)) != 1 )
243 0 : poDS->pszProjection = VSIStrdup("");
244 : else
245 : {
246 : const char *pszDatumLong, *pszDatumShort;
247 : const char *pszUnits;
248 : int nZone;
249 :
250 1 : nZone = (int) DOQGetField(poOpenInfo->pabyHeader + 198, 6);
251 :
252 1 : if( ((int) DOQGetField(poOpenInfo->pabyHeader + 204, 3)) == 1 )
253 0 : pszUnits = "UNIT[\"US survey foot\",0.304800609601219]";
254 : else
255 1 : pszUnits = "UNIT[\"metre\",1]";
256 :
257 1 : switch( (int) DOQGetField(poOpenInfo->pabyHeader + 167, 2) )
258 : {
259 : case 1:
260 0 : pszDatumLong = NAD27_DATUM;
261 0 : pszDatumShort = "NAD 27";
262 0 : break;
263 :
264 : case 2:
265 0 : pszDatumLong = WGS72_DATUM;
266 0 : pszDatumShort = "WGS 72";
267 0 : break;
268 :
269 : case 3:
270 1 : pszDatumLong = WGS84_DATUM;
271 1 : pszDatumShort = "WGS 84";
272 1 : break;
273 :
274 : case 4:
275 0 : pszDatumLong = NAD83_DATUM;
276 0 : pszDatumShort = "NAD 83";
277 0 : break;
278 :
279 : default:
280 0 : pszDatumLong = "DATUM[\"unknown\"]";
281 0 : pszDatumShort = "unknown";
282 : break;
283 : }
284 :
285 : poDS->pszProjection =
286 : CPLStrdup(CPLSPrintf( UTM_FORMAT, pszDatumShort, nZone,
287 1 : pszDatumLong, nZone * 6 - 183, pszUnits ));
288 : }
289 :
290 : /* -------------------------------------------------------------------- */
291 : /* Read the georeferencing information. */
292 : /* -------------------------------------------------------------------- */
293 : unsigned char abyRecordData[500];
294 :
295 1 : if( VSIFSeekL( poDS->fpImage, nBytesPerLine * 2, SEEK_SET ) != 0
296 : || VSIFReadL(abyRecordData,sizeof(abyRecordData),1,poDS->fpImage) != 1 )
297 : {
298 : CPLError( CE_Failure, CPLE_FileIO,
299 : "Header read error on %s.\n",
300 0 : poOpenInfo->pszFilename );
301 0 : delete poDS;
302 0 : return NULL;
303 : }
304 :
305 1 : poDS->dfULX = DOQGetField( abyRecordData + 288, 24 );
306 1 : poDS->dfULY = DOQGetField( abyRecordData + 312, 24 );
307 :
308 1 : if( VSIFSeekL( poDS->fpImage, nBytesPerLine * 3, SEEK_SET ) != 0
309 : || VSIFReadL(abyRecordData,sizeof(abyRecordData),1,poDS->fpImage) != 1 )
310 : {
311 : CPLError( CE_Failure, CPLE_FileIO,
312 : "Header read error on %s.\n",
313 0 : poOpenInfo->pszFilename );
314 0 : delete poDS;
315 0 : return NULL;
316 : }
317 :
318 1 : poDS->dfXPixelSize = DOQGetField( abyRecordData + 59, 12 );
319 1 : poDS->dfYPixelSize = DOQGetField( abyRecordData + 71, 12 );
320 :
321 : /* -------------------------------------------------------------------- */
322 : /* Initialize any PAM information. */
323 : /* -------------------------------------------------------------------- */
324 1 : poDS->SetDescription( poOpenInfo->pszFilename );
325 1 : poDS->TryLoadXML();
326 :
327 : /* -------------------------------------------------------------------- */
328 : /* Check for overviews. */
329 : /* -------------------------------------------------------------------- */
330 1 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
331 :
332 1 : return( poDS );
333 : }
334 :
335 : /************************************************************************/
336 : /* GDALRegister_DOQ1() */
337 : /************************************************************************/
338 :
339 582 : void GDALRegister_DOQ1()
340 :
341 : {
342 : GDALDriver *poDriver;
343 :
344 582 : if( GDALGetDriverByName( "DOQ1" ) == NULL )
345 : {
346 561 : poDriver = new GDALDriver();
347 :
348 561 : poDriver->SetDescription( "DOQ1" );
349 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
350 561 : "USGS DOQ (Old Style)" );
351 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
352 561 : "frmt_various.html#DOQ1" );
353 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
354 :
355 561 : poDriver->pfnOpen = DOQ1Dataset::Open;
356 :
357 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
358 : }
359 582 : }
360 :
361 : /************************************************************************/
362 : /* DOQGetField() */
363 : /************************************************************************/
364 :
365 2760 : static double DOQGetField( unsigned char *pabyData, int nBytes )
366 :
367 : {
368 : char szWork[128];
369 : int i;
370 :
371 2760 : strncpy( szWork, (const char *) pabyData, nBytes );
372 2760 : szWork[nBytes] = '\0';
373 :
374 15230 : for( i = 0; i < nBytes; i++ )
375 : {
376 12470 : if( szWork[i] == 'D' || szWork[i] == 'd' )
377 59 : szWork[i] = 'E';
378 : }
379 :
380 2760 : return atof(szWork);
381 : }
382 :
383 : /************************************************************************/
384 : /* DOQGetDescription() */
385 : /************************************************************************/
386 :
387 1 : static void DOQGetDescription( GDALDataset *poDS, unsigned char *pabyData )
388 :
389 : {
390 : char szWork[128];
391 1 : int i = 0;
392 :
393 1 : memset( szWork, ' ', 128 );
394 1 : strncpy( szWork, "USGS GeoTIFF DOQ 1:12000 Q-Quad of ", 35 );
395 1 : strncpy( szWork + 35, (const char *) pabyData + 0, 38 );
396 2 : while ( *(szWork + 72 - i) == ' ' ) {
397 0 : i++;
398 : }
399 1 : i--;
400 1 : strncpy( szWork + 73 - i, (const char *) pabyData + 38, 2 );
401 1 : strncpy( szWork + 76 - i, (const char *) pabyData + 44, 2 );
402 1 : szWork[77-i] = '\0';
403 :
404 1 : poDS->SetMetadataItem( "DOQ_DESC", szWork );
405 1 : }
|
