Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: USGS DOQ Driver (Second Generation Format)
4 : * Purpose: Implementation of DOQ2Dataset
5 : * Author: Derrick J Brashear, shadow@dementia.org
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2000, Derrick J Brashear
9 : * Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_string.h"
15 : #include "gdal_frmts.h"
16 : #include "gdal_priv.h"
17 : #include "rawdataset.h"
18 :
19 : #ifndef UTM_FORMAT_defined
20 : #define UTM_FORMAT_defined
21 :
22 : static const char UTM_FORMAT[] =
23 : "PROJCS[\"%s / UTM zone %dN\",GEOGCS[%s,PRIMEM[\"Greenwich\",0],"
24 : "UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],"
25 : "PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%d],"
26 : "PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],"
27 : "PARAMETER[\"false_northing\",0],%s]";
28 :
29 : static const char WGS84_DATUM[] =
30 : "\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]]";
31 :
32 : static const char WGS72_DATUM[] =
33 : "\"WGS 72\",DATUM[\"WGS_1972\",SPHEROID[\"NWL 10D\",6378135,298.26]]";
34 :
35 : static const char NAD27_DATUM[] =
36 : "\"NAD27\",DATUM[\"North_American_Datum_1927\","
37 : "SPHEROID[\"Clarke 1866\",6378206.4,294.978698213901]]";
38 :
39 : static const char NAD83_DATUM[] =
40 : "\"NAD83\",DATUM[\"North_American_Datum_1983\","
41 : "SPHEROID[\"GRS 1980\",6378137,298.257222101]]";
42 :
43 : #endif
44 :
45 : /************************************************************************/
46 : /* ==================================================================== */
47 : /* DOQ2Dataset */
48 : /* ==================================================================== */
49 : /************************************************************************/
50 :
51 : class DOQ2Dataset final : public RawDataset
52 : {
53 : VSILFILE *fpImage = nullptr; // Image data file.
54 :
55 : double dfULX = 0;
56 : double dfULY = 0;
57 : double dfXPixelSize = 0;
58 : double dfYPixelSize = 0;
59 :
60 : OGRSpatialReference m_oSRS{};
61 :
62 : CPL_DISALLOW_COPY_ASSIGN(DOQ2Dataset)
63 :
64 : CPLErr Close() override;
65 :
66 : public:
67 : DOQ2Dataset();
68 : ~DOQ2Dataset() override;
69 :
70 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
71 :
72 0 : const OGRSpatialReference *GetSpatialRef() const override
73 : {
74 0 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
75 : }
76 :
77 : static GDALDataset *Open(GDALOpenInfo *);
78 : };
79 :
80 : /************************************************************************/
81 : /* DOQ2Dataset() */
82 : /************************************************************************/
83 :
84 1 : DOQ2Dataset::DOQ2Dataset()
85 : {
86 1 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
87 1 : }
88 :
89 : /************************************************************************/
90 : /* ~DOQ2Dataset() */
91 : /************************************************************************/
92 :
93 2 : DOQ2Dataset::~DOQ2Dataset()
94 :
95 : {
96 1 : DOQ2Dataset::Close();
97 2 : }
98 :
99 : /************************************************************************/
100 : /* Close() */
101 : /************************************************************************/
102 :
103 2 : CPLErr DOQ2Dataset::Close()
104 : {
105 2 : CPLErr eErr = CE_None;
106 2 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
107 : {
108 1 : if (DOQ2Dataset::FlushCache(true) != CE_None)
109 0 : eErr = CE_Failure;
110 :
111 1 : if (fpImage)
112 : {
113 1 : if (VSIFCloseL(fpImage) != 0)
114 : {
115 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
116 0 : eErr = CE_Failure;
117 : }
118 : }
119 :
120 1 : if (GDALPamDataset::Close() != CE_None)
121 0 : eErr = CE_Failure;
122 : }
123 2 : return eErr;
124 : }
125 :
126 : /************************************************************************/
127 : /* GetGeoTransform() */
128 : /************************************************************************/
129 :
130 1 : CPLErr DOQ2Dataset::GetGeoTransform(GDALGeoTransform >) const
131 :
132 : {
133 1 : gt[0] = dfULX;
134 1 : gt[1] = dfXPixelSize;
135 1 : gt[2] = 0.0;
136 1 : gt[3] = dfULY;
137 1 : gt[4] = 0.0;
138 1 : gt[5] = -1 * dfYPixelSize;
139 :
140 1 : return CE_None;
141 : }
142 :
143 : /************************************************************************/
144 : /* Open() */
145 : /************************************************************************/
146 :
147 34388 : GDALDataset *DOQ2Dataset::Open(GDALOpenInfo *poOpenInfo)
148 :
149 : {
150 : /* -------------------------------------------------------------------- */
151 : /* We assume the user is pointing to the binary (i.e. .bil) file. */
152 : /* -------------------------------------------------------------------- */
153 34388 : if (poOpenInfo->nHeaderBytes < 212 || poOpenInfo->fpL == nullptr)
154 31329 : return nullptr;
155 :
156 3059 : if (!STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
157 : "BEGIN_USGS_DOQ_HEADER"))
158 3058 : return nullptr;
159 :
160 : /* -------------------------------------------------------------------- */
161 : /* Confirm the requested access is supported. */
162 : /* -------------------------------------------------------------------- */
163 1 : if (poOpenInfo->eAccess == GA_Update)
164 : {
165 0 : ReportUpdateNotSupportedByDriver("DOQ2");
166 0 : return nullptr;
167 : }
168 :
169 1 : int nBytesPerPixel = 0;
170 1 : int nWidth = 0;
171 1 : int nHeight = 0;
172 1 : int nBandStorage = 0;
173 1 : int nBandTypes = 0;
174 1 : const char *pszDatumLong = nullptr;
175 1 : const char *pszDatumShort = nullptr;
176 1 : const char *pszUnits = nullptr;
177 1 : int nZone = 0;
178 1 : int nProjType = 0;
179 1 : int nSkipBytes = 0;
180 1 : int nBandCount = 0;
181 1 : double dfULXMap = 0.0;
182 1 : double dfULYMap = 0.0;
183 1 : double dfXDim = 0.0;
184 1 : double dfYDim = 0.0;
185 1 : char **papszMetadata = nullptr;
186 :
187 : /* read and discard the first line */
188 1 : CPL_IGNORE_RET_VAL(CPLReadLineL(poOpenInfo->fpL));
189 :
190 1 : const char *pszLine = nullptr;
191 46 : while ((pszLine = CPLReadLineL(poOpenInfo->fpL)) != nullptr)
192 : {
193 46 : if (EQUAL(pszLine, "END_USGS_DOQ_HEADER"))
194 0 : break;
195 :
196 46 : char **papszTokens = CSLTokenizeString(pszLine);
197 46 : if (CSLCount(papszTokens) < 2)
198 : {
199 1 : CSLDestroy(papszTokens);
200 1 : break;
201 : }
202 :
203 46 : if (EQUAL(papszTokens[0], "SAMPLES_AND_LINES") &&
204 1 : CSLCount(papszTokens) >= 3)
205 : {
206 1 : nWidth = atoi(papszTokens[1]);
207 1 : nHeight = atoi(papszTokens[2]);
208 : }
209 44 : else if (EQUAL(papszTokens[0], "BYTE_COUNT"))
210 : {
211 1 : nSkipBytes = atoi(papszTokens[1]);
212 : }
213 44 : else if (EQUAL(papszTokens[0], "XY_ORIGIN") &&
214 1 : CSLCount(papszTokens) >= 3)
215 : {
216 1 : dfULXMap = CPLAtof(papszTokens[1]);
217 1 : dfULYMap = CPLAtof(papszTokens[2]);
218 : }
219 42 : else if (EQUAL(papszTokens[0], "HORIZONTAL_RESOLUTION"))
220 : {
221 1 : dfXDim = CPLAtof(papszTokens[1]);
222 1 : dfYDim = dfXDim;
223 : }
224 41 : else if (EQUAL(papszTokens[0], "BAND_ORGANIZATION"))
225 : {
226 1 : if (EQUAL(papszTokens[1], "SINGLE FILE"))
227 0 : nBandStorage = 1;
228 1 : if (EQUAL(papszTokens[1], "BSQ"))
229 0 : nBandStorage = 1;
230 1 : if (EQUAL(papszTokens[1], "BIL"))
231 0 : nBandStorage = 1;
232 1 : if (EQUAL(papszTokens[1], "BIP"))
233 1 : nBandStorage = 4;
234 : }
235 40 : else if (EQUAL(papszTokens[0], "BAND_CONTENT"))
236 : {
237 3 : if (EQUAL(papszTokens[1], "BLACK&WHITE"))
238 0 : nBandTypes = 1;
239 3 : else if (EQUAL(papszTokens[1], "COLOR"))
240 0 : nBandTypes = 5;
241 3 : else if (EQUAL(papszTokens[1], "RGB"))
242 0 : nBandTypes = 5;
243 3 : else if (EQUAL(papszTokens[1], "RED"))
244 1 : nBandTypes = 5;
245 2 : else if (EQUAL(papszTokens[1], "GREEN"))
246 1 : nBandTypes = 5;
247 1 : else if (EQUAL(papszTokens[1], "BLUE"))
248 1 : nBandTypes = 5;
249 :
250 3 : nBandCount++;
251 : }
252 37 : else if (EQUAL(papszTokens[0], "BITS_PER_PIXEL"))
253 : {
254 1 : nBytesPerPixel = atoi(papszTokens[1]) / 8;
255 : }
256 36 : else if (EQUAL(papszTokens[0], "HORIZONTAL_COORDINATE_SYSTEM"))
257 : {
258 1 : if (EQUAL(papszTokens[1], "UTM"))
259 1 : nProjType = 1;
260 0 : else if (EQUAL(papszTokens[1], "SPCS"))
261 0 : nProjType = 2;
262 0 : else if (EQUAL(papszTokens[1], "GEOGRAPHIC"))
263 0 : nProjType = 0;
264 : }
265 35 : else if (EQUAL(papszTokens[0], "COORDINATE_ZONE"))
266 : {
267 1 : nZone = atoi(papszTokens[1]);
268 : }
269 34 : else if (EQUAL(papszTokens[0], "HORIZONTAL_UNITS"))
270 : {
271 1 : if (EQUAL(papszTokens[1], "METERS"))
272 1 : pszUnits = "UNIT[\"metre\",1]";
273 0 : else if (EQUAL(papszTokens[1], "FEET"))
274 0 : pszUnits = "UNIT[\"US survey foot\",0.304800609601219]";
275 : }
276 33 : else if (EQUAL(papszTokens[0], "HORIZONTAL_DATUM"))
277 : {
278 1 : if (EQUAL(papszTokens[1], "NAD27"))
279 : {
280 0 : pszDatumLong = NAD27_DATUM;
281 0 : pszDatumShort = "NAD 27";
282 : }
283 1 : else if (EQUAL(papszTokens[1], " WGS72"))
284 : {
285 0 : pszDatumLong = WGS72_DATUM;
286 0 : pszDatumShort = "WGS 72";
287 : }
288 1 : else if (EQUAL(papszTokens[1], "WGS84"))
289 : {
290 0 : pszDatumLong = WGS84_DATUM;
291 0 : pszDatumShort = "WGS 84";
292 : }
293 1 : else if (EQUAL(papszTokens[1], "NAD83"))
294 : {
295 1 : pszDatumLong = NAD83_DATUM;
296 1 : pszDatumShort = "NAD 83";
297 : }
298 : else
299 : {
300 0 : pszDatumLong = "DATUM[\"unknown\"]";
301 0 : pszDatumShort = "unknown";
302 : }
303 : }
304 : else
305 : {
306 : /* we want to generically capture all the other metadata */
307 32 : CPLString osMetaDataValue;
308 :
309 250 : for (int iToken = 1; papszTokens[iToken] != nullptr; iToken++)
310 : {
311 218 : if (EQUAL(papszTokens[iToken], "*"))
312 1 : continue;
313 :
314 217 : if (iToken > 1)
315 186 : osMetaDataValue += " ";
316 217 : osMetaDataValue += papszTokens[iToken];
317 : }
318 : papszMetadata =
319 32 : CSLAddNameValue(papszMetadata, papszTokens[0], osMetaDataValue);
320 : }
321 :
322 45 : CSLDestroy(papszTokens);
323 : }
324 :
325 1 : CPLReadLineL(nullptr);
326 :
327 : /* -------------------------------------------------------------------- */
328 : /* Do these values look coherent for a DOQ file? It would be */
329 : /* nice to do a more comprehensive test than this! */
330 : /* -------------------------------------------------------------------- */
331 1 : if (nWidth < 500 || nWidth > 25000 || nHeight < 500 || nHeight > 25000 ||
332 1 : nBandStorage < 0 || nBandStorage > 4 || nBandTypes < 1 ||
333 1 : nBandTypes > 9 || nBytesPerPixel < 0)
334 : {
335 0 : CSLDestroy(papszMetadata);
336 0 : return nullptr;
337 : }
338 :
339 : /* -------------------------------------------------------------------- */
340 : /* Check the configuration. We don't currently handle all */
341 : /* variations, only the common ones. */
342 : /* -------------------------------------------------------------------- */
343 1 : if (nBandTypes > 5)
344 : {
345 0 : CPLError(CE_Failure, CPLE_OpenFailed,
346 : "DOQ Data Type (%d) is not a supported configuration.",
347 : nBandTypes);
348 0 : CSLDestroy(papszMetadata);
349 0 : return nullptr;
350 : }
351 :
352 : /* -------------------------------------------------------------------- */
353 : /* Confirm the requested access is supported. */
354 : /* -------------------------------------------------------------------- */
355 1 : if (poOpenInfo->eAccess == GA_Update)
356 : {
357 0 : CSLDestroy(papszMetadata);
358 0 : ReportUpdateNotSupportedByDriver("DOQ2");
359 0 : return nullptr;
360 : }
361 : /* -------------------------------------------------------------------- */
362 : /* Create a corresponding GDALDataset. */
363 : /* -------------------------------------------------------------------- */
364 2 : auto poDS = std::make_unique<DOQ2Dataset>();
365 :
366 1 : poDS->nRasterXSize = nWidth;
367 1 : poDS->nRasterYSize = nHeight;
368 :
369 1 : poDS->SetMetadata(papszMetadata);
370 1 : CSLDestroy(papszMetadata);
371 1 : papszMetadata = nullptr;
372 :
373 1 : poDS->fpImage = poOpenInfo->fpL;
374 1 : poOpenInfo->fpL = nullptr;
375 :
376 : /* -------------------------------------------------------------------- */
377 : /* Compute layout of data. */
378 : /* -------------------------------------------------------------------- */
379 1 : if (nBandCount < 2)
380 : {
381 0 : nBandCount = nBytesPerPixel;
382 0 : if (!GDALCheckBandCount(nBandCount, FALSE))
383 : {
384 0 : return nullptr;
385 : }
386 : }
387 : else
388 : {
389 1 : if (nBytesPerPixel > INT_MAX / nBandCount)
390 : {
391 0 : return nullptr;
392 : }
393 1 : nBytesPerPixel *= nBandCount;
394 : }
395 :
396 1 : if (nBytesPerPixel > INT_MAX / nWidth)
397 : {
398 0 : return nullptr;
399 : }
400 1 : const int nBytesPerLine = nBytesPerPixel * nWidth;
401 :
402 : /* -------------------------------------------------------------------- */
403 : /* Create band information objects. */
404 : /* -------------------------------------------------------------------- */
405 4 : for (int i = 0; i < nBandCount; i++)
406 : {
407 : auto poBand = RawRasterBand::Create(
408 6 : poDS.get(), i + 1, poDS->fpImage, nSkipBytes + i, nBytesPerPixel,
409 : nBytesPerLine, GDT_Byte,
410 : RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
411 3 : RawRasterBand::OwnFP::NO);
412 3 : if (!poBand)
413 0 : return nullptr;
414 3 : poDS->SetBand(i + 1, std::move(poBand));
415 : }
416 :
417 1 : if (nProjType == 1)
418 : {
419 2 : poDS->m_oSRS.importFromWkt(
420 : CPLSPrintf(UTM_FORMAT, pszDatumShort ? pszDatumShort : "", nZone,
421 : pszDatumLong ? pszDatumLong : "",
422 1 : nZone >= 1 && nZone <= 60 ? nZone * 6 - 183 : 0,
423 : pszUnits ? pszUnits : ""));
424 : }
425 :
426 1 : poDS->dfULX = dfULXMap;
427 1 : poDS->dfULY = dfULYMap;
428 :
429 1 : poDS->dfXPixelSize = dfXDim;
430 1 : poDS->dfYPixelSize = dfYDim;
431 :
432 : /* -------------------------------------------------------------------- */
433 : /* Initialize any PAM information. */
434 : /* -------------------------------------------------------------------- */
435 1 : poDS->SetDescription(poOpenInfo->pszFilename);
436 1 : poDS->TryLoadXML();
437 :
438 : /* -------------------------------------------------------------------- */
439 : /* Check for overviews. */
440 : /* -------------------------------------------------------------------- */
441 1 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
442 :
443 1 : return poDS.release();
444 : }
445 :
446 : /************************************************************************/
447 : /* GDALRegister_DOQ1() */
448 : /************************************************************************/
449 :
450 2033 : void GDALRegister_DOQ2()
451 :
452 : {
453 2033 : if (GDALGetDriverByName("DOQ2") != nullptr)
454 283 : return;
455 :
456 1750 : GDALDriver *poDriver = new GDALDriver();
457 :
458 1750 : poDriver->SetDescription("DOQ2");
459 1750 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
460 1750 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "USGS DOQ (New Style)");
461 1750 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/doq2.html");
462 1750 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
463 :
464 1750 : poDriver->pfnOpen = DOQ2Dataset::Open;
465 :
466 1750 : GetGDALDriverManager()->RegisterDriver(poDriver);
467 : }
|
