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 30620 : 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 30620 : if (poOpenInfo->nHeaderBytes < ERD_HEADER_SIZE ||
365 3021 : poOpenInfo->fpL == nullptr)
366 27637 : return nullptr;
367 :
368 2983 : if (!STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
369 2983 : "HEADER") &&
370 2983 : !STARTS_WITH_CI(reinterpret_cast<char *>(poOpenInfo->pabyHeader),
371 : "HEAD74"))
372 2957 : 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(CPLGetPath(poOpenInfo->pszFilename));
580 24 : char *pszBasename = CPLStrdup(CPLGetBasename(poOpenInfo->pszFilename));
581 24 : const char *pszTRLFilename = CPLFormCIFilename(pszPath, pszBasename, "trl");
582 24 : VSILFILE *fpTRL = VSIFOpenL(pszTRLFilename, "rb");
583 24 : if (fpTRL != nullptr)
584 : {
585 0 : char szTRLData[896] = {'\0'};
586 :
587 0 : CPL_IGNORE_RET_VAL(VSIFReadL(szTRLData, 1, 896, fpTRL));
588 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpTRL));
589 :
590 0 : GDALColorTable oCT;
591 0 : for (int iColor = 0; iColor < 256; iColor++)
592 : {
593 0 : GDALColorEntry sEntry = {0, 0, 0, 0};
594 :
595 0 : sEntry.c2 = reinterpret_cast<GByte *>(szTRLData)[iColor + 128];
596 0 : sEntry.c1 =
597 0 : reinterpret_cast<GByte *>(szTRLData)[iColor + 128 + 256];
598 0 : sEntry.c3 =
599 0 : reinterpret_cast<GByte *>(szTRLData)[iColor + 128 + 512];
600 0 : sEntry.c4 = 255;
601 0 : oCT.SetColorEntry(iColor, &sEntry);
602 :
603 : // Only 16 colors in 4bit files.
604 0 : if (nPixelOffset == -1 && iColor == 15)
605 0 : break;
606 : }
607 :
608 0 : poDS->GetRasterBand(1)->SetColorTable(&oCT);
609 0 : poDS->GetRasterBand(1)->SetColorInterpretation(GCI_PaletteIndex);
610 : }
611 :
612 24 : CPLFree(pszPath);
613 24 : CPLFree(pszBasename);
614 :
615 24 : return poDS.release();
616 : }
617 :
618 : /************************************************************************/
619 : /* GetGeoTransform() */
620 : /************************************************************************/
621 :
622 7 : CPLErr LANDataset::GetGeoTransform(double *padfTransform)
623 :
624 : {
625 7 : if (adfGeoTransform[1] != 0.0 && adfGeoTransform[5] != 0.0)
626 : {
627 7 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
628 7 : return CE_None;
629 : }
630 :
631 0 : return GDALPamDataset::GetGeoTransform(padfTransform);
632 : }
633 :
634 : /************************************************************************/
635 : /* SetGeoTransform() */
636 : /************************************************************************/
637 :
638 13 : CPLErr LANDataset::SetGeoTransform(double *padfTransform)
639 :
640 : {
641 13 : unsigned char abyHeader[128] = {'\0'};
642 :
643 13 : memcpy(adfGeoTransform, padfTransform, sizeof(double) * 6);
644 :
645 13 : CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));
646 13 : CPL_IGNORE_RET_VAL(VSIFReadL(abyHeader, 128, 1, fpImage));
647 :
648 : // Upper Left X.
649 13 : float f32Val =
650 13 : static_cast<float>(adfGeoTransform[0] + 0.5 * adfGeoTransform[1]);
651 13 : memcpy(abyHeader + 112, &f32Val, 4);
652 :
653 : // Upper Left Y.
654 13 : f32Val = static_cast<float>(adfGeoTransform[3] + 0.5 * adfGeoTransform[5]);
655 13 : memcpy(abyHeader + 116, &f32Val, 4);
656 :
657 : // Width of pixel.
658 13 : f32Val = static_cast<float>(adfGeoTransform[1]);
659 13 : memcpy(abyHeader + 120, &f32Val, 4);
660 :
661 : // Height of pixel.
662 13 : f32Val = static_cast<float>(std::abs(adfGeoTransform[5]));
663 13 : memcpy(abyHeader + 124, &f32Val, 4);
664 :
665 26 : if (VSIFSeekL(fpImage, 0, SEEK_SET) != 0 ||
666 13 : VSIFWriteL(abyHeader, 128, 1, fpImage) != 1)
667 : {
668 0 : CPLError(CE_Failure, CPLE_FileIO,
669 : "File IO Error writing header with new geotransform.");
670 0 : return CE_Failure;
671 : }
672 :
673 13 : return CE_None;
674 : }
675 :
676 : /************************************************************************/
677 : /* GetSpatialRef() */
678 : /* */
679 : /* Use PAM coordinate system if available in preference to the */
680 : /* generally poor value derived from the file itself. */
681 : /************************************************************************/
682 :
683 0 : const OGRSpatialReference *LANDataset::GetSpatialRef() const
684 :
685 : {
686 0 : const auto poSRS = GDALPamDataset::GetSpatialRef();
687 0 : if (poSRS)
688 0 : return poSRS;
689 :
690 0 : return m_poSRS;
691 : }
692 :
693 : /************************************************************************/
694 : /* SetSpatialRef() */
695 : /************************************************************************/
696 :
697 13 : CPLErr LANDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
698 :
699 : {
700 13 : if (poSRS == nullptr)
701 0 : return GDALPamDataset::SetSpatialRef(poSRS);
702 :
703 13 : unsigned char abyHeader[128] = {'\0'};
704 :
705 13 : CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));
706 13 : CPL_IGNORE_RET_VAL(VSIFReadL(abyHeader, 128, 1, fpImage));
707 :
708 13 : GUInt16 nProjCode = 0;
709 :
710 13 : if (poSRS->IsGeographic())
711 : {
712 13 : nProjCode = 0;
713 : }
714 0 : else if (poSRS->GetUTMZone() != 0)
715 : {
716 0 : nProjCode = 1;
717 : }
718 : // Too bad we have no way of recognising state plane projections.
719 : else
720 : {
721 0 : const char *l_pszProjection = poSRS->GetAttrValue("PROJECTION");
722 :
723 0 : if (l_pszProjection == nullptr)
724 : ;
725 0 : else if (EQUAL(l_pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA))
726 0 : nProjCode = 3;
727 0 : else if (EQUAL(l_pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
728 0 : nProjCode = 4;
729 0 : else if (EQUAL(l_pszProjection, SRS_PT_MERCATOR_1SP))
730 0 : nProjCode = 5;
731 0 : else if (EQUAL(l_pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
732 0 : nProjCode = 6;
733 0 : else if (EQUAL(l_pszProjection, SRS_PT_POLYCONIC))
734 0 : nProjCode = 7;
735 0 : else if (EQUAL(l_pszProjection, SRS_PT_EQUIDISTANT_CONIC))
736 0 : nProjCode = 8;
737 0 : else if (EQUAL(l_pszProjection, SRS_PT_TRANSVERSE_MERCATOR))
738 0 : nProjCode = 9;
739 0 : else if (EQUAL(l_pszProjection, SRS_PT_STEREOGRAPHIC))
740 0 : nProjCode = 10;
741 0 : else if (EQUAL(l_pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
742 0 : nProjCode = 11;
743 0 : else if (EQUAL(l_pszProjection, SRS_PT_AZIMUTHAL_EQUIDISTANT))
744 0 : nProjCode = 12;
745 0 : else if (EQUAL(l_pszProjection, SRS_PT_GNOMONIC))
746 0 : nProjCode = 13;
747 0 : else if (EQUAL(l_pszProjection, SRS_PT_ORTHOGRAPHIC))
748 0 : nProjCode = 14;
749 : // We do not have GVNP.
750 0 : else if (EQUAL(l_pszProjection, SRS_PT_SINUSOIDAL))
751 0 : nProjCode = 16;
752 0 : else if (EQUAL(l_pszProjection, SRS_PT_EQUIRECTANGULAR))
753 0 : nProjCode = 17;
754 0 : else if (EQUAL(l_pszProjection, SRS_PT_MILLER_CYLINDRICAL))
755 0 : nProjCode = 18;
756 0 : else if (EQUAL(l_pszProjection, SRS_PT_VANDERGRINTEN))
757 0 : nProjCode = 19;
758 0 : else if (EQUAL(l_pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR))
759 0 : nProjCode = 20;
760 : }
761 :
762 13 : memcpy(abyHeader + 88, &nProjCode, 2);
763 :
764 13 : CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));
765 13 : CPL_IGNORE_RET_VAL(VSIFWriteL(abyHeader, 128, 1, fpImage));
766 :
767 13 : return GDALPamDataset::SetSpatialRef(poSRS);
768 : }
769 :
770 : /************************************************************************/
771 : /* GetFileList() */
772 : /************************************************************************/
773 :
774 2 : char **LANDataset::GetFileList()
775 :
776 : {
777 : // Main data file, etc.
778 2 : char **papszFileList = GDALPamDataset::GetFileList();
779 :
780 2 : if (!osSTAFilename.empty())
781 0 : papszFileList = CSLAddString(papszFileList, osSTAFilename);
782 :
783 2 : return papszFileList;
784 : }
785 :
786 : /************************************************************************/
787 : /* CheckForStatistics() */
788 : /************************************************************************/
789 :
790 24 : void LANDataset::CheckForStatistics()
791 :
792 : {
793 : /* -------------------------------------------------------------------- */
794 : /* Do we have a statistics file? */
795 : /* -------------------------------------------------------------------- */
796 24 : osSTAFilename = CPLResetExtension(GetDescription(), "sta");
797 :
798 24 : VSILFILE *fpSTA = VSIFOpenL(osSTAFilename, "r");
799 :
800 24 : if (fpSTA == nullptr && VSIIsCaseSensitiveFS(osSTAFilename))
801 : {
802 24 : osSTAFilename = CPLResetExtension(GetDescription(), "STA");
803 24 : fpSTA = VSIFOpenL(osSTAFilename, "r");
804 : }
805 :
806 24 : if (fpSTA == nullptr)
807 : {
808 24 : osSTAFilename = "";
809 24 : return;
810 : }
811 :
812 : /* -------------------------------------------------------------------- */
813 : /* Read it one band at a time. */
814 : /* -------------------------------------------------------------------- */
815 0 : GByte abyBandInfo[1152] = {'\0'};
816 :
817 0 : for (int iBand = 0; iBand < nBands; iBand++)
818 : {
819 0 : if (VSIFReadL(abyBandInfo, 1152, 1, fpSTA) != 1)
820 0 : break;
821 :
822 0 : const int nBandNumber = abyBandInfo[7];
823 0 : GDALRasterBand *poBand = GetRasterBand(nBandNumber);
824 0 : if (poBand == nullptr)
825 0 : break;
826 :
827 0 : GInt16 nMin = 0;
828 0 : GInt16 nMax = 0;
829 :
830 0 : if (poBand->GetRasterDataType() != GDT_Byte)
831 : {
832 0 : memcpy(&nMin, abyBandInfo + 28, 2);
833 0 : memcpy(&nMax, abyBandInfo + 30, 2);
834 0 : CPL_LSBPTR16(&nMin);
835 0 : CPL_LSBPTR16(&nMax);
836 : }
837 : else
838 : {
839 0 : nMin = abyBandInfo[9];
840 0 : nMax = abyBandInfo[8];
841 : }
842 :
843 0 : float fMean = 0.0;
844 0 : float fStdDev = 0.0;
845 0 : memcpy(&fMean, abyBandInfo + 12, 4);
846 0 : memcpy(&fStdDev, abyBandInfo + 24, 4);
847 0 : CPL_LSBPTR32(&fMean);
848 0 : CPL_LSBPTR32(&fStdDev);
849 :
850 0 : poBand->SetStatistics(nMin, nMax, fMean, fStdDev);
851 : }
852 :
853 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpSTA));
854 : }
855 :
856 : /************************************************************************/
857 : /* Create() */
858 : /************************************************************************/
859 :
860 60 : GDALDataset *LANDataset::Create(const char *pszFilename, int nXSize, int nYSize,
861 : int nBandsIn, GDALDataType eType,
862 : char ** /* papszOptions */)
863 : {
864 60 : if (eType != GDT_Byte && eType != GDT_Int16)
865 : {
866 33 : CPLError(CE_Failure, CPLE_AppDefined,
867 : "Attempt to create .GIS file with unsupported data type '%s'.",
868 : GDALGetDataTypeName(eType));
869 33 : return nullptr;
870 : }
871 :
872 : /* -------------------------------------------------------------------- */
873 : /* Try to create the file. */
874 : /* -------------------------------------------------------------------- */
875 27 : VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
876 27 : if (fp == nullptr)
877 : {
878 3 : CPLError(CE_Failure, CPLE_OpenFailed,
879 : "Attempt to create file `%s' failed.\n", pszFilename);
880 3 : return nullptr;
881 : }
882 :
883 : /* -------------------------------------------------------------------- */
884 : /* Write out the header. */
885 : /* -------------------------------------------------------------------- */
886 24 : unsigned char abyHeader[128] = {'\0'};
887 :
888 24 : memset(abyHeader, 0, sizeof(abyHeader));
889 :
890 24 : memcpy(abyHeader + 0, "HEAD74", 6);
891 :
892 : // Pixel type.
893 24 : GInt16 n16Val = 0;
894 24 : if (eType == GDT_Byte) // Do we want 4bit?
895 21 : n16Val = 0;
896 : else
897 3 : n16Val = 2;
898 24 : memcpy(abyHeader + 6, &n16Val, 2);
899 :
900 : // Number of Bands.
901 24 : n16Val = static_cast<GInt16>(nBandsIn);
902 24 : memcpy(abyHeader + 8, &n16Val, 2);
903 :
904 : // Unknown (6).
905 :
906 : // Width.
907 24 : GInt32 n32Val = nXSize;
908 24 : memcpy(abyHeader + 16, &n32Val, 4);
909 :
910 : // Height.
911 24 : n32Val = nYSize;
912 24 : memcpy(abyHeader + 20, &n32Val, 4);
913 :
914 : // X Start (4).
915 : // Y Start (4).
916 :
917 : // Unknown (56).
918 :
919 : // Coordinate System.
920 24 : n16Val = 0;
921 24 : memcpy(abyHeader + 88, &n16Val, 2);
922 :
923 : // Classes in coverage.
924 24 : n16Val = 0;
925 24 : memcpy(abyHeader + 90, &n16Val, 2);
926 :
927 : // Unknown (14).
928 :
929 : // Area Unit.
930 24 : n16Val = 0;
931 24 : memcpy(abyHeader + 106, &n16Val, 2);
932 :
933 : // Pixel Area.
934 24 : float f32Val = 0.0f;
935 24 : memcpy(abyHeader + 108, &f32Val, 4);
936 :
937 : // Upper Left X.
938 24 : f32Val = 0.5f;
939 24 : memcpy(abyHeader + 112, &f32Val, 4);
940 :
941 : // Upper Left Y
942 24 : f32Val = static_cast<float>(nYSize - 0.5);
943 24 : memcpy(abyHeader + 116, &f32Val, 4);
944 :
945 : // Width of pixel.
946 24 : f32Val = 1.0f;
947 24 : memcpy(abyHeader + 120, &f32Val, 4);
948 :
949 : // Height of pixel.
950 24 : f32Val = 1.0f;
951 24 : memcpy(abyHeader + 124, &f32Val, 4);
952 :
953 24 : CPL_IGNORE_RET_VAL(VSIFWriteL(abyHeader, sizeof(abyHeader), 1, fp));
954 :
955 : /* -------------------------------------------------------------------- */
956 : /* Extend the file to the target size. */
957 : /* -------------------------------------------------------------------- */
958 24 : vsi_l_offset nImageBytes = 0;
959 :
960 24 : if (eType != GDT_Byte)
961 3 : nImageBytes = nXSize * static_cast<vsi_l_offset>(nYSize) * 2;
962 : else
963 21 : nImageBytes = nXSize * static_cast<vsi_l_offset>(nYSize);
964 :
965 24 : memset(abyHeader, 0, sizeof(abyHeader));
966 :
967 820 : while (nImageBytes > 0)
968 : {
969 : const vsi_l_offset nWriteThisTime =
970 805 : std::min(static_cast<size_t>(nImageBytes), sizeof(abyHeader));
971 :
972 805 : if (VSIFWriteL(abyHeader, 1, static_cast<size_t>(nWriteThisTime), fp) !=
973 : nWriteThisTime)
974 : {
975 9 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
976 9 : CPLError(CE_Failure, CPLE_FileIO,
977 : "Failed to write whole Istar file.");
978 9 : return nullptr;
979 : }
980 796 : nImageBytes -= nWriteThisTime;
981 : }
982 :
983 15 : if (VSIFCloseL(fp) != 0)
984 : {
985 0 : CPLError(CE_Failure, CPLE_FileIO, "Failed to write whole Istar file.");
986 0 : return nullptr;
987 : }
988 :
989 15 : return static_cast<GDALDataset *>(GDALOpen(pszFilename, GA_Update));
990 : }
991 :
992 : /************************************************************************/
993 : /* GDALRegister_LAN() */
994 : /************************************************************************/
995 :
996 1595 : void GDALRegister_LAN()
997 :
998 : {
999 1595 : if (GDALGetDriverByName("LAN") != nullptr)
1000 302 : return;
1001 :
1002 1293 : GDALDriver *poDriver = new GDALDriver();
1003 :
1004 1293 : poDriver->SetDescription("LAN");
1005 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1006 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "Erdas .LAN/.GIS");
1007 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/lan.html");
1008 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1009 1293 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Byte Int16");
1010 :
1011 1293 : poDriver->pfnOpen = LANDataset::Open;
1012 1293 : poDriver->pfnCreate = LANDataset::Create;
1013 :
1014 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1015 : }
|
