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