Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: USGS DOQ Driver (First Generation Format)
4 : * Purpose: Implementation of DOQ1Dataset
5 : * Author: Frank Warmerdam, warmerda@home.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 1999, Frank Warmerdam
9 : * Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "gdal_frmts.h"
15 : #include "cpl_string.h"
16 : #include "gdal_priv.h"
17 : #include "rawdataset.h"
18 :
19 : #include <algorithm>
20 : #include <cmath>
21 :
22 : #ifndef UTM_FORMAT_defined
23 : #define UTM_FORMAT_defined
24 :
25 : static const char UTM_FORMAT[] =
26 : "PROJCS[\"%s / UTM zone %dN\",GEOGCS[%s,PRIMEM[\"Greenwich\",0],"
27 : "UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],"
28 : "PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",%d],"
29 : "PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],"
30 : "PARAMETER[\"false_northing\",0],%s]";
31 :
32 : static const char WGS84_DATUM[] =
33 : "\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563]]";
34 :
35 : static const char WGS72_DATUM[] =
36 : "\"WGS 72\",DATUM[\"WGS_1972\",SPHEROID[\"NWL 10D\",6378135,298.26]]";
37 :
38 : static const char NAD27_DATUM[] =
39 : "\"NAD27\",DATUM[\"North_American_Datum_1927\","
40 : "SPHEROID[\"Clarke 1866\",6378206.4,294.978698213901]]";
41 :
42 : static const char NAD83_DATUM[] =
43 : "\"NAD83\",DATUM[\"North_American_Datum_1983\","
44 : "SPHEROID[\"GRS 1980\",6378137,298.257222101]]";
45 :
46 : #endif
47 :
48 : /************************************************************************/
49 : /* DOQGetField() */
50 : /************************************************************************/
51 :
52 12268 : static double DOQGetField(unsigned char *pabyData, int nBytes)
53 :
54 : {
55 12268 : char szWork[128] = {'\0'};
56 :
57 12268 : memcpy(szWork, reinterpret_cast<const char *>(pabyData), nBytes);
58 12268 : szWork[nBytes] = '\0';
59 :
60 67574 : for (int i = 0; i < nBytes; i++)
61 : {
62 55306 : if (szWork[i] == 'D' || szWork[i] == 'd')
63 430 : szWork[i] = 'E';
64 : }
65 :
66 24536 : return CPLAtof(szWork);
67 : }
68 :
69 : /************************************************************************/
70 : /* DOQGetDescription() */
71 : /************************************************************************/
72 :
73 2 : static void DOQGetDescription(GDALDataset *poDS, unsigned char *pabyData)
74 :
75 : {
76 2 : char szWork[128] = {' '};
77 :
78 2 : const char *pszDescBegin = "USGS GeoTIFF DOQ 1:12000 Q-Quad of ";
79 2 : memcpy(szWork, pszDescBegin, strlen(pszDescBegin));
80 2 : memcpy(szWork + strlen(pszDescBegin),
81 : reinterpret_cast<const char *>(pabyData + 0), 38);
82 :
83 2 : int i = 0;
84 2 : while (*(szWork + 72 - i) == ' ')
85 : {
86 0 : i++;
87 : }
88 2 : i--;
89 :
90 2 : memcpy(szWork + 73 - i, reinterpret_cast<const char *>(pabyData + 38), 2);
91 2 : memcpy(szWork + 76 - i, reinterpret_cast<const char *>(pabyData + 44), 2);
92 2 : szWork[77 - i] = '\0';
93 :
94 2 : poDS->SetMetadataItem("DOQ_DESC", szWork);
95 2 : }
96 :
97 : /************************************************************************/
98 : /* ==================================================================== */
99 : /* DOQ1Dataset */
100 : /* ==================================================================== */
101 : /************************************************************************/
102 :
103 : class DOQ1Dataset final : public RawDataset
104 : {
105 : VSILFILE *fpImage = nullptr; // Image data file.
106 :
107 : double dfULX = 0;
108 : double dfULY = 0;
109 : double dfXPixelSize = 0;
110 : double dfYPixelSize = 0;
111 :
112 : OGRSpatialReference m_oSRS{};
113 :
114 : CPL_DISALLOW_COPY_ASSIGN(DOQ1Dataset)
115 :
116 : CPLErr Close() override;
117 :
118 : public:
119 : DOQ1Dataset();
120 : ~DOQ1Dataset() override;
121 :
122 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
123 :
124 0 : const OGRSpatialReference *GetSpatialRef() const override
125 : {
126 0 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
127 : }
128 :
129 : static GDALDataset *Open(GDALOpenInfo *);
130 : };
131 :
132 : /************************************************************************/
133 : /* DOQ1Dataset() */
134 : /************************************************************************/
135 :
136 2 : DOQ1Dataset::DOQ1Dataset()
137 : {
138 2 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
139 2 : }
140 :
141 : /************************************************************************/
142 : /* ~DOQ1Dataset() */
143 : /************************************************************************/
144 :
145 4 : DOQ1Dataset::~DOQ1Dataset()
146 :
147 : {
148 2 : DOQ1Dataset::Close();
149 4 : }
150 :
151 : /************************************************************************/
152 : /* Close() */
153 : /************************************************************************/
154 :
155 4 : CPLErr DOQ1Dataset::Close()
156 : {
157 4 : CPLErr eErr = CE_None;
158 4 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
159 : {
160 2 : if (DOQ1Dataset::FlushCache(true) != CE_None)
161 0 : eErr = CE_Failure;
162 :
163 2 : if (fpImage)
164 : {
165 2 : if (VSIFCloseL(fpImage) != 0)
166 : {
167 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
168 0 : eErr = CE_Failure;
169 : }
170 : }
171 :
172 2 : if (GDALPamDataset::Close() != CE_None)
173 0 : eErr = CE_Failure;
174 : }
175 4 : return eErr;
176 : }
177 :
178 : /************************************************************************/
179 : /* GetGeoTransform() */
180 : /************************************************************************/
181 :
182 0 : CPLErr DOQ1Dataset::GetGeoTransform(GDALGeoTransform >) const
183 :
184 : {
185 0 : gt[0] = dfULX;
186 0 : gt[1] = dfXPixelSize;
187 0 : gt[2] = 0.0;
188 0 : gt[3] = dfULY;
189 0 : gt[4] = 0.0;
190 0 : gt[5] = -1 * dfYPixelSize;
191 :
192 0 : return CE_None;
193 : }
194 :
195 : /************************************************************************/
196 : /* Open() */
197 : /************************************************************************/
198 :
199 34432 : GDALDataset *DOQ1Dataset::Open(GDALOpenInfo *poOpenInfo)
200 :
201 : {
202 : /* -------------------------------------------------------------------- */
203 : /* We assume the user is pointing to the binary (i.e. .bil) file. */
204 : /* -------------------------------------------------------------------- */
205 34432 : if (poOpenInfo->nHeaderBytes < 212 || poOpenInfo->fpL == nullptr)
206 31369 : return nullptr;
207 :
208 : /* -------------------------------------------------------------------- */
209 : /* Attempt to extract a few key values from the header. */
210 : /* -------------------------------------------------------------------- */
211 3063 : const double dfWidth = DOQGetField(poOpenInfo->pabyHeader + 150, 6);
212 3063 : const double dfHeight = DOQGetField(poOpenInfo->pabyHeader + 144, 6);
213 3063 : const double dfBandStorage = DOQGetField(poOpenInfo->pabyHeader + 162, 3);
214 3063 : const double dfBandTypes = DOQGetField(poOpenInfo->pabyHeader + 156, 3);
215 :
216 : /* -------------------------------------------------------------------- */
217 : /* Do these values look coherent for a DOQ file? It would be */
218 : /* nice to do a more comprehensive test than this! */
219 : /* -------------------------------------------------------------------- */
220 57 : if (dfWidth < 500 || dfWidth > 25000 || std::isnan(dfWidth) ||
221 2 : dfHeight < 500 || dfHeight > 25000 || std::isnan(dfHeight) ||
222 2 : dfBandStorage < 0 || dfBandStorage > 4 || std::isnan(dfBandStorage) ||
223 3120 : dfBandTypes < 1 || dfBandTypes > 9 || std::isnan(dfBandTypes))
224 3061 : return nullptr;
225 :
226 2 : const int nWidth = static_cast<int>(dfWidth);
227 2 : const int nHeight = static_cast<int>(dfHeight);
228 : /*const int nBandStorage = static_cast<int>(dfBandStorage);*/
229 2 : const int nBandTypes = static_cast<int>(dfBandTypes);
230 :
231 : /* -------------------------------------------------------------------- */
232 : /* Check the configuration. We don't currently handle all */
233 : /* variations, only the common ones. */
234 : /* -------------------------------------------------------------------- */
235 2 : if (nBandTypes > 5)
236 : {
237 0 : CPLError(CE_Failure, CPLE_OpenFailed,
238 : "DOQ Data Type (%d) is not a supported configuration.",
239 : nBandTypes);
240 0 : return nullptr;
241 : }
242 :
243 : /* -------------------------------------------------------------------- */
244 : /* Confirm the requested access is supported. */
245 : /* -------------------------------------------------------------------- */
246 2 : if (poOpenInfo->eAccess == GA_Update)
247 : {
248 0 : ReportUpdateNotSupportedByDriver("DOQ1");
249 0 : return nullptr;
250 : }
251 :
252 : /* -------------------------------------------------------------------- */
253 : /* Create a corresponding GDALDataset. */
254 : /* -------------------------------------------------------------------- */
255 4 : auto poDS = std::make_unique<DOQ1Dataset>();
256 :
257 : /* -------------------------------------------------------------------- */
258 : /* Capture some information from the file that is of interest. */
259 : /* -------------------------------------------------------------------- */
260 2 : poDS->nRasterXSize = nWidth;
261 2 : poDS->nRasterYSize = nHeight;
262 :
263 2 : std::swap(poDS->fpImage, poOpenInfo->fpL);
264 :
265 : /* -------------------------------------------------------------------- */
266 : /* Compute layout of data. */
267 : /* -------------------------------------------------------------------- */
268 2 : int nBytesPerPixel = 0;
269 :
270 2 : if (nBandTypes < 5)
271 2 : nBytesPerPixel = 1;
272 : else /* if( nBandTypes == 5 ) */
273 0 : nBytesPerPixel = 3;
274 :
275 2 : const int nBytesPerLine = nBytesPerPixel * nWidth;
276 2 : const int nSkipBytes = 4 * nBytesPerLine;
277 :
278 : /* -------------------------------------------------------------------- */
279 : /* Create band information objects. */
280 : /* -------------------------------------------------------------------- */
281 4 : for (int i = 0; i < nBytesPerPixel; i++)
282 : {
283 : auto poBand = RawRasterBand::Create(
284 4 : poDS.get(), i + 1, poDS->fpImage, nSkipBytes + i, nBytesPerPixel,
285 : nBytesPerLine, GDT_Byte,
286 : RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
287 2 : RawRasterBand::OwnFP::NO);
288 2 : if (!poBand)
289 0 : return nullptr;
290 2 : poDS->SetBand(i + 1, std::move(poBand));
291 : }
292 :
293 : /* -------------------------------------------------------------------- */
294 : /* Set the description. */
295 : /* -------------------------------------------------------------------- */
296 2 : DOQGetDescription(poDS.get(), poOpenInfo->pabyHeader);
297 :
298 : /* -------------------------------------------------------------------- */
299 : /* Establish the projection string. */
300 : /* -------------------------------------------------------------------- */
301 2 : if (static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 195, 3)) == 1)
302 : {
303 : int nZone =
304 2 : static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 198, 6));
305 2 : if (nZone < 0 || nZone > 60)
306 0 : nZone = 0;
307 :
308 2 : const char *pszUnits = nullptr;
309 2 : if (static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 204, 3)) == 1)
310 0 : pszUnits = "UNIT[\"US survey foot\",0.304800609601219]";
311 : else
312 2 : pszUnits = "UNIT[\"metre\",1]";
313 :
314 2 : const char *pszDatumLong = nullptr;
315 2 : const char *pszDatumShort = nullptr;
316 2 : switch (static_cast<int>(DOQGetField(poOpenInfo->pabyHeader + 167, 2)))
317 : {
318 0 : case 1:
319 0 : pszDatumLong = NAD27_DATUM;
320 0 : pszDatumShort = "NAD 27";
321 0 : break;
322 :
323 0 : case 2:
324 0 : pszDatumLong = WGS72_DATUM;
325 0 : pszDatumShort = "WGS 72";
326 0 : break;
327 :
328 2 : case 3:
329 2 : pszDatumLong = WGS84_DATUM;
330 2 : pszDatumShort = "WGS 84";
331 2 : break;
332 :
333 0 : case 4:
334 0 : pszDatumLong = NAD83_DATUM;
335 0 : pszDatumShort = "NAD 83";
336 0 : break;
337 :
338 0 : default:
339 0 : pszDatumLong = "DATUM[\"unknown\"]";
340 0 : pszDatumShort = "unknown";
341 0 : break;
342 : }
343 :
344 4 : poDS->m_oSRS.importFromWkt(CPLSPrintf(UTM_FORMAT, pszDatumShort, nZone,
345 2 : pszDatumLong, nZone * 6 - 183,
346 : pszUnits));
347 : }
348 :
349 : /* -------------------------------------------------------------------- */
350 : /* Read the georeferencing information. */
351 : /* -------------------------------------------------------------------- */
352 2 : unsigned char abyRecordData[500] = {'\0'};
353 :
354 4 : if (VSIFSeekL(poDS->fpImage, nBytesPerLine * 2, SEEK_SET) != 0 ||
355 2 : VSIFReadL(abyRecordData, sizeof(abyRecordData), 1, poDS->fpImage) != 1)
356 : {
357 0 : CPLError(CE_Failure, CPLE_FileIO, "Header read error on %s.",
358 : poOpenInfo->pszFilename);
359 0 : return nullptr;
360 : }
361 :
362 2 : poDS->dfULX = DOQGetField(abyRecordData + 288, 24);
363 2 : poDS->dfULY = DOQGetField(abyRecordData + 312, 24);
364 :
365 4 : if (VSIFSeekL(poDS->fpImage, nBytesPerLine * 3, SEEK_SET) != 0 ||
366 2 : VSIFReadL(abyRecordData, sizeof(abyRecordData), 1, poDS->fpImage) != 1)
367 : {
368 0 : CPLError(CE_Failure, CPLE_FileIO, "Header read error on %s.",
369 : poOpenInfo->pszFilename);
370 0 : return nullptr;
371 : }
372 :
373 2 : poDS->dfXPixelSize = DOQGetField(abyRecordData + 59, 12);
374 2 : poDS->dfYPixelSize = DOQGetField(abyRecordData + 71, 12);
375 :
376 : /* -------------------------------------------------------------------- */
377 : /* Initialize any PAM information. */
378 : /* -------------------------------------------------------------------- */
379 2 : poDS->SetDescription(poOpenInfo->pszFilename);
380 2 : poDS->TryLoadXML();
381 :
382 : /* -------------------------------------------------------------------- */
383 : /* Check for overviews. */
384 : /* -------------------------------------------------------------------- */
385 2 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
386 :
387 2 : return poDS.release();
388 : }
389 :
390 : /************************************************************************/
391 : /* GDALRegister_DOQ1() */
392 : /************************************************************************/
393 :
394 2038 : void GDALRegister_DOQ1()
395 :
396 : {
397 2038 : if (GDALGetDriverByName("DOQ1") != nullptr)
398 283 : return;
399 :
400 1755 : GDALDriver *poDriver = new GDALDriver();
401 :
402 1755 : poDriver->SetDescription("DOQ1");
403 1755 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
404 1755 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "USGS DOQ (Old Style)");
405 1755 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/doq1.html");
406 1755 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
407 :
408 1755 : poDriver->pfnOpen = DOQ1Dataset::Open;
409 :
410 1755 : GetGDALDriverManager()->RegisterDriver(poDriver);
411 : }
|
