1 : /******************************************************************************
2 : * $Id: doq2dataset.cpp 21717 2011-02-13 20:16:30Z rouault $
3 : *
4 : * Project: USGS DOQ Driver (Second Generation Format)
5 : * Purpose: Implementation of DOQ2Dataset
6 : * Author: Derrick J Brashear, shadow@dementia.org
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2000, Derrick J Brashear
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: doq2dataset.cpp 21717 2011-02-13 20:16:30Z rouault $");
34 :
35 : CPL_C_START
36 : void GDALRegister_DOQ2(void);
37 : CPL_C_END
38 :
39 : #define UTM_FORMAT \
40 : "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]"
41 :
42 : #define WGS84_DATUM \
43 : "\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]]"
44 :
45 : #define WGS72_DATUM \
46 : "\"WGS 72\",DATUM[\"WGS_1972\",SPHEROID[\"NWL 10D\",6378135,298.26]]"
47 :
48 : #define NAD27_DATUM \
49 : "\"NAD27\",DATUM[\"North_American_Datum_1927\",SPHEROID[\"Clarke 1866\",6378206.4,294.978698213901]]"
50 :
51 : #define NAD83_DATUM \
52 : "\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.257222101]]"
53 :
54 : /************************************************************************/
55 : /* ==================================================================== */
56 : /* DOQ2Dataset */
57 : /* ==================================================================== */
58 : /************************************************************************/
59 :
60 : class DOQ2Dataset : public RawDataset
61 : {
62 : VSILFILE *fpImage; // image data file.
63 :
64 : double dfULX, dfULY;
65 : double dfXPixelSize, dfYPixelSize;
66 :
67 : char *pszProjection;
68 :
69 : public:
70 : DOQ2Dataset();
71 : ~DOQ2Dataset();
72 :
73 : CPLErr GetGeoTransform( double * padfTransform );
74 : const char *GetProjectionRef( void );
75 :
76 : static GDALDataset *Open( GDALOpenInfo * );
77 : };
78 :
79 : /************************************************************************/
80 : /* DOQ2Dataset() */
81 : /************************************************************************/
82 :
83 1 : DOQ2Dataset::DOQ2Dataset()
84 : {
85 1 : pszProjection = NULL;
86 1 : fpImage = NULL;
87 1 : }
88 :
89 : /************************************************************************/
90 : /* ~DOQ2Dataset() */
91 : /************************************************************************/
92 :
93 1 : DOQ2Dataset::~DOQ2Dataset()
94 :
95 : {
96 1 : FlushCache();
97 :
98 1 : CPLFree( pszProjection );
99 1 : if( fpImage != NULL )
100 1 : VSIFCloseL( fpImage );
101 1 : }
102 :
103 : /************************************************************************/
104 : /* GetGeoTransform() */
105 : /************************************************************************/
106 :
107 1 : CPLErr DOQ2Dataset::GetGeoTransform( double * padfTransform )
108 :
109 : {
110 1 : padfTransform[0] = dfULX;
111 1 : padfTransform[1] = dfXPixelSize;
112 1 : padfTransform[2] = 0.0;
113 1 : padfTransform[3] = dfULY;
114 1 : padfTransform[4] = 0.0;
115 1 : padfTransform[5] = -1 * dfYPixelSize;
116 :
117 1 : return( CE_None );
118 : }
119 :
120 : /************************************************************************/
121 : /* GetProjectionString() */
122 : /************************************************************************/
123 :
124 0 : const char *DOQ2Dataset::GetProjectionRef()
125 :
126 : {
127 0 : return pszProjection;
128 : }
129 :
130 : /************************************************************************/
131 : /* Open() */
132 : /************************************************************************/
133 :
134 12339 : GDALDataset *DOQ2Dataset::Open( GDALOpenInfo * poOpenInfo )
135 :
136 : {
137 12339 : int nWidth=0, nHeight=0, nBandStorage=0, nBandTypes=0;
138 :
139 : /* -------------------------------------------------------------------- */
140 : /* We assume the user is pointing to the binary (ie. .bil) file. */
141 : /* -------------------------------------------------------------------- */
142 12339 : if( poOpenInfo->nHeaderBytes < 212 )
143 11652 : return NULL;
144 :
145 687 : int nLineCount = 0;
146 : const char *pszLine;
147 687 : int nBytesPerPixel=0;
148 687 : const char *pszDatumLong=NULL, *pszDatumShort=NULL;
149 687 : const char *pszUnits=NULL;
150 687 : char *pszQuadname = NULL;
151 687 : char *pszQuadquad = NULL;
152 687 : char *pszState = NULL;
153 687 : int nZone=0, nProjType=0;
154 687 : int nSkipBytes=0, nBytesPerLine, i, nBandCount = 0;
155 687 : double dfULXMap=0.0, dfULYMap = 0.0;
156 687 : double dfXDim=0.0, dfYDim=0.0;
157 687 : char **papszMetadata = NULL;
158 :
159 687 : if(! EQUALN((const char *) poOpenInfo->pabyHeader,
160 : "BEGIN_USGS_DOQ_HEADER", 21) )
161 686 : return NULL;
162 :
163 1 : VSILFILE* fp = VSIFOpenL(poOpenInfo->pszFilename, "rb");
164 1 : if (fp == NULL)
165 0 : return NULL;
166 :
167 : /* read and discard the first line */
168 1 : pszLine = CPLReadLineL( fp );
169 :
170 47 : while( (pszLine = CPLReadLineL( fp )) != NULL )
171 : {
172 : char **papszTokens;
173 :
174 46 : nLineCount++;
175 :
176 46 : if( EQUAL(pszLine,"END_USGS_DOQ_HEADER") )
177 0 : break;
178 :
179 46 : papszTokens = CSLTokenizeString( pszLine );
180 46 : if( CSLCount( papszTokens ) < 2 )
181 : {
182 1 : CSLDestroy( papszTokens );
183 1 : break;
184 : }
185 :
186 45 : if( EQUAL(papszTokens[0],"SAMPLES_AND_LINES") && CSLCount(papszTokens) >= 3 )
187 : {
188 1 : nWidth = atoi(papszTokens[1]);
189 1 : nHeight = atoi(papszTokens[2]);
190 : }
191 44 : else if( EQUAL(papszTokens[0],"BYTE_COUNT") )
192 : {
193 1 : nSkipBytes = atoi(papszTokens[1]);
194 : }
195 43 : else if( EQUAL(papszTokens[0],"XY_ORIGIN") && CSLCount(papszTokens) >= 3 )
196 : {
197 1 : dfULXMap = atof(papszTokens[1]);
198 1 : dfULYMap = atof(papszTokens[2]);
199 : }
200 42 : else if( EQUAL(papszTokens[0],"HORIZONTAL_RESOLUTION") )
201 : {
202 1 : dfXDim = dfYDim = atof(papszTokens[1]);
203 : }
204 41 : else if( EQUAL(papszTokens[0],"BAND_ORGANIZATION") )
205 : {
206 1 : if( EQUAL(papszTokens[1],"SINGLE FILE") )
207 0 : nBandStorage = 1;
208 1 : if( EQUAL(papszTokens[1],"BSQ") )
209 0 : nBandStorage = 1;
210 1 : if( EQUAL(papszTokens[1],"BIL") )
211 0 : nBandStorage = 1;
212 1 : if( EQUAL(papszTokens[1],"BIP") )
213 1 : nBandStorage = 4;
214 : }
215 40 : else if( EQUAL(papszTokens[0],"BAND_CONTENT") )
216 : {
217 3 : if( EQUAL(papszTokens[1],"BLACK&WHITE") )
218 0 : nBandTypes = 1;
219 3 : else if( EQUAL(papszTokens[1],"COLOR") )
220 0 : nBandTypes = 5;
221 3 : else if( EQUAL(papszTokens[1],"RGB") )
222 0 : nBandTypes = 5;
223 3 : else if( EQUAL(papszTokens[1],"RED") )
224 1 : nBandTypes = 5;
225 2 : else if( EQUAL(papszTokens[1],"GREEN") )
226 1 : nBandTypes = 5;
227 1 : else if( EQUAL(papszTokens[1],"BLUE") )
228 1 : nBandTypes = 5;
229 :
230 3 : nBandCount++;
231 : }
232 37 : else if( EQUAL(papszTokens[0],"BITS_PER_PIXEL") )
233 : {
234 1 : nBytesPerPixel = (atoi(papszTokens[1]) / 8);
235 : }
236 36 : else if( EQUAL(papszTokens[0],"HORIZONTAL_COORDINATE_SYSTEM") )
237 : {
238 1 : if( EQUAL(papszTokens[1],"UTM") )
239 1 : nProjType = 1;
240 0 : else if( EQUAL(papszTokens[1],"SPCS") )
241 0 : nProjType = 2;
242 0 : else if( EQUAL(papszTokens[1],"GEOGRAPHIC") )
243 0 : nProjType = 0;
244 : }
245 35 : else if( EQUAL(papszTokens[0],"COORDINATE_ZONE") )
246 : {
247 1 : nZone = atoi(papszTokens[1]);
248 : }
249 34 : else if( EQUAL(papszTokens[0],"HORIZONTAL_UNITS") )
250 : {
251 1 : if( EQUAL(papszTokens[1],"METERS") )
252 1 : pszUnits = "UNIT[\"metre\",1]";
253 0 : else if( EQUAL(papszTokens[1],"FEET") )
254 0 : pszUnits = "UNIT[\"US survey foot\",0.304800609601219]";
255 : }
256 33 : else if( EQUAL(papszTokens[0],"HORIZONTAL_DATUM") )
257 : {
258 1 : if( EQUAL(papszTokens[1],"NAD27") )
259 : {
260 0 : pszDatumLong = NAD27_DATUM;
261 0 : pszDatumShort = "NAD 27";
262 : }
263 1 : else if( EQUAL(papszTokens[1],"WGS72") )
264 : {
265 0 : pszDatumLong = WGS72_DATUM;
266 0 : pszDatumShort = "WGS 72";
267 : }
268 1 : else if( EQUAL(papszTokens[1],"WGS84") )
269 : {
270 0 : pszDatumLong = WGS84_DATUM;
271 0 : pszDatumShort = "WGS 84";
272 : }
273 1 : else if( EQUAL(papszTokens[1],"NAD83") )
274 : {
275 1 : pszDatumLong = NAD83_DATUM;
276 1 : pszDatumShort = "NAD 83";
277 : }
278 : else
279 : {
280 0 : pszDatumLong = "DATUM[\"unknown\"]";
281 0 : pszDatumShort = "unknown";
282 : }
283 : }
284 : else
285 : {
286 : /* we want to generically capture all the other metadata */
287 32 : CPLString osMetaDataValue;
288 : int iToken;
289 :
290 250 : for( iToken = 1; papszTokens[iToken] != NULL; iToken++ )
291 : {
292 218 : if( EQUAL(papszTokens[iToken],"*") )
293 1 : continue;
294 :
295 217 : if( iToken > 1 )
296 186 : osMetaDataValue += " " ;
297 217 : osMetaDataValue += papszTokens[iToken];
298 : }
299 : papszMetadata = CSLAddNameValue( papszMetadata,
300 32 : papszTokens[0],
301 64 : osMetaDataValue );
302 : }
303 :
304 45 : CSLDestroy( papszTokens );
305 : }
306 :
307 1 : CPLReadLineL( NULL );
308 :
309 : /* -------------------------------------------------------------------- */
310 : /* Do these values look coherent for a DOQ file? It would be */
311 : /* nice to do a more comprehensive test than this! */
312 : /* -------------------------------------------------------------------- */
313 1 : if( nWidth < 500 || nWidth > 25000
314 : || nHeight < 500 || nHeight > 25000
315 : || nBandStorage < 0 || nBandStorage > 4
316 : || nBandTypes < 1 || nBandTypes > 9 )
317 : {
318 0 : CSLDestroy( papszMetadata );
319 0 : VSIFCloseL(fp);
320 0 : return NULL;
321 : }
322 :
323 : /* -------------------------------------------------------------------- */
324 : /* Check the configuration. We don't currently handle all */
325 : /* variations, only the common ones. */
326 : /* -------------------------------------------------------------------- */
327 1 : if( nBandTypes > 5 )
328 : {
329 : CPLError( CE_Failure, CPLE_OpenFailed,
330 : "DOQ Data Type (%d) is not a supported configuration.\n",
331 0 : nBandTypes );
332 0 : CSLDestroy( papszMetadata );
333 0 : VSIFCloseL(fp);
334 0 : return NULL;
335 : }
336 :
337 : /* -------------------------------------------------------------------- */
338 : /* Confirm the requested access is supported. */
339 : /* -------------------------------------------------------------------- */
340 1 : if( poOpenInfo->eAccess == GA_Update )
341 : {
342 0 : CSLDestroy( papszMetadata );
343 : CPLError( CE_Failure, CPLE_NotSupported,
344 : "The DOQ2 driver does not support update access to existing"
345 0 : " datasets.\n" );
346 0 : VSIFCloseL(fp);
347 0 : return NULL;
348 : }
349 : /* -------------------------------------------------------------------- */
350 : /* Create a corresponding GDALDataset. */
351 : /* -------------------------------------------------------------------- */
352 : DOQ2Dataset *poDS;
353 :
354 1 : poDS = new DOQ2Dataset();
355 :
356 1 : poDS->nRasterXSize = nWidth;
357 1 : poDS->nRasterYSize = nHeight;
358 :
359 1 : poDS->SetMetadata( papszMetadata );
360 1 : CSLDestroy( papszMetadata );
361 :
362 1 : poDS->fpImage = fp;
363 :
364 : /* -------------------------------------------------------------------- */
365 : /* Compute layout of data. */
366 : /* -------------------------------------------------------------------- */
367 1 : if( nBandCount < 2 )
368 0 : nBandCount = nBytesPerPixel;
369 : else
370 1 : nBytesPerPixel *= nBandCount;
371 :
372 1 : nBytesPerLine = nBytesPerPixel * nWidth;
373 :
374 : /* -------------------------------------------------------------------- */
375 : /* Create band information objects. */
376 : /* -------------------------------------------------------------------- */
377 8 : for( i = 0; i < nBandCount; i++ )
378 : {
379 : poDS->SetBand( i+1,
380 : new RawRasterBand( poDS, i+1, poDS->fpImage,
381 : nSkipBytes + i, nBytesPerPixel, nBytesPerLine,
382 3 : GDT_Byte, TRUE, TRUE ) );
383 : }
384 :
385 1 : if (nProjType == 1)
386 : {
387 : poDS->pszProjection =
388 : CPLStrdup(CPLSPrintf( UTM_FORMAT, pszDatumShort, nZone,
389 1 : pszDatumLong, nZone * 6 - 183, pszUnits ));
390 : }
391 : else
392 : {
393 0 : poDS->pszProjection = CPLStrdup("");
394 : }
395 :
396 1 : poDS->dfULX = dfULXMap;
397 1 : poDS->dfULY = dfULYMap;
398 :
399 1 : poDS->dfXPixelSize = dfXDim;
400 1 : poDS->dfYPixelSize = dfYDim;
401 :
402 1 : if ( pszQuadname ) CPLFree( pszQuadname );
403 1 : if ( pszQuadquad) CPLFree( pszQuadquad );
404 1 : if ( pszState) CPLFree( pszState );
405 :
406 : /* -------------------------------------------------------------------- */
407 : /* Initialize any PAM information. */
408 : /* -------------------------------------------------------------------- */
409 1 : poDS->SetDescription( poOpenInfo->pszFilename );
410 1 : poDS->TryLoadXML();
411 :
412 : /* -------------------------------------------------------------------- */
413 : /* Check for overviews. */
414 : /* -------------------------------------------------------------------- */
415 1 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
416 :
417 1 : return( poDS );
418 : }
419 :
420 : /************************************************************************/
421 : /* GDALRegister_DOQ1() */
422 : /************************************************************************/
423 :
424 582 : void GDALRegister_DOQ2()
425 :
426 : {
427 : GDALDriver *poDriver;
428 :
429 582 : if( GDALGetDriverByName( "DOQ2" ) == NULL )
430 : {
431 561 : poDriver = new GDALDriver();
432 :
433 561 : poDriver->SetDescription( "DOQ2" );
434 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
435 561 : "USGS DOQ (New Style)" );
436 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
437 561 : "frmt_various.html#DOQ2" );
438 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
439 :
440 561 : poDriver->pfnOpen = DOQ2Dataset::Open;
441 :
442 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
443 : }
444 582 : }
445 :
|
