1 : /******************************************************************************
2 : * $Id: doq1dataset.cpp 17664 2009-09-21 21:16:45Z 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 17664 2009-09-21 21:16:45Z 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 : FILE *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 2 : DOQ1Dataset::~DOQ1Dataset()
97 :
98 : {
99 1 : FlushCache();
100 :
101 1 : CPLFree( pszProjection );
102 1 : if( fpImage != NULL )
103 1 : VSIFClose( fpImage );
104 2 : }
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 8736 : 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 8736 : if( poOpenInfo->nHeaderBytes < 212 || poOpenInfo->fp == NULL )
146 8485 : return NULL;
147 :
148 : /* -------------------------------------------------------------------- */
149 : /* Attempt to extract a few key values from the header. */
150 : /* -------------------------------------------------------------------- */
151 251 : nWidth = (int) DOQGetField(poOpenInfo->pabyHeader + 150, 6);
152 251 : nHeight = (int) DOQGetField(poOpenInfo->pabyHeader + 144, 6);
153 251 : nBandStorage = (int) DOQGetField(poOpenInfo->pabyHeader + 162, 3);
154 251 : 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 251 : if( nWidth < 500 || nWidth > 25000
161 : || nHeight < 500 || nHeight > 25000
162 : || nBandStorage < 0 || nBandStorage > 4
163 : || nBandTypes < 1 || nBandTypes > 9 )
164 250 : 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 : /* -------------------------------------------------------------------- */
203 : /* Assume ownership of the file handled from the GDALOpenInfo. */
204 : /* -------------------------------------------------------------------- */
205 1 : poDS->fpImage = poOpenInfo->fp;
206 1 : poOpenInfo->fp = NULL;
207 :
208 : /* -------------------------------------------------------------------- */
209 : /* Compute layout of data. */
210 : /* -------------------------------------------------------------------- */
211 1 : int nSkipBytes, nBytesPerPixel=0, nBytesPerLine, i;
212 :
213 1 : if( nBandTypes < 5 )
214 1 : nBytesPerPixel = 1;
215 0 : else if( nBandTypes == 5 )
216 0 : nBytesPerPixel = 3;
217 :
218 1 : nBytesPerLine = nBytesPerPixel * nWidth;
219 1 : nSkipBytes = 4 * nBytesPerLine;
220 :
221 : /* -------------------------------------------------------------------- */
222 : /* Create band information objects. */
223 : /* -------------------------------------------------------------------- */
224 1 : poDS->nBands = nBytesPerPixel;
225 4 : for( i = 0; i < poDS->nBands; i++ )
226 : {
227 : poDS->SetBand( i+1,
228 : new RawRasterBand( poDS, i+1, poDS->fpImage,
229 : nSkipBytes + i, nBytesPerPixel, nBytesPerLine,
230 1 : GDT_Byte, TRUE ) );
231 : }
232 :
233 : /* -------------------------------------------------------------------- */
234 : /* Set the description. */
235 : /* -------------------------------------------------------------------- */
236 1 : DOQGetDescription(poDS, poOpenInfo->pabyHeader);
237 :
238 : /* -------------------------------------------------------------------- */
239 : /* Establish the projection string. */
240 : /* -------------------------------------------------------------------- */
241 1 : if( ((int) DOQGetField(poOpenInfo->pabyHeader + 195, 3)) != 1 )
242 0 : poDS->pszProjection = VSIStrdup("");
243 : else
244 : {
245 : const char *pszDatumLong, *pszDatumShort;
246 : const char *pszUnits;
247 : int nZone;
248 :
249 1 : nZone = (int) DOQGetField(poOpenInfo->pabyHeader + 198, 6);
250 :
251 1 : if( ((int) DOQGetField(poOpenInfo->pabyHeader + 204, 3)) == 1 )
252 0 : pszUnits = "UNIT[\"US survey foot\",0.304800609601219]";
253 : else
254 1 : pszUnits = "UNIT[\"metre\",1]";
255 :
256 1 : switch( (int) DOQGetField(poOpenInfo->pabyHeader + 167, 2) )
257 : {
258 : case 1:
259 0 : pszDatumLong = NAD27_DATUM;
260 0 : pszDatumShort = "NAD 27";
261 0 : break;
262 :
263 : case 2:
264 0 : pszDatumLong = WGS72_DATUM;
265 0 : pszDatumShort = "WGS 72";
266 0 : break;
267 :
268 : case 3:
269 1 : pszDatumLong = WGS84_DATUM;
270 1 : pszDatumShort = "WGS 84";
271 1 : break;
272 :
273 : case 4:
274 0 : pszDatumLong = NAD83_DATUM;
275 0 : pszDatumShort = "NAD 83";
276 0 : break;
277 :
278 : default:
279 0 : pszDatumLong = "DATUM[\"unknown\"]";
280 0 : pszDatumShort = "unknown";
281 : break;
282 : }
283 :
284 : poDS->pszProjection =
285 : CPLStrdup(CPLSPrintf( UTM_FORMAT, pszDatumShort, nZone,
286 1 : pszDatumLong, nZone * 6 - 183, pszUnits ));
287 : }
288 :
289 : /* -------------------------------------------------------------------- */
290 : /* Read the georeferencing information. */
291 : /* -------------------------------------------------------------------- */
292 : unsigned char abyRecordData[500];
293 :
294 1 : if( VSIFSeek( poDS->fpImage, nBytesPerLine * 2, SEEK_SET ) != 0
295 : || VSIFRead(abyRecordData,sizeof(abyRecordData),1,poDS->fpImage) != 1 )
296 : {
297 : CPLError( CE_Failure, CPLE_FileIO,
298 : "Header read error on %s.\n",
299 0 : poOpenInfo->pszFilename );
300 0 : delete poDS;
301 0 : return NULL;
302 : }
303 :
304 1 : poDS->dfULX = DOQGetField( abyRecordData + 288, 24 );
305 1 : poDS->dfULY = DOQGetField( abyRecordData + 312, 24 );
306 :
307 1 : if( VSIFSeek( poDS->fpImage, nBytesPerLine * 3, SEEK_SET ) != 0
308 : || VSIFRead(abyRecordData,sizeof(abyRecordData),1,poDS->fpImage) != 1 )
309 : {
310 : CPLError( CE_Failure, CPLE_FileIO,
311 : "Header read error on %s.\n",
312 0 : poOpenInfo->pszFilename );
313 0 : delete poDS;
314 0 : return NULL;
315 : }
316 :
317 1 : poDS->dfXPixelSize = DOQGetField( abyRecordData + 59, 12 );
318 1 : poDS->dfYPixelSize = DOQGetField( abyRecordData + 71, 12 );
319 :
320 : /* -------------------------------------------------------------------- */
321 : /* Initialize any PAM information. */
322 : /* -------------------------------------------------------------------- */
323 1 : poDS->SetDescription( poOpenInfo->pszFilename );
324 1 : poDS->TryLoadXML();
325 :
326 : /* -------------------------------------------------------------------- */
327 : /* Check for overviews. */
328 : /* -------------------------------------------------------------------- */
329 1 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
330 :
331 1 : return( poDS );
332 : }
333 :
334 : /************************************************************************/
335 : /* GDALRegister_DOQ1() */
336 : /************************************************************************/
337 :
338 338 : void GDALRegister_DOQ1()
339 :
340 : {
341 : GDALDriver *poDriver;
342 :
343 338 : if( GDALGetDriverByName( "DOQ1" ) == NULL )
344 : {
345 336 : poDriver = new GDALDriver();
346 :
347 336 : poDriver->SetDescription( "DOQ1" );
348 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
349 336 : "USGS DOQ (Old Style)" );
350 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
351 336 : "frmt_various.html#DOQ1" );
352 :
353 336 : poDriver->pfnOpen = DOQ1Dataset::Open;
354 :
355 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
356 : }
357 338 : }
358 :
359 : /************************************************************************/
360 : /* DOQGetField() */
361 : /************************************************************************/
362 :
363 1012 : static double DOQGetField( unsigned char *pabyData, int nBytes )
364 :
365 : {
366 : char szWork[128];
367 : int i;
368 :
369 1012 : strncpy( szWork, (const char *) pabyData, nBytes );
370 1012 : szWork[nBytes] = '\0';
371 :
372 5616 : for( i = 0; i < nBytes; i++ )
373 : {
374 4604 : if( szWork[i] == 'D' || szWork[i] == 'd' )
375 19 : szWork[i] = 'E';
376 : }
377 :
378 1012 : return atof(szWork);
379 : }
380 :
381 : /************************************************************************/
382 : /* DOQGetDescription() */
383 : /************************************************************************/
384 :
385 1 : static void DOQGetDescription( GDALDataset *poDS, unsigned char *pabyData )
386 :
387 : {
388 : char szWork[128];
389 1 : int i = 0;
390 :
391 1 : memset( szWork, ' ', 128 );
392 1 : strncpy( szWork, "USGS GeoTIFF DOQ 1:12000 Q-Quad of ", 35 );
393 1 : strncpy( szWork + 35, (const char *) pabyData + 0, 38 );
394 2 : while ( *(szWork + 72 - i) == ' ' ) {
395 0 : i++;
396 : }
397 1 : i--;
398 1 : strncpy( szWork + 73 - i, (const char *) pabyData + 38, 2 );
399 1 : strncpy( szWork + 76 - i, (const char *) pabyData + 44, 2 );
400 1 : szWork[77-i] = '\0';
401 :
402 1 : poDS->SetMetadataItem( "DOQ_DESC", szWork );
403 1 : }
|
