Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: eCognition
4 : * Purpose: Implementation of Erdas .LAN / .GIS format.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2004, Frank Warmerdam
9 : * Copyright (c) 2008-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 "ogr_spatialref.h"
17 : #include "rawdataset.h"
18 :
19 : #include <cmath>
20 :
21 : #include <algorithm>
22 :
23 : /**
24 :
25 : Erdas Header format: "HEAD74"
26 :
27 : Offset Size Type Description
28 : ------ ---- ---- -----------
29 : 0 6 char magic cookie / version (i.e. HEAD74).
30 : 6 2 Int16 Pixel type, 0=8bit, 1=4bit, 2=16bit
31 : 8 2 Int16 Number of Bands.
32 : 10 6 char Unknown.
33 : 16 4 Int32 Width
34 : 20 4 Int32 Height
35 : 24 4 Int32 X Start (offset in original file?)
36 : 28 4 Int32 Y Start (offset in original file?)
37 : 32 56 char Unknown.
38 : 88 2 Int16 0=LAT, 1=UTM, 2=StatePlane, 3- are projections?
39 : 90 2 Int16 Classes in coverage.
40 : 92 14 char Unknown.
41 : 106 2 Int16 Area Unit (0=none, 1=Acre, 2=Hectare, 3=Other)
42 : 108 4 Float32 Pixel area.
43 : 112 4 Float32 Upper Left corner X (center of pixel?)
44 : 116 4 Float32 Upper Left corner Y (center of pixel?)
45 : 120 4 Float32 Width of a pixel.
46 : 124 4 Float32 Height of a pixel.
47 :
48 : Erdas Header format: "HEADER"
49 :
50 : Offset Size Type Description
51 : ------ ---- ---- -----------
52 : 0 6 char magic cookie / version (i.e. HEAD74).
53 : 6 2 Int16 Pixel type, 0=8bit, 1=4bit, 2=16bit
54 : 8 2 Int16 Number of Bands.
55 : 10 6 char Unknown.
56 : 16 4 Float32 Width
57 : 20 4 Float32 Height
58 : 24 4 Int32 X Start (offset in original file?)
59 : 28 4 Int32 Y Start (offset in original file?)
60 : 32 56 char Unknown.
61 : 88 2 Int16 0=LAT, 1=UTM, 2=StatePlane, 3- are projections?
62 : 90 2 Int16 Classes in coverage.
63 : 92 14 char Unknown.
64 : 106 2 Int16 Area Unit (0=none, 1=Acre, 2=Hectare, 3=Other)
65 : 108 4 Float32 Pixel area.
66 : 112 4 Float32 Upper Left corner X (center of pixel?)
67 : 116 4 Float32 Upper Left corner Y (center of pixel?)
68 : 120 4 Float32 Width of a pixel.
69 : 124 4 Float32 Height of a pixel.
70 :
71 : All binary fields are in the same byte order but it may be big endian or
72 : little endian depending on what platform the file was written on. Usually
73 : this can be checked against the number of bands though this test won't work
74 : if there are more than 255 bands.
75 :
76 : There is also some information on .STA and .TRL files at:
77 :
78 : http://www.pcigeomatics.com/cgi-bin/pcihlp/ERDASWR%7CTRAILER+FORMAT
79 :
80 : **/
81 :
82 : constexpr int ERD_HEADER_SIZE = 128;
83 :
84 : /************************************************************************/
85 : /* ==================================================================== */
86 : /* LAN4BitRasterBand */
87 : /* ==================================================================== */
88 : /************************************************************************/
89 :
90 : class LANDataset;
91 :
92 : class LAN4BitRasterBand final : public GDALPamRasterBand
93 : {
94 : GDALColorTable *poCT;
95 : GDALColorInterp eInterp;
96 :
97 : CPL_DISALLOW_COPY_ASSIGN(LAN4BitRasterBand)
98 :
99 : public:
100 : LAN4BitRasterBand(LANDataset *, int);
101 : ~LAN4BitRasterBand() override;
102 :
103 : GDALColorTable *GetColorTable() override;
104 : GDALColorInterp GetColorInterpretation() override;
105 : CPLErr SetColorTable(GDALColorTable *) override;
106 : CPLErr SetColorInterpretation(GDALColorInterp) override;
107 :
108 : CPLErr IReadBlock(int, int, void *) override;
109 : };
110 :
111 : /************************************************************************/
112 : /* ==================================================================== */
113 : /* LANDataset */
114 : /* ==================================================================== */
115 : /************************************************************************/
116 :
117 : class LANDataset final : public RawDataset
118 : {
119 : CPL_DISALLOW_COPY_ASSIGN(LANDataset)
120 :
121 : public:
122 : VSILFILE *fpImage; // Image data file.
123 :
124 : char pachHeader[ERD_HEADER_SIZE];
125 :
126 : OGRSpatialReference *m_poSRS = nullptr;
127 :
128 : double adfGeoTransform[6];
129 :
130 : CPLString osSTAFilename{};
131 : void CheckForStatistics(void);
132 :
133 : char **GetFileList() override;
134 :
135 : CPLErr Close() override;
136 :
137 : public:
138 : LANDataset();
139 : ~LANDataset() override;
140 :
141 : CPLErr GetGeoTransform(double *padfTransform) override;
142 : CPLErr SetGeoTransform(double *padfTransform) override;
143 :
144 : const OGRSpatialReference *GetSpatialRef() const override;
145 : CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
146 :
147 : static GDALDataset *Open(GDALOpenInfo *);
148 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
149 : int nBandsIn, GDALDataType eType,
150 : char **papszOptions);
151 : };
152 :
153 : /************************************************************************/
154 : /* ==================================================================== */
155 : /* LAN4BitRasterBand */
156 : /* ==================================================================== */
157 : /************************************************************************/
158 :
159 : /************************************************************************/
160 : /* LAN4BitRasterBand() */
161 : /************************************************************************/
162 :
163 2 : LAN4BitRasterBand::LAN4BitRasterBand(LANDataset *poDSIn, int nBandIn)
164 2 : : poCT(nullptr), eInterp(GCI_Undefined)
165 : {
166 2 : poDS = poDSIn;
167 2 : nBand = nBandIn;
168 2 : eDataType = GDT_Byte;
169 :
170 2 : nBlockXSize = poDSIn->GetRasterXSize();
171 2 : nBlockYSize = 1;
172 2 : }
173 :
174 : /************************************************************************/
175 : /* ~LAN4BitRasterBand() */
176 : /************************************************************************/
177 :
178 4 : LAN4BitRasterBand::~LAN4BitRasterBand()
179 :
180 : {
181 2 : if (poCT)
182 0 : delete poCT;
183 4 : }
184 :
185 : /************************************************************************/
186 : /* IReadBlock() */
187 : /************************************************************************/
188 :
189 2 : CPLErr LAN4BitRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
190 : void *pImage)
191 :
192 : {
193 2 : LANDataset *poLAN_DS = reinterpret_cast<LANDataset *>(poDS);
194 2 : CPLAssert(nBlockXOff == 0);
195 :
196 : /* -------------------------------------------------------------------- */
197 : /* Seek to profile. */
198 : /* -------------------------------------------------------------------- */
199 : const vsi_l_offset nOffset =
200 : ERD_HEADER_SIZE +
201 4 : (static_cast<vsi_l_offset>(nBlockYOff) * nRasterXSize *
202 2 : poLAN_DS->GetRasterCount()) /
203 2 : 2 +
204 2 : (static_cast<vsi_l_offset>(nBand - 1) * nRasterXSize) / 2;
205 :
206 2 : if (VSIFSeekL(poLAN_DS->fpImage, nOffset, SEEK_SET) != 0)
207 : {
208 0 : CPLError(CE_Failure, CPLE_FileIO, "LAN Seek failed:%s",
209 0 : VSIStrerror(errno));
210 0 : return CE_Failure;
211 : }
212 :
213 : /* -------------------------------------------------------------------- */
214 : /* Read the profile. */
215 : /* -------------------------------------------------------------------- */
216 2 : if (VSIFReadL(pImage, 1, nRasterXSize / 2, poLAN_DS->fpImage) !=
217 2 : static_cast<size_t>(nRasterXSize) / 2)
218 : {
219 0 : CPLError(CE_Failure, CPLE_FileIO, "LAN Read failed:%s",
220 0 : VSIStrerror(errno));
221 0 : return CE_Failure;
222 : }
223 :
224 : /* -------------------------------------------------------------------- */
225 : /* Convert 4bit to 8bit. */
226 : /* -------------------------------------------------------------------- */
227 6 : for (int i = nRasterXSize - 1; i >= 0; i--)
228 : {
229 4 : if ((i & 0x01) != 0)
230 2 : reinterpret_cast<GByte *>(pImage)[i] =
231 2 : reinterpret_cast<GByte *>(pImage)[i / 2] & 0x0f;
232 : else
233 2 : reinterpret_cast<GByte *>(pImage)[i] =
234 2 : (reinterpret_cast<GByte *>(pImage)[i / 2] & 0xf0) / 16;
235 : }
236 :
237 2 : return CE_None;
238 : }
239 :
240 : /************************************************************************/
241 : /* SetColorTable() */
242 : /************************************************************************/
243 :
244 0 : CPLErr LAN4BitRasterBand::SetColorTable(GDALColorTable *poNewCT)
245 :
246 : {
247 0 : if (poCT)
248 0 : delete poCT;
249 0 : if (poNewCT == nullptr)
250 0 : poCT = nullptr;
251 : else
252 0 : poCT = poNewCT->Clone();
253 :
254 0 : return CE_None;
255 : }
256 :
257 : /************************************************************************/
258 : /* GetColorTable() */
259 : /************************************************************************/
260 :
261 0 : GDALColorTable *LAN4BitRasterBand::GetColorTable()
262 :
263 : {
264 0 : if (poCT != nullptr)
265 0 : return poCT;
266 :
267 0 : return GDALPamRasterBand::GetColorTable();
268 : }
269 :
270 : /************************************************************************/
271 : /* SetColorInterpretation() */
272 : /************************************************************************/
273 :
274 0 : CPLErr LAN4BitRasterBand::SetColorInterpretation(GDALColorInterp eNewInterp)
275 :
276 : {
277 0 : eInterp = eNewInterp;
278 :
279 0 : return CE_None;
280 : }
281 :
282 : /************************************************************************/
283 : /* GetColorInterpretation() */
284 : /************************************************************************/
285 :
286 0 : GDALColorInterp LAN4BitRasterBand::GetColorInterpretation()
287 :
288 : {
289 0 : return eInterp;
290 : }
291 :
292 : /************************************************************************/
293 : /* ==================================================================== */
294 : /* LANDataset */
295 : /* ==================================================================== */
296 : /************************************************************************/
297 :
298 : /************************************************************************/
299 : /* LANDataset() */
300 : /************************************************************************/
301 :
302 26 : LANDataset::LANDataset() : fpImage(nullptr)
303 : {
304 26 : memset(pachHeader, 0, sizeof(pachHeader));
305 26 : adfGeoTransform[0] = 0.0;
306 26 : adfGeoTransform[1] = 0.0; // TODO(schwehr): Should this be 1.0?
307 26 : adfGeoTransform[2] = 0.0;
308 26 : adfGeoTransform[3] = 0.0;
309 26 : adfGeoTransform[4] = 0.0;
310 26 : adfGeoTransform[5] = 0.0; // TODO(schwehr): Should this be 1.0?
311 26 : }
312 :
313 : /************************************************************************/
314 : /* ~LANDataset() */
315 : /************************************************************************/
316 :
317 52 : LANDataset::~LANDataset()
318 :
319 : {
320 26 : LANDataset::Close();
321 52 : }
322 :
323 : /************************************************************************/
324 : /* Close() */
325 : /************************************************************************/
326 :
327 50 : CPLErr LANDataset::Close()
328 : {
329 50 : CPLErr eErr = CE_None;
330 50 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
331 : {
332 26 : if (LANDataset::FlushCache(true) != CE_None)
333 0 : eErr = CE_Failure;
334 :
335 26 : if (fpImage)
336 : {
337 26 : if (VSIFCloseL(fpImage) != 0)
338 : {
339 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
340 0 : eErr = CE_Failure;
341 : }
342 : }
343 :
344 26 : if (m_poSRS)
345 24 : m_poSRS->Release();
346 :
347 26 : if (GDALPamDataset::Close() != CE_None)
348 0 : eErr = CE_Failure;
349 : }
350 50 : return eErr;
351 : }
352 :
353 : /************************************************************************/
354 : /* Open() */
355 : /************************************************************************/
356 :
357 30959 : GDALDataset *LANDataset::Open(GDALOpenInfo *poOpenInfo)
358 :
359 : {
360 : /* -------------------------------------------------------------------- */
361 : /* We assume the user is pointing to the header (.pcb) file. */
362 : /* Does this appear to be a pcb file? */
363 : /* -------------------------------------------------------------------- */
364 30959 : if (poOpenInfo->nHeaderBytes < ERD_HEADER_SIZE ||
365 3094 : poOpenInfo->fpL == nullptr)
366 27903 : return nullptr;
367 :
368 3056 : if (!STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
369 3056 : "HEADER") &&
370 3056 : !STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
371 : "HEAD74"))
372 3030 : return nullptr;
373 :
374 26 : if (memcmp(poOpenInfo->pabyHeader + 16, "S LAT ", 8) == 0)
375 : {
376 : // NTV1 format
377 0 : return nullptr;
378 : }
379 :
380 : /* -------------------------------------------------------------------- */
381 : /* Create a corresponding GDALDataset. */
382 : /* -------------------------------------------------------------------- */
383 52 : auto poDS = std::make_unique<LANDataset>();
384 :
385 26 : poDS->eAccess = poOpenInfo->eAccess;
386 26 : std::swap(poDS->fpImage, poOpenInfo->fpL);
387 :
388 : /* -------------------------------------------------------------------- */
389 : /* Do we need to byte swap the headers to local machine order? */
390 : /* -------------------------------------------------------------------- */
391 26 : const RawRasterBand::ByteOrder eByteOrder =
392 26 : poOpenInfo->pabyHeader[8] == 0
393 26 : ? RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN
394 : : RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN;
395 :
396 26 : memcpy(poDS->pachHeader, poOpenInfo->pabyHeader, ERD_HEADER_SIZE);
397 :
398 26 : if (eByteOrder != RawRasterBand::NATIVE_BYTE_ORDER)
399 : {
400 2 : CPL_SWAP16PTR(poDS->pachHeader + 6);
401 2 : CPL_SWAP16PTR(poDS->pachHeader + 8);
402 :
403 2 : CPL_SWAP32PTR(poDS->pachHeader + 16);
404 2 : CPL_SWAP32PTR(poDS->pachHeader + 20);
405 2 : CPL_SWAP32PTR(poDS->pachHeader + 24);
406 2 : CPL_SWAP32PTR(poDS->pachHeader + 28);
407 :
408 2 : CPL_SWAP16PTR(poDS->pachHeader + 88);
409 2 : CPL_SWAP16PTR(poDS->pachHeader + 90);
410 :
411 2 : CPL_SWAP16PTR(poDS->pachHeader + 106);
412 2 : CPL_SWAP32PTR(poDS->pachHeader + 108);
413 2 : CPL_SWAP32PTR(poDS->pachHeader + 112);
414 2 : CPL_SWAP32PTR(poDS->pachHeader + 116);
415 2 : CPL_SWAP32PTR(poDS->pachHeader + 120);
416 2 : CPL_SWAP32PTR(poDS->pachHeader + 124);
417 : }
418 :
419 : /* -------------------------------------------------------------------- */
420 : /* Capture some information from the file that is of interest. */
421 : /* -------------------------------------------------------------------- */
422 26 : if (STARTS_WITH_CI(poDS->pachHeader, "HEADER"))
423 : {
424 0 : float fTmp = 0.0;
425 0 : memcpy(&fTmp, poDS->pachHeader + 16, 4);
426 0 : poDS->nRasterXSize = static_cast<int>(fTmp);
427 0 : memcpy(&fTmp, poDS->pachHeader + 20, 4);
428 0 : poDS->nRasterYSize = static_cast<int>(fTmp);
429 : }
430 : else
431 : {
432 26 : GInt32 nTmp = 0;
433 26 : memcpy(&nTmp, poDS->pachHeader + 16, 4);
434 26 : poDS->nRasterXSize = nTmp;
435 26 : memcpy(&nTmp, poDS->pachHeader + 20, 4);
436 26 : poDS->nRasterYSize = nTmp;
437 : }
438 :
439 26 : GInt16 nTmp16 = 0;
440 26 : memcpy(&nTmp16, poDS->pachHeader + 6, 2);
441 :
442 26 : int nPixelOffset = 0;
443 26 : GDALDataType eDataType = GDT_Unknown;
444 26 : if (nTmp16 == 0)
445 : {
446 19 : eDataType = GDT_Byte;
447 19 : nPixelOffset = 1;
448 : }
449 7 : else if (nTmp16 == 1) // 4 bit
450 : {
451 2 : eDataType = GDT_Byte;
452 2 : nPixelOffset = -1;
453 : }
454 5 : else if (nTmp16 == 2)
455 : {
456 5 : nPixelOffset = 2;
457 5 : eDataType = GDT_Int16;
458 : }
459 : else
460 : {
461 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unsupported pixel type (%d).",
462 : nTmp16);
463 0 : return nullptr;
464 : }
465 :
466 26 : memcpy(&nTmp16, poDS->pachHeader + 8, 2);
467 26 : const int nBandCount = nTmp16;
468 :
469 52 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
470 26 : !GDALCheckBandCount(nBandCount, FALSE))
471 : {
472 2 : return nullptr;
473 : }
474 :
475 : // cppcheck-suppress knownConditionTrueFalse
476 46 : if (nPixelOffset != -1 &&
477 22 : poDS->nRasterXSize > INT_MAX / (nPixelOffset * nBandCount))
478 : {
479 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow occurred.");
480 0 : return nullptr;
481 : }
482 :
483 : /* -------------------------------------------------------------------- */
484 : /* Create band information object. */
485 : /* -------------------------------------------------------------------- */
486 82 : for (int iBand = 1; iBand <= nBandCount; iBand++)
487 : {
488 58 : if (nPixelOffset == -1) /* 4 bit case */
489 2 : poDS->SetBand(iBand, new LAN4BitRasterBand(poDS.get(), iBand));
490 : else
491 : {
492 : auto poBand = RawRasterBand::Create(
493 112 : poDS.get(), iBand, poDS->fpImage,
494 56 : ERD_HEADER_SIZE +
495 112 : (iBand - 1) * nPixelOffset * poDS->nRasterXSize,
496 56 : nPixelOffset, poDS->nRasterXSize * nPixelOffset * nBandCount,
497 56 : eDataType, eByteOrder, RawRasterBand::OwnFP::NO);
498 56 : if (!poBand)
499 0 : return nullptr;
500 56 : poDS->SetBand(iBand, std::move(poBand));
501 : }
502 : }
503 :
504 : /* -------------------------------------------------------------------- */
505 : /* Initialize any PAM information. */
506 : /* -------------------------------------------------------------------- */
507 24 : poDS->SetDescription(poOpenInfo->pszFilename);
508 24 : poDS->CheckForStatistics();
509 24 : poDS->TryLoadXML();
510 :
511 : /* -------------------------------------------------------------------- */
512 : /* Check for overviews. */
513 : /* -------------------------------------------------------------------- */
514 24 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
515 :
516 : /* -------------------------------------------------------------------- */
517 : /* Try to interpret georeferencing. */
518 : /* -------------------------------------------------------------------- */
519 24 : float fTmp = 0.0;
520 :
521 24 : memcpy(&fTmp, poDS->pachHeader + 112, 4);
522 24 : poDS->adfGeoTransform[0] = fTmp;
523 24 : memcpy(&fTmp, poDS->pachHeader + 120, 4);
524 24 : poDS->adfGeoTransform[1] = fTmp;
525 24 : poDS->adfGeoTransform[2] = 0.0;
526 24 : memcpy(&fTmp, poDS->pachHeader + 116, 4);
527 24 : poDS->adfGeoTransform[3] = fTmp;
528 24 : poDS->adfGeoTransform[4] = 0.0;
529 24 : memcpy(&fTmp, poDS->pachHeader + 124, 4);
530 24 : poDS->adfGeoTransform[5] = -fTmp;
531 :
532 : // adjust for center of pixel vs. top left corner of pixel.
533 24 : poDS->adfGeoTransform[0] -= poDS->adfGeoTransform[1] * 0.5;
534 24 : poDS->adfGeoTransform[3] -= poDS->adfGeoTransform[5] * 0.5;
535 :
536 : /* -------------------------------------------------------------------- */
537 : /* If we didn't get any georeferencing, try for a worldfile. */
538 : /* -------------------------------------------------------------------- */
539 24 : if (poDS->adfGeoTransform[1] == 0.0 || poDS->adfGeoTransform[5] == 0.0)
540 : {
541 0 : if (!GDALReadWorldFile(poOpenInfo->pszFilename, nullptr,
542 0 : poDS->adfGeoTransform))
543 0 : GDALReadWorldFile(poOpenInfo->pszFilename, ".wld",
544 0 : poDS->adfGeoTransform);
545 : }
546 :
547 : /* -------------------------------------------------------------------- */
548 : /* Try to come up with something for the coordinate system. */
549 : /* -------------------------------------------------------------------- */
550 24 : memcpy(&nTmp16, poDS->pachHeader + 88, 2);
551 24 : int nCoordSys = nTmp16;
552 :
553 24 : poDS->m_poSRS = new OGRSpatialReference();
554 24 : poDS->m_poSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
555 24 : if (nCoordSys == 0)
556 : {
557 24 : poDS->m_poSRS->SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
558 : }
559 0 : else if (nCoordSys == 1)
560 : {
561 0 : poDS->m_poSRS->SetFromUserInput(
562 : "LOCAL_CS[\"UTM - Zone Unknown\",UNIT[\"Meter\",1]]");
563 : }
564 0 : else if (nCoordSys == 2)
565 : {
566 0 : poDS->m_poSRS->SetFromUserInput(
567 : "LOCAL_CS[\"State Plane - Zone Unknown\","
568 : "UNIT[\"US survey foot\",0.3048006096012192]]");
569 : }
570 : else
571 : {
572 0 : poDS->m_poSRS->SetFromUserInput(
573 : "LOCAL_CS[\"Unknown\",UNIT[\"Meter\",1]]");
574 : }
575 :
576 : /* -------------------------------------------------------------------- */
577 : /* Check for a trailer file with a colormap in it. */
578 : /* -------------------------------------------------------------------- */
579 24 : char *pszPath = CPLStrdup(CPLGetPathSafe(poOpenInfo->pszFilename).c_str());
580 : char *pszBasename =
581 24 : CPLStrdup(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str());
582 : const std::string osTRLFilename =
583 48 : CPLFormCIFilenameSafe(pszPath, pszBasename, "trl");
584 24 : VSILFILE *fpTRL = VSIFOpenL(osTRLFilename.c_str(), "rb");
585 24 : if (fpTRL != nullptr)
586 : {
587 0 : char szTRLData[896] = {'\0'};
588 :
589 0 : CPL_IGNORE_RET_VAL(VSIFReadL(szTRLData, 1, 896, fpTRL));
590 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpTRL));
591 :
592 0 : GDALColorTable oCT;
593 0 : for (int iColor = 0; iColor < 256; iColor++)
594 : {
595 0 : GDALColorEntry sEntry = {0, 0, 0, 0};
596 :
597 0 : sEntry.c2 = reinterpret_cast<GByte *>(szTRLData)[iColor + 128];
598 0 : sEntry.c1 =
599 0 : reinterpret_cast<GByte *>(szTRLData)[iColor + 128 + 256];
600 0 : sEntry.c3 =
601 0 : reinterpret_cast<GByte *>(szTRLData)[iColor + 128 + 512];
602 0 : sEntry.c4 = 255;
603 0 : oCT.SetColorEntry(iColor, &sEntry);
604 :
605 : // Only 16 colors in 4bit files.
606 0 : if (nPixelOffset == -1 && iColor == 15)
607 0 : break;
608 : }
609 :
610 0 : poDS->GetRasterBand(1)->SetColorTable(&oCT);
611 0 : poDS->GetRasterBand(1)->SetColorInterpretation(GCI_PaletteIndex);
612 : }
613 :
614 24 : CPLFree(pszPath);
615 24 : CPLFree(pszBasename);
616 :
617 24 : return poDS.release();
618 : }
619 :
620 : /************************************************************************/
621 : /* GetGeoTransform() */
622 : /************************************************************************/
623 :
624 7 : CPLErr LANDataset::GetGeoTransform(double *padfTransform)
625 :
626 : {
627 7 : if (adfGeoTransform[1] != 0.0 && adfGeoTransform[5] != 0.0)
628 : {
629 7 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
630 7 : return CE_None;
631 : }
632 :
633 0 : return GDALPamDataset::GetGeoTransform(padfTransform);
634 : }
635 :
636 : /************************************************************************/
637 : /* SetGeoTransform() */
638 : /************************************************************************/
639 :
640 13 : CPLErr LANDataset::SetGeoTransform(double *padfTransform)
641 :
642 : {
643 13 : unsigned char abyHeader[128] = {'\0'};
644 :
645 13 : memcpy(adfGeoTransform, padfTransform, sizeof(double) * 6);
646 :
647 13 : CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));
648 13 : CPL_IGNORE_RET_VAL(VSIFReadL(abyHeader, 128, 1, fpImage));
649 :
650 : // Upper Left X.
651 13 : float f32Val =
652 13 : static_cast<float>(adfGeoTransform[0] + 0.5 * adfGeoTransform[1]);
653 13 : memcpy(abyHeader + 112, &f32Val, 4);
654 :
655 : // Upper Left Y.
656 13 : f32Val = static_cast<float>(adfGeoTransform[3] + 0.5 * adfGeoTransform[5]);
657 13 : memcpy(abyHeader + 116, &f32Val, 4);
658 :
659 : // Width of pixel.
660 13 : f32Val = static_cast<float>(adfGeoTransform[1]);
661 13 : memcpy(abyHeader + 120, &f32Val, 4);
662 :
663 : // Height of pixel.
664 13 : f32Val = static_cast<float>(std::abs(adfGeoTransform[5]));
665 13 : memcpy(abyHeader + 124, &f32Val, 4);
666 :
667 26 : if (VSIFSeekL(fpImage, 0, SEEK_SET) != 0 ||
668 13 : VSIFWriteL(abyHeader, 128, 1, fpImage) != 1)
669 : {
670 0 : CPLError(CE_Failure, CPLE_FileIO,
671 : "File IO Error writing header with new geotransform.");
672 0 : return CE_Failure;
673 : }
674 :
675 13 : return CE_None;
676 : }
677 :
678 : /************************************************************************/
679 : /* GetSpatialRef() */
680 : /* */
681 : /* Use PAM coordinate system if available in preference to the */
682 : /* generally poor value derived from the file itself. */
683 : /************************************************************************/
684 :
685 0 : const OGRSpatialReference *LANDataset::GetSpatialRef() const
686 :
687 : {
688 0 : const auto poSRS = GDALPamDataset::GetSpatialRef();
689 0 : if (poSRS)
690 0 : return poSRS;
691 :
692 0 : return m_poSRS;
693 : }
694 :
695 : /************************************************************************/
696 : /* SetSpatialRef() */
697 : /************************************************************************/
698 :
699 13 : CPLErr LANDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
700 :
701 : {
702 13 : if (poSRS == nullptr)
703 0 : return GDALPamDataset::SetSpatialRef(poSRS);
704 :
705 13 : unsigned char abyHeader[128] = {'\0'};
706 :
707 13 : CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));
708 13 : CPL_IGNORE_RET_VAL(VSIFReadL(abyHeader, 128, 1, fpImage));
709 :
710 13 : GUInt16 nProjCode = 0;
711 :
712 13 : if (poSRS->IsGeographic())
713 : {
714 13 : nProjCode = 0;
715 : }
716 0 : else if (poSRS->GetUTMZone() != 0)
717 : {
718 0 : nProjCode = 1;
719 : }
720 : // Too bad we have no way of recognising state plane projections.
721 : else
722 : {
723 0 : const char *l_pszProjection = poSRS->GetAttrValue("PROJECTION");
724 :
725 0 : if (l_pszProjection == nullptr)
726 : ;
727 0 : else if (EQUAL(l_pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA))
728 0 : nProjCode = 3;
729 0 : else if (EQUAL(l_pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
730 0 : nProjCode = 4;
731 0 : else if (EQUAL(l_pszProjection, SRS_PT_MERCATOR_1SP))
732 0 : nProjCode = 5;
733 0 : else if (EQUAL(l_pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
734 0 : nProjCode = 6;
735 0 : else if (EQUAL(l_pszProjection, SRS_PT_POLYCONIC))
736 0 : nProjCode = 7;
737 0 : else if (EQUAL(l_pszProjection, SRS_PT_EQUIDISTANT_CONIC))
738 0 : nProjCode = 8;
739 0 : else if (EQUAL(l_pszProjection, SRS_PT_TRANSVERSE_MERCATOR))
740 0 : nProjCode = 9;
741 0 : else if (EQUAL(l_pszProjection, SRS_PT_STEREOGRAPHIC))
742 0 : nProjCode = 10;
743 0 : else if (EQUAL(l_pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
744 0 : nProjCode = 11;
745 0 : else if (EQUAL(l_pszProjection, SRS_PT_AZIMUTHAL_EQUIDISTANT))
746 0 : nProjCode = 12;
747 0 : else if (EQUAL(l_pszProjection, SRS_PT_GNOMONIC))
748 0 : nProjCode = 13;
749 0 : else if (EQUAL(l_pszProjection, SRS_PT_ORTHOGRAPHIC))
750 0 : nProjCode = 14;
751 : // We do not have GVNP.
752 0 : else if (EQUAL(l_pszProjection, SRS_PT_SINUSOIDAL))
753 0 : nProjCode = 16;
754 0 : else if (EQUAL(l_pszProjection, SRS_PT_EQUIRECTANGULAR))
755 0 : nProjCode = 17;
756 0 : else if (EQUAL(l_pszProjection, SRS_PT_MILLER_CYLINDRICAL))
757 0 : nProjCode = 18;
758 0 : else if (EQUAL(l_pszProjection, SRS_PT_VANDERGRINTEN))
759 0 : nProjCode = 19;
760 0 : else if (EQUAL(l_pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR))
761 0 : nProjCode = 20;
762 : }
763 :
764 13 : memcpy(abyHeader + 88, &nProjCode, 2);
765 :
766 13 : CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));
767 13 : CPL_IGNORE_RET_VAL(VSIFWriteL(abyHeader, 128, 1, fpImage));
768 :
769 13 : return GDALPamDataset::SetSpatialRef(poSRS);
770 : }
771 :
772 : /************************************************************************/
773 : /* GetFileList() */
774 : /************************************************************************/
775 :
776 2 : char **LANDataset::GetFileList()
777 :
778 : {
779 : // Main data file, etc.
780 2 : char **papszFileList = GDALPamDataset::GetFileList();
781 :
782 2 : if (!osSTAFilename.empty())
783 0 : papszFileList = CSLAddString(papszFileList, osSTAFilename);
784 :
785 2 : return papszFileList;
786 : }
787 :
788 : /************************************************************************/
789 : /* CheckForStatistics() */
790 : /************************************************************************/
791 :
792 24 : void LANDataset::CheckForStatistics()
793 :
794 : {
795 : /* -------------------------------------------------------------------- */
796 : /* Do we have a statistics file? */
797 : /* -------------------------------------------------------------------- */
798 24 : osSTAFilename = CPLResetExtensionSafe(GetDescription(), "sta");
799 :
800 24 : VSILFILE *fpSTA = VSIFOpenL(osSTAFilename, "r");
801 :
802 24 : if (fpSTA == nullptr && VSIIsCaseSensitiveFS(osSTAFilename))
803 : {
804 24 : osSTAFilename = CPLResetExtensionSafe(GetDescription(), "STA");
805 24 : fpSTA = VSIFOpenL(osSTAFilename, "r");
806 : }
807 :
808 24 : if (fpSTA == nullptr)
809 : {
810 24 : osSTAFilename = "";
811 24 : return;
812 : }
813 :
814 : /* -------------------------------------------------------------------- */
815 : /* Read it one band at a time. */
816 : /* -------------------------------------------------------------------- */
817 0 : GByte abyBandInfo[1152] = {'\0'};
818 :
819 0 : for (int iBand = 0; iBand < nBands; iBand++)
820 : {
821 0 : if (VSIFReadL(abyBandInfo, 1152, 1, fpSTA) != 1)
822 0 : break;
823 :
824 0 : const int nBandNumber = abyBandInfo[7];
825 0 : GDALRasterBand *poBand = GetRasterBand(nBandNumber);
826 0 : if (poBand == nullptr)
827 0 : break;
828 :
829 0 : GInt16 nMin = 0;
830 0 : GInt16 nMax = 0;
831 :
832 0 : if (poBand->GetRasterDataType() != GDT_Byte)
833 : {
834 0 : memcpy(&nMin, abyBandInfo + 28, 2);
835 0 : memcpy(&nMax, abyBandInfo + 30, 2);
836 0 : CPL_LSBPTR16(&nMin);
837 0 : CPL_LSBPTR16(&nMax);
838 : }
839 : else
840 : {
841 0 : nMin = abyBandInfo[9];
842 0 : nMax = abyBandInfo[8];
843 : }
844 :
845 0 : float fMean = 0.0;
846 0 : float fStdDev = 0.0;
847 0 : memcpy(&fMean, abyBandInfo + 12, 4);
848 0 : memcpy(&fStdDev, abyBandInfo + 24, 4);
849 0 : CPL_LSBPTR32(&fMean);
850 0 : CPL_LSBPTR32(&fStdDev);
851 :
852 0 : poBand->SetStatistics(nMin, nMax, fMean, fStdDev);
853 : }
854 :
855 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpSTA));
856 : }
857 :
858 : /************************************************************************/
859 : /* Create() */
860 : /************************************************************************/
861 :
862 60 : GDALDataset *LANDataset::Create(const char *pszFilename, int nXSize, int nYSize,
863 : int nBandsIn, GDALDataType eType,
864 : char ** /* papszOptions */)
865 : {
866 60 : if (eType != GDT_Byte && eType != GDT_Int16)
867 : {
868 33 : CPLError(CE_Failure, CPLE_AppDefined,
869 : "Attempt to create .GIS file with unsupported data type '%s'.",
870 : GDALGetDataTypeName(eType));
871 33 : return nullptr;
872 : }
873 :
874 : /* -------------------------------------------------------------------- */
875 : /* Try to create the file. */
876 : /* -------------------------------------------------------------------- */
877 27 : VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
878 27 : if (fp == nullptr)
879 : {
880 3 : CPLError(CE_Failure, CPLE_OpenFailed,
881 : "Attempt to create file `%s' failed.\n", pszFilename);
882 3 : return nullptr;
883 : }
884 :
885 : /* -------------------------------------------------------------------- */
886 : /* Write out the header. */
887 : /* -------------------------------------------------------------------- */
888 24 : unsigned char abyHeader[128] = {'\0'};
889 :
890 24 : memset(abyHeader, 0, sizeof(abyHeader));
891 :
892 24 : memcpy(abyHeader + 0, "HEAD74", 6);
893 :
894 : // Pixel type.
895 24 : GInt16 n16Val = 0;
896 24 : if (eType == GDT_Byte) // Do we want 4bit?
897 21 : n16Val = 0;
898 : else
899 3 : n16Val = 2;
900 24 : memcpy(abyHeader + 6, &n16Val, 2);
901 :
902 : // Number of Bands.
903 24 : n16Val = static_cast<GInt16>(nBandsIn);
904 24 : memcpy(abyHeader + 8, &n16Val, 2);
905 :
906 : // Unknown (6).
907 :
908 : // Width.
909 24 : GInt32 n32Val = nXSize;
910 24 : memcpy(abyHeader + 16, &n32Val, 4);
911 :
912 : // Height.
913 24 : n32Val = nYSize;
914 24 : memcpy(abyHeader + 20, &n32Val, 4);
915 :
916 : // X Start (4).
917 : // Y Start (4).
918 :
919 : // Unknown (56).
920 :
921 : // Coordinate System.
922 24 : n16Val = 0;
923 24 : memcpy(abyHeader + 88, &n16Val, 2);
924 :
925 : // Classes in coverage.
926 24 : n16Val = 0;
927 24 : memcpy(abyHeader + 90, &n16Val, 2);
928 :
929 : // Unknown (14).
930 :
931 : // Area Unit.
932 24 : n16Val = 0;
933 24 : memcpy(abyHeader + 106, &n16Val, 2);
934 :
935 : // Pixel Area.
936 24 : float f32Val = 0.0f;
937 24 : memcpy(abyHeader + 108, &f32Val, 4);
938 :
939 : // Upper Left X.
940 24 : f32Val = 0.5f;
941 24 : memcpy(abyHeader + 112, &f32Val, 4);
942 :
943 : // Upper Left Y
944 24 : f32Val = static_cast<float>(nYSize - 0.5);
945 24 : memcpy(abyHeader + 116, &f32Val, 4);
946 :
947 : // Width of pixel.
948 24 : f32Val = 1.0f;
949 24 : memcpy(abyHeader + 120, &f32Val, 4);
950 :
951 : // Height of pixel.
952 24 : f32Val = 1.0f;
953 24 : memcpy(abyHeader + 124, &f32Val, 4);
954 :
955 24 : CPL_IGNORE_RET_VAL(VSIFWriteL(abyHeader, sizeof(abyHeader), 1, fp));
956 :
957 : /* -------------------------------------------------------------------- */
958 : /* Extend the file to the target size. */
959 : /* -------------------------------------------------------------------- */
960 24 : vsi_l_offset nImageBytes = 0;
961 :
962 24 : if (eType != GDT_Byte)
963 3 : nImageBytes = nXSize * static_cast<vsi_l_offset>(nYSize) * 2;
964 : else
965 21 : nImageBytes = nXSize * static_cast<vsi_l_offset>(nYSize);
966 :
967 24 : memset(abyHeader, 0, sizeof(abyHeader));
968 :
969 820 : while (nImageBytes > 0)
970 : {
971 : const vsi_l_offset nWriteThisTime =
972 805 : std::min(static_cast<size_t>(nImageBytes), sizeof(abyHeader));
973 :
974 805 : if (VSIFWriteL(abyHeader, 1, static_cast<size_t>(nWriteThisTime), fp) !=
975 : nWriteThisTime)
976 : {
977 9 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
978 9 : CPLError(CE_Failure, CPLE_FileIO,
979 : "Failed to write whole Istar file.");
980 9 : return nullptr;
981 : }
982 796 : nImageBytes -= nWriteThisTime;
983 : }
984 :
985 15 : if (VSIFCloseL(fp) != 0)
986 : {
987 0 : CPLError(CE_Failure, CPLE_FileIO, "Failed to write whole Istar file.");
988 0 : return nullptr;
989 : }
990 :
991 15 : return static_cast<GDALDataset *>(GDALOpen(pszFilename, GA_Update));
992 : }
993 :
994 : /************************************************************************/
995 : /* GDALRegister_LAN() */
996 : /************************************************************************/
997 :
998 1682 : void GDALRegister_LAN()
999 :
1000 : {
1001 1682 : if (GDALGetDriverByName("LAN") != nullptr)
1002 301 : return;
1003 :
1004 1381 : GDALDriver *poDriver = new GDALDriver();
1005 :
1006 1381 : poDriver->SetDescription("LAN");
1007 1381 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1008 1381 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Erdas .LAN/.GIS");
1009 1381 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/lan.html");
1010 1381 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1011 1381 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Byte Int16");
1012 :
1013 1381 : poDriver->pfnOpen = LANDataset::Open;
1014 1381 : poDriver->pfnCreate = LANDataset::Create;
1015 :
1016 1381 : GetGDALDriverManager()->RegisterDriver(poDriver);
1017 : }
|
