1 : /******************************************************************************
2 : * $Id: doq2dataset.cpp 17664 2009-09-21 21:16:45Z 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 17664 2009-09-21 21:16:45Z 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 : FILE *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 2 : DOQ2Dataset::~DOQ2Dataset()
94 :
95 : {
96 1 : FlushCache();
97 :
98 1 : CPLFree( pszProjection );
99 1 : if( fpImage != NULL )
100 1 : VSIFClose( fpImage );
101 2 : }
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 8735 : GDALDataset *DOQ2Dataset::Open( GDALOpenInfo * poOpenInfo )
135 :
136 : {
137 8735 : 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 8735 : if( poOpenInfo->nHeaderBytes < 212 || poOpenInfo->fp == NULL )
143 8485 : return NULL;
144 :
145 250 : int nLineCount = 0;
146 : const char *pszLine;
147 250 : int nBytesPerPixel=0;
148 250 : const char *pszDatumLong=NULL, *pszDatumShort=NULL;
149 250 : const char *pszUnits=NULL;
150 250 : char *pszQuadname = NULL;
151 250 : char *pszQuadquad = NULL;
152 250 : char *pszState = NULL;
153 250 : int nZone=0, nProjType=0;
154 250 : int nSkipBytes=0, nBytesPerLine, i, nBandCount = 0;
155 250 : double dfULXMap=0.0, dfULYMap = 0.0;
156 250 : double dfXDim=0.0, dfYDim=0.0;
157 250 : char **papszMetadata = NULL;
158 :
159 250 : if(! EQUALN((const char *) poOpenInfo->pabyHeader,
160 : "BEGIN_USGS_DOQ_HEADER", 21) )
161 249 : return NULL;
162 :
163 : /* read and discard the first line */
164 1 : pszLine = CPLReadLine( poOpenInfo->fp );
165 :
166 47 : while( (pszLine = CPLReadLine( poOpenInfo->fp )) )
167 : {
168 : char **papszTokens;
169 :
170 46 : nLineCount++;
171 :
172 46 : if( EQUAL(pszLine,"END_USGS_DOQ_HEADER") )
173 0 : break;
174 :
175 46 : papszTokens = CSLTokenizeString( pszLine );
176 46 : if( CSLCount( papszTokens ) < 2 )
177 : {
178 1 : CSLDestroy( papszTokens );
179 1 : break;
180 : }
181 :
182 45 : if( EQUAL(papszTokens[0],"SAMPLES_AND_LINES") && CSLCount(papszTokens) >= 3 )
183 : {
184 1 : nWidth = atoi(papszTokens[1]);
185 1 : nHeight = atoi(papszTokens[2]);
186 : }
187 44 : else if( EQUAL(papszTokens[0],"BYTE_COUNT") )
188 : {
189 1 : nSkipBytes = atoi(papszTokens[1]);
190 : }
191 43 : else if( EQUAL(papszTokens[0],"XY_ORIGIN") && CSLCount(papszTokens) >= 3 )
192 : {
193 1 : dfULXMap = atof(papszTokens[1]);
194 1 : dfULYMap = atof(papszTokens[2]);
195 : }
196 42 : else if( EQUAL(papszTokens[0],"HORIZONTAL_RESOLUTION") )
197 : {
198 1 : dfXDim = dfYDim = atof(papszTokens[1]);
199 : }
200 41 : else if( EQUAL(papszTokens[0],"BAND_ORGANIZATION") )
201 : {
202 1 : if( EQUAL(papszTokens[1],"SINGLE FILE") )
203 0 : nBandStorage = 1;
204 1 : if( EQUAL(papszTokens[1],"BSQ") )
205 0 : nBandStorage = 1;
206 1 : if( EQUAL(papszTokens[1],"BIL") )
207 0 : nBandStorage = 1;
208 1 : if( EQUAL(papszTokens[1],"BIP") )
209 1 : nBandStorage = 4;
210 : }
211 40 : else if( EQUAL(papszTokens[0],"BAND_CONTENT") )
212 : {
213 3 : if( EQUAL(papszTokens[1],"BLACK&WHITE") )
214 0 : nBandTypes = 1;
215 3 : else if( EQUAL(papszTokens[1],"COLOR") )
216 0 : nBandTypes = 5;
217 3 : else if( EQUAL(papszTokens[1],"RGB") )
218 0 : nBandTypes = 5;
219 3 : else if( EQUAL(papszTokens[1],"RED") )
220 1 : nBandTypes = 5;
221 2 : else if( EQUAL(papszTokens[1],"GREEN") )
222 1 : nBandTypes = 5;
223 1 : else if( EQUAL(papszTokens[1],"BLUE") )
224 1 : nBandTypes = 5;
225 :
226 3 : nBandCount++;
227 : }
228 37 : else if( EQUAL(papszTokens[0],"BITS_PER_PIXEL") )
229 : {
230 1 : nBytesPerPixel = (atoi(papszTokens[1]) / 8);
231 : }
232 36 : else if( EQUAL(papszTokens[0],"HORIZONTAL_COORDINATE_SYSTEM") )
233 : {
234 1 : if( EQUAL(papszTokens[1],"UTM") )
235 1 : nProjType = 1;
236 0 : else if( EQUAL(papszTokens[1],"SPCS") )
237 0 : nProjType = 2;
238 0 : else if( EQUAL(papszTokens[1],"GEOGRAPHIC") )
239 0 : nProjType = 0;
240 : }
241 35 : else if( EQUAL(papszTokens[0],"COORDINATE_ZONE") )
242 : {
243 1 : nZone = atoi(papszTokens[1]);
244 : }
245 34 : else if( EQUAL(papszTokens[0],"HORIZONTAL_UNITS") )
246 : {
247 1 : if( EQUAL(papszTokens[1],"METERS") )
248 1 : pszUnits = "UNIT[\"metre\",1]";
249 0 : else if( EQUAL(papszTokens[1],"FEET") )
250 0 : pszUnits = "UNIT[\"US survey foot\",0.304800609601219]";
251 : }
252 33 : else if( EQUAL(papszTokens[0],"HORIZONTAL_DATUM") )
253 : {
254 1 : if( EQUAL(papszTokens[1],"NAD27") )
255 : {
256 0 : pszDatumLong = NAD27_DATUM;
257 0 : pszDatumShort = "NAD 27";
258 : }
259 1 : else if( EQUAL(papszTokens[1],"WGS72") )
260 : {
261 0 : pszDatumLong = WGS72_DATUM;
262 0 : pszDatumShort = "WGS 72";
263 : }
264 1 : else if( EQUAL(papszTokens[1],"WGS84") )
265 : {
266 0 : pszDatumLong = WGS84_DATUM;
267 0 : pszDatumShort = "WGS 84";
268 : }
269 1 : else if( EQUAL(papszTokens[1],"NAD83") )
270 : {
271 1 : pszDatumLong = NAD83_DATUM;
272 1 : pszDatumShort = "NAD 83";
273 : }
274 : else
275 : {
276 0 : pszDatumLong = "DATUM[\"unknown\"]";
277 0 : pszDatumShort = "unknown";
278 : }
279 : }
280 : else
281 : {
282 : /* we want to generically capture all the other metadata */
283 32 : CPLString osMetaDataValue;
284 : int iToken;
285 :
286 250 : for( iToken = 1; papszTokens[iToken] != NULL; iToken++ )
287 : {
288 218 : if( EQUAL(papszTokens[iToken],"*") )
289 1 : continue;
290 :
291 217 : if( iToken > 1 )
292 186 : osMetaDataValue += " " ;
293 217 : osMetaDataValue += papszTokens[iToken];
294 : }
295 : papszMetadata = CSLAddNameValue( papszMetadata,
296 32 : papszTokens[0],
297 64 : osMetaDataValue );
298 : }
299 :
300 45 : CSLDestroy( papszTokens );
301 : }
302 :
303 1 : CPLReadLine( NULL );
304 :
305 : /* -------------------------------------------------------------------- */
306 : /* Do these values look coherent for a DOQ file? It would be */
307 : /* nice to do a more comprehensive test than this! */
308 : /* -------------------------------------------------------------------- */
309 1 : if( nWidth < 500 || nWidth > 25000
310 : || nHeight < 500 || nHeight > 25000
311 : || nBandStorage < 0 || nBandStorage > 4
312 : || nBandTypes < 1 || nBandTypes > 9 )
313 : {
314 0 : CSLDestroy( papszMetadata );
315 0 : return NULL;
316 : }
317 :
318 : /* -------------------------------------------------------------------- */
319 : /* Check the configuration. We don't currently handle all */
320 : /* variations, only the common ones. */
321 : /* -------------------------------------------------------------------- */
322 1 : if( nBandTypes > 5 )
323 : {
324 : CPLError( CE_Failure, CPLE_OpenFailed,
325 : "DOQ Data Type (%d) is not a supported configuration.\n",
326 0 : nBandTypes );
327 0 : CSLDestroy( papszMetadata );
328 0 : return NULL;
329 : }
330 :
331 : /* -------------------------------------------------------------------- */
332 : /* Confirm the requested access is supported. */
333 : /* -------------------------------------------------------------------- */
334 1 : if( poOpenInfo->eAccess == GA_Update )
335 : {
336 0 : CSLDestroy( papszMetadata );
337 : CPLError( CE_Failure, CPLE_NotSupported,
338 : "The DOQ2 driver does not support update access to existing"
339 0 : " datasets.\n" );
340 0 : return NULL;
341 : }
342 : /* -------------------------------------------------------------------- */
343 : /* Create a corresponding GDALDataset. */
344 : /* -------------------------------------------------------------------- */
345 : DOQ2Dataset *poDS;
346 :
347 1 : poDS = new DOQ2Dataset();
348 :
349 1 : poDS->nRasterXSize = nWidth;
350 1 : poDS->nRasterYSize = nHeight;
351 :
352 1 : poDS->SetMetadata( papszMetadata );
353 1 : CSLDestroy( papszMetadata );
354 :
355 : /* -------------------------------------------------------------------- */
356 : /* Assume ownership of the file handled from the GDALOpenInfo. */
357 : /* -------------------------------------------------------------------- */
358 1 : poDS->fpImage = poOpenInfo->fp;
359 1 : poOpenInfo->fp = NULL;
360 :
361 : /* -------------------------------------------------------------------- */
362 : /* Compute layout of data. */
363 : /* -------------------------------------------------------------------- */
364 1 : if( nBandCount < 2 )
365 0 : nBandCount = nBytesPerPixel;
366 : else
367 1 : nBytesPerPixel *= nBandCount;
368 :
369 1 : nBytesPerLine = nBytesPerPixel * nWidth;
370 :
371 : /* -------------------------------------------------------------------- */
372 : /* Create band information objects. */
373 : /* -------------------------------------------------------------------- */
374 8 : for( i = 0; i < nBandCount; i++ )
375 : {
376 : poDS->SetBand( i+1,
377 : new RawRasterBand( poDS, i+1, poDS->fpImage,
378 : nSkipBytes + i, nBytesPerPixel, nBytesPerLine,
379 3 : GDT_Byte, TRUE ) );
380 : }
381 :
382 1 : if (nProjType == 1)
383 : {
384 : poDS->pszProjection =
385 : CPLStrdup(CPLSPrintf( UTM_FORMAT, pszDatumShort, nZone,
386 1 : pszDatumLong, nZone * 6 - 183, pszUnits ));
387 : }
388 : else
389 : {
390 0 : poDS->pszProjection = CPLStrdup("");
391 : }
392 :
393 1 : poDS->dfULX = dfULXMap;
394 1 : poDS->dfULY = dfULYMap;
395 :
396 1 : poDS->dfXPixelSize = dfXDim;
397 1 : poDS->dfYPixelSize = dfYDim;
398 :
399 1 : if ( pszQuadname ) CPLFree( pszQuadname );
400 1 : if ( pszQuadquad) CPLFree( pszQuadquad );
401 1 : if ( pszState) CPLFree( pszState );
402 :
403 : /* -------------------------------------------------------------------- */
404 : /* Initialize any PAM information. */
405 : /* -------------------------------------------------------------------- */
406 1 : poDS->SetDescription( poOpenInfo->pszFilename );
407 1 : poDS->TryLoadXML();
408 :
409 : /* -------------------------------------------------------------------- */
410 : /* Check for overviews. */
411 : /* -------------------------------------------------------------------- */
412 1 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
413 :
414 1 : return( poDS );
415 : }
416 :
417 : /************************************************************************/
418 : /* GDALRegister_DOQ1() */
419 : /************************************************************************/
420 :
421 338 : void GDALRegister_DOQ2()
422 :
423 : {
424 : GDALDriver *poDriver;
425 :
426 338 : if( GDALGetDriverByName( "DOQ2" ) == NULL )
427 : {
428 336 : poDriver = new GDALDriver();
429 :
430 336 : poDriver->SetDescription( "DOQ2" );
431 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
432 336 : "USGS DOQ (New Style)" );
433 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
434 336 : "frmt_various.html#DOQ2" );
435 :
436 336 : poDriver->pfnOpen = DOQ2Dataset::Open;
437 :
438 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
439 : }
440 338 : }
441 :
|
