Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: LCP Driver
4 : * Purpose: FARSITE v.4 Landscape file (.lcp) reader for GDAL
5 : * Author: Chris Toney
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2008, Chris Toney
9 : * Copyright (c) 2008-2011, Even Rouault <even dot rouault at spatialys.com>
10 : * Copyright (c) 2013, Kyle Shannon <kyle at pobox dot com>
11 : *
12 : * SPDX-License-Identifier: MIT
13 : ****************************************************************************/
14 :
15 : #include "cpl_port.h"
16 : #include "cpl_string.h"
17 : #include "gdal_frmts.h"
18 : #include "gdal_priv.h"
19 : #include "ogr_spatialref.h"
20 : #include "rawdataset.h"
21 :
22 : #include <algorithm>
23 : #include <limits>
24 :
25 : constexpr size_t LCP_HEADER_SIZE = 7316;
26 : constexpr int LCP_MAX_BANDS = 10;
27 : constexpr int LCP_MAX_PATH = 256;
28 : constexpr int LCP_MAX_DESC = 512;
29 : constexpr int LCP_MAX_CLASSES = 100;
30 :
31 : /************************************************************************/
32 : /* ==================================================================== */
33 : /* LCPDataset */
34 : /* ==================================================================== */
35 : /************************************************************************/
36 :
37 : class LCPDataset final : public RawDataset
38 : {
39 : VSILFILE *fpImage; // image data file.
40 : char pachHeader[LCP_HEADER_SIZE];
41 :
42 : CPLString osPrjFilename{};
43 : OGRSpatialReference m_oSRS{};
44 :
45 : static CPLErr ClassifyBandData(GDALRasterBand *poBand, GInt32 &nNumClasses,
46 : GInt32 *panClasses);
47 :
48 : CPL_DISALLOW_COPY_ASSIGN(LCPDataset)
49 :
50 : CPLErr Close() override;
51 :
52 : public:
53 : LCPDataset();
54 : ~LCPDataset() override;
55 :
56 : char **GetFileList(void) override;
57 :
58 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
59 :
60 : static int Identify(GDALOpenInfo *);
61 : static GDALDataset *Open(GDALOpenInfo *);
62 : static GDALDataset *CreateCopy(const char *pszFilename,
63 : GDALDataset *poSrcDS, int bStrict,
64 : char **papszOptions,
65 : GDALProgressFunc pfnProgress,
66 : void *pProgressData);
67 :
68 3 : const OGRSpatialReference *GetSpatialRef() const override
69 : {
70 3 : return &m_oSRS;
71 : }
72 : };
73 :
74 : /************************************************************************/
75 : /* LCPDataset() */
76 : /************************************************************************/
77 :
78 75 : LCPDataset::LCPDataset() : fpImage(nullptr)
79 : {
80 75 : memset(pachHeader, 0, sizeof(pachHeader));
81 75 : }
82 :
83 : /************************************************************************/
84 : /* ~LCPDataset() */
85 : /************************************************************************/
86 :
87 150 : LCPDataset::~LCPDataset()
88 :
89 : {
90 75 : LCPDataset::Close();
91 150 : }
92 :
93 : /************************************************************************/
94 : /* Close() */
95 : /************************************************************************/
96 :
97 132 : CPLErr LCPDataset::Close()
98 : {
99 132 : CPLErr eErr = CE_None;
100 132 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
101 : {
102 75 : if (LCPDataset::FlushCache(true) != CE_None)
103 0 : eErr = CE_Failure;
104 :
105 75 : if (fpImage)
106 : {
107 75 : if (VSIFCloseL(fpImage) != 0)
108 : {
109 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
110 0 : eErr = CE_Failure;
111 : }
112 : }
113 :
114 75 : if (GDALPamDataset::Close() != CE_None)
115 0 : eErr = CE_Failure;
116 : }
117 132 : return eErr;
118 : }
119 :
120 : /************************************************************************/
121 : /* GetGeoTransform() */
122 : /************************************************************************/
123 :
124 2 : CPLErr LCPDataset::GetGeoTransform(GDALGeoTransform >) const
125 : {
126 2 : double dfEast = 0.0;
127 2 : double dfWest = 0.0;
128 2 : double dfNorth = 0.0;
129 2 : double dfSouth = 0.0;
130 2 : double dfCellX = 0.0;
131 2 : double dfCellY = 0.0;
132 :
133 2 : memcpy(&dfEast, pachHeader + 4172, sizeof(double));
134 2 : memcpy(&dfWest, pachHeader + 4180, sizeof(double));
135 2 : memcpy(&dfNorth, pachHeader + 4188, sizeof(double));
136 2 : memcpy(&dfSouth, pachHeader + 4196, sizeof(double));
137 2 : memcpy(&dfCellX, pachHeader + 4208, sizeof(double));
138 2 : memcpy(&dfCellY, pachHeader + 4216, sizeof(double));
139 2 : CPL_LSBPTR64(&dfEast);
140 2 : CPL_LSBPTR64(&dfWest);
141 2 : CPL_LSBPTR64(&dfNorth);
142 2 : CPL_LSBPTR64(&dfSouth);
143 2 : CPL_LSBPTR64(&dfCellX);
144 2 : CPL_LSBPTR64(&dfCellY);
145 :
146 2 : gt[0] = dfWest;
147 2 : gt[3] = dfNorth;
148 2 : gt[1] = dfCellX;
149 2 : gt[2] = 0.0;
150 :
151 2 : gt[4] = 0.0;
152 2 : gt[5] = -1 * dfCellY;
153 :
154 2 : return CE_None;
155 : }
156 :
157 : /************************************************************************/
158 : /* Identify() */
159 : /************************************************************************/
160 :
161 59274 : int LCPDataset::Identify(GDALOpenInfo *poOpenInfo)
162 :
163 : {
164 : /* -------------------------------------------------------------------- */
165 : /* Verify that this is a FARSITE v.4 LCP file */
166 : /* -------------------------------------------------------------------- */
167 59274 : if (poOpenInfo->nHeaderBytes < 50)
168 55413 : return FALSE;
169 :
170 : /* check if first three fields have valid data */
171 3861 : if ((CPL_LSBSINT32PTR(poOpenInfo->pabyHeader) != 20 &&
172 3848 : CPL_LSBSINT32PTR(poOpenInfo->pabyHeader) != 21) ||
173 167 : (CPL_LSBSINT32PTR(poOpenInfo->pabyHeader + 4) != 20 &&
174 142 : CPL_LSBSINT32PTR(poOpenInfo->pabyHeader + 4) != 21) ||
175 167 : (CPL_LSBSINT32PTR(poOpenInfo->pabyHeader + 8) < -90 ||
176 168 : CPL_LSBSINT32PTR(poOpenInfo->pabyHeader + 8) > 90))
177 : {
178 3693 : return FALSE;
179 : }
180 :
181 : /* -------------------------------------------------------------------- */
182 : /* Check file extension */
183 : /* -------------------------------------------------------------------- */
184 : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
185 168 : const char *pszFileExtension = poOpenInfo->osExtension.c_str();
186 168 : if (!EQUAL(pszFileExtension, "lcp"))
187 : {
188 0 : return FALSE;
189 : }
190 : #endif
191 :
192 168 : return TRUE;
193 : }
194 :
195 : /************************************************************************/
196 : /* GetFileList() */
197 : /************************************************************************/
198 :
199 38 : char **LCPDataset::GetFileList()
200 :
201 : {
202 38 : char **papszFileList = GDALPamDataset::GetFileList();
203 :
204 38 : if (!m_oSRS.IsEmpty())
205 : {
206 1 : papszFileList = CSLAddString(papszFileList, osPrjFilename);
207 : }
208 :
209 38 : return papszFileList;
210 : }
211 :
212 : /************************************************************************/
213 : /* Open() */
214 : /************************************************************************/
215 :
216 75 : GDALDataset *LCPDataset::Open(GDALOpenInfo *poOpenInfo)
217 :
218 : {
219 : /* -------------------------------------------------------------------- */
220 : /* Verify that this is a FARSITE LCP file */
221 : /* -------------------------------------------------------------------- */
222 75 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
223 0 : return nullptr;
224 :
225 : /* -------------------------------------------------------------------- */
226 : /* Confirm the requested access is supported. */
227 : /* -------------------------------------------------------------------- */
228 75 : if (poOpenInfo->eAccess == GA_Update)
229 : {
230 0 : ReportUpdateNotSupportedByDriver("LCP");
231 0 : return nullptr;
232 : }
233 : /* -------------------------------------------------------------------- */
234 : /* Create a corresponding GDALDataset. */
235 : /* -------------------------------------------------------------------- */
236 150 : auto poDS = std::make_unique<LCPDataset>();
237 75 : std::swap(poDS->fpImage, poOpenInfo->fpL);
238 :
239 : /* -------------------------------------------------------------------- */
240 : /* Read the header and extract some information. */
241 : /* -------------------------------------------------------------------- */
242 150 : if (VSIFSeekL(poDS->fpImage, 0, SEEK_SET) < 0 ||
243 75 : VSIFReadL(poDS->pachHeader, 1, LCP_HEADER_SIZE, poDS->fpImage) !=
244 : LCP_HEADER_SIZE)
245 : {
246 0 : CPLError(CE_Failure, CPLE_FileIO, "File too short");
247 0 : return nullptr;
248 : }
249 :
250 75 : const int nWidth = CPL_LSBSINT32PTR(poDS->pachHeader + 4164);
251 75 : const int nHeight = CPL_LSBSINT32PTR(poDS->pachHeader + 4168);
252 :
253 75 : poDS->nRasterXSize = nWidth;
254 75 : poDS->nRasterYSize = nHeight;
255 :
256 75 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
257 : {
258 0 : return nullptr;
259 : }
260 :
261 : // Crown fuels = canopy height, canopy base height, canopy bulk density.
262 : // 21 = have them, 20 = do not have them
263 : const bool bHaveCrownFuels =
264 75 : CPL_TO_BOOL(CPL_LSBSINT32PTR(poDS->pachHeader + 0) - 20);
265 : // Ground fuels = duff loading, coarse woody.
266 : const bool bHaveGroundFuels =
267 75 : CPL_TO_BOOL(CPL_LSBSINT32PTR(poDS->pachHeader + 4) - 20);
268 :
269 75 : int nBands = 0;
270 75 : if (bHaveCrownFuels)
271 : {
272 69 : if (bHaveGroundFuels)
273 60 : nBands = 10;
274 : else
275 9 : nBands = 8;
276 : }
277 : else
278 : {
279 6 : if (bHaveGroundFuels)
280 3 : nBands = 7;
281 : else
282 3 : nBands = 5;
283 : }
284 :
285 : // Add dataset-level metadata.
286 :
287 75 : int nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 8);
288 75 : char szTemp[32] = {'\0'};
289 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
290 75 : poDS->SetMetadataItem("LATITUDE", szTemp);
291 :
292 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 4204);
293 75 : if (nTemp == 0)
294 74 : poDS->SetMetadataItem("LINEAR_UNIT", "Meters");
295 75 : if (nTemp == 1)
296 1 : poDS->SetMetadataItem("LINEAR_UNIT", "Feet");
297 :
298 75 : poDS->pachHeader[LCP_HEADER_SIZE - 1] = '\0';
299 75 : poDS->SetMetadataItem("DESCRIPTION", poDS->pachHeader + 6804);
300 :
301 : /* -------------------------------------------------------------------- */
302 : /* Create band information objects. */
303 : /* -------------------------------------------------------------------- */
304 :
305 75 : const int iPixelSize = nBands * 2;
306 :
307 75 : if (nWidth > INT_MAX / iPixelSize)
308 : {
309 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow occurred");
310 0 : return nullptr;
311 : }
312 :
313 783 : for (int iBand = 1; iBand <= nBands; iBand++)
314 : {
315 : auto poBand =
316 1416 : RawRasterBand::Create(poDS.get(), iBand, poDS->fpImage,
317 708 : LCP_HEADER_SIZE + ((iBand - 1) * 2),
318 : iPixelSize, iPixelSize * nWidth, GDT_Int16,
319 : RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN,
320 708 : RawRasterBand::OwnFP::NO);
321 708 : if (!poBand)
322 0 : return nullptr;
323 :
324 708 : switch (iBand)
325 : {
326 75 : case 1:
327 75 : poBand->SetDescription("Elevation");
328 :
329 75 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4224);
330 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
331 75 : poBand->SetMetadataItem("ELEVATION_UNIT", szTemp);
332 :
333 75 : if (nTemp == 0)
334 74 : poBand->SetMetadataItem("ELEVATION_UNIT_NAME", "Meters");
335 75 : if (nTemp == 1)
336 1 : poBand->SetMetadataItem("ELEVATION_UNIT_NAME", "Feet");
337 :
338 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 44);
339 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
340 75 : poBand->SetMetadataItem("ELEVATION_MIN", szTemp);
341 :
342 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 48);
343 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
344 75 : poBand->SetMetadataItem("ELEVATION_MAX", szTemp);
345 :
346 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 52);
347 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
348 75 : poBand->SetMetadataItem("ELEVATION_NUM_CLASSES", szTemp);
349 :
350 75 : *(poDS->pachHeader + 4244 + 255) = '\0';
351 75 : poBand->SetMetadataItem("ELEVATION_FILE",
352 75 : poDS->pachHeader + 4244);
353 :
354 75 : break;
355 :
356 75 : case 2:
357 75 : poBand->SetDescription("Slope");
358 :
359 75 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4226);
360 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
361 75 : poBand->SetMetadataItem("SLOPE_UNIT", szTemp);
362 :
363 75 : if (nTemp == 0)
364 74 : poBand->SetMetadataItem("SLOPE_UNIT_NAME", "Degrees");
365 75 : if (nTemp == 1)
366 1 : poBand->SetMetadataItem("SLOPE_UNIT_NAME", "Percent");
367 :
368 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 456);
369 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
370 75 : poBand->SetMetadataItem("SLOPE_MIN", szTemp);
371 :
372 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 460);
373 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
374 75 : poBand->SetMetadataItem("SLOPE_MAX", szTemp);
375 :
376 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 464);
377 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
378 75 : poBand->SetMetadataItem("SLOPE_NUM_CLASSES", szTemp);
379 :
380 75 : *(poDS->pachHeader + 4500 + 255) = '\0';
381 75 : poBand->SetMetadataItem("SLOPE_FILE", poDS->pachHeader + 4500);
382 :
383 75 : break;
384 :
385 75 : case 3:
386 75 : poBand->SetDescription("Aspect");
387 :
388 75 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4228);
389 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
390 75 : poBand->SetMetadataItem("ASPECT_UNIT", szTemp);
391 :
392 75 : if (nTemp == 0)
393 3 : poBand->SetMetadataItem("ASPECT_UNIT_NAME",
394 3 : "Grass categories");
395 75 : if (nTemp == 1)
396 1 : poBand->SetMetadataItem("ASPECT_UNIT_NAME",
397 1 : "Grass degrees");
398 75 : if (nTemp == 2)
399 71 : poBand->SetMetadataItem("ASPECT_UNIT_NAME",
400 71 : "Azimuth degrees");
401 :
402 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 868);
403 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
404 75 : poBand->SetMetadataItem("ASPECT_MIN", szTemp);
405 :
406 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 872);
407 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
408 75 : poBand->SetMetadataItem("ASPECT_MAX", szTemp);
409 :
410 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 876);
411 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
412 75 : poBand->SetMetadataItem("ASPECT_NUM_CLASSES", szTemp);
413 :
414 75 : *(poDS->pachHeader + 4756 + 255) = '\0';
415 75 : poBand->SetMetadataItem("ASPECT_FILE", poDS->pachHeader + 4756);
416 :
417 75 : break;
418 :
419 75 : case 4:
420 : {
421 75 : poBand->SetDescription("Fuel models");
422 :
423 75 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4230);
424 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
425 75 : poBand->SetMetadataItem("FUEL_MODEL_OPTION", szTemp);
426 :
427 75 : if (nTemp == 0)
428 75 : poBand->SetMetadataItem(
429 : "FUEL_MODEL_OPTION_DESC",
430 75 : "no custom models AND no conversion file needed");
431 75 : if (nTemp == 1)
432 0 : poBand->SetMetadataItem(
433 : "FUEL_MODEL_OPTION_DESC",
434 0 : "custom models BUT no conversion file needed");
435 75 : if (nTemp == 2)
436 0 : poBand->SetMetadataItem(
437 : "FUEL_MODEL_OPTION_DESC",
438 0 : "no custom models BUT conversion file needed");
439 75 : if (nTemp == 3)
440 0 : poBand->SetMetadataItem(
441 : "FUEL_MODEL_OPTION_DESC",
442 0 : "custom models AND conversion file needed");
443 :
444 75 : const int nMinFM = CPL_LSBSINT32PTR(poDS->pachHeader + 1280);
445 75 : snprintf(szTemp, sizeof(szTemp), "%d", nMinFM);
446 75 : poBand->SetMetadataItem("FUEL_MODEL_MIN", szTemp);
447 :
448 75 : const int nMaxFM = CPL_LSBSINT32PTR(poDS->pachHeader + 1284);
449 75 : snprintf(szTemp, sizeof(szTemp), "%d", nMaxFM);
450 75 : poBand->SetMetadataItem("FUEL_MODEL_MAX", szTemp);
451 :
452 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 1288);
453 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
454 75 : poBand->SetMetadataItem("FUEL_MODEL_NUM_CLASSES", szTemp);
455 :
456 150 : std::string osValues;
457 75 : if (nTemp > 0 && nTemp <= 100)
458 : {
459 268 : for (int i = 0; i <= nTemp; i++)
460 : {
461 193 : const int nTemp2 = CPL_LSBSINT32PTR(poDS->pachHeader +
462 : (1292 + (i * 4)));
463 193 : if (nTemp2 >= nMinFM && nTemp2 <= nMaxFM)
464 : {
465 100 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp2);
466 100 : if (!osValues.empty())
467 27 : osValues += ',';
468 100 : osValues += szTemp;
469 : }
470 : }
471 : }
472 75 : poBand->SetMetadataItem("FUEL_MODEL_VALUES", osValues.c_str());
473 :
474 75 : *(poDS->pachHeader + 5012 + 255) = '\0';
475 75 : poBand->SetMetadataItem("FUEL_MODEL_FILE",
476 75 : poDS->pachHeader + 5012);
477 :
478 75 : break;
479 : }
480 75 : case 5:
481 75 : poBand->SetDescription("Canopy cover");
482 :
483 75 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4232);
484 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
485 75 : poBand->SetMetadataItem("CANOPY_COV_UNIT", szTemp);
486 :
487 75 : if (nTemp == 0)
488 4 : poBand->SetMetadataItem("CANOPY_COV_UNIT_NAME",
489 4 : "Categories (0-4)");
490 75 : if (nTemp == 1)
491 71 : poBand->SetMetadataItem("CANOPY_COV_UNIT_NAME", "Percent");
492 :
493 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 1692);
494 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
495 75 : poBand->SetMetadataItem("CANOPY_COV_MIN", szTemp);
496 :
497 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 1696);
498 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
499 75 : poBand->SetMetadataItem("CANOPY_COV_MAX", szTemp);
500 :
501 75 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 1700);
502 75 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
503 75 : poBand->SetMetadataItem("CANOPY_COV_NUM_CLASSES", szTemp);
504 :
505 75 : *(poDS->pachHeader + 5268 + 255) = '\0';
506 75 : poBand->SetMetadataItem("CANOPY_COV_FILE",
507 75 : poDS->pachHeader + 5268);
508 :
509 75 : break;
510 :
511 72 : case 6:
512 72 : if (bHaveCrownFuels)
513 : {
514 69 : poBand->SetDescription("Canopy height");
515 :
516 69 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4234);
517 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
518 69 : poBand->SetMetadataItem("CANOPY_HT_UNIT", szTemp);
519 :
520 69 : if (nTemp == 1)
521 3 : poBand->SetMetadataItem("CANOPY_HT_UNIT_NAME",
522 3 : "Meters");
523 69 : if (nTemp == 2)
524 3 : poBand->SetMetadataItem("CANOPY_HT_UNIT_NAME", "Feet");
525 69 : if (nTemp == 3)
526 62 : poBand->SetMetadataItem("CANOPY_HT_UNIT_NAME",
527 62 : "Meters x 10");
528 69 : if (nTemp == 4)
529 1 : poBand->SetMetadataItem("CANOPY_HT_UNIT_NAME",
530 1 : "Feet x 10");
531 :
532 69 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 2104);
533 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
534 69 : poBand->SetMetadataItem("CANOPY_HT_MIN", szTemp);
535 :
536 69 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 2108);
537 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
538 69 : poBand->SetMetadataItem("CANOPY_HT_MAX", szTemp);
539 :
540 69 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 2112);
541 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
542 69 : poBand->SetMetadataItem("CANOPY_HT_NUM_CLASSES", szTemp);
543 :
544 69 : *(poDS->pachHeader + 5524 + 255) = '\0';
545 69 : poBand->SetMetadataItem("CANOPY_HT_FILE",
546 69 : poDS->pachHeader + 5524);
547 : }
548 : else
549 : {
550 3 : poBand->SetDescription("Duff");
551 :
552 3 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4240);
553 3 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
554 3 : poBand->SetMetadataItem("DUFF_UNIT", szTemp);
555 :
556 3 : if (nTemp == 1)
557 3 : poBand->SetMetadataItem("DUFF_UNIT_NAME", "Mg/ha");
558 3 : if (nTemp == 2)
559 0 : poBand->SetMetadataItem("DUFF_UNIT_NAME", "t/ac");
560 :
561 3 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3340);
562 3 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
563 3 : poBand->SetMetadataItem("DUFF_MIN", szTemp);
564 :
565 3 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3344);
566 3 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
567 3 : poBand->SetMetadataItem("DUFF_MAX", szTemp);
568 :
569 3 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3348);
570 3 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
571 3 : poBand->SetMetadataItem("DUFF_NUM_CLASSES", szTemp);
572 :
573 3 : *(poDS->pachHeader + 6292 + 255) = '\0';
574 3 : poBand->SetMetadataItem("DUFF_FILE",
575 3 : poDS->pachHeader + 6292);
576 : }
577 72 : break;
578 :
579 72 : case 7:
580 72 : if (bHaveCrownFuels)
581 : {
582 69 : poBand->SetDescription("Canopy base height");
583 :
584 69 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4236);
585 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
586 69 : poBand->SetMetadataItem("CBH_UNIT", szTemp);
587 :
588 69 : if (nTemp == 1)
589 3 : poBand->SetMetadataItem("CBH_UNIT_NAME", "Meters");
590 69 : if (nTemp == 2)
591 3 : poBand->SetMetadataItem("CBH_UNIT_NAME", "Feet");
592 69 : if (nTemp == 3)
593 62 : poBand->SetMetadataItem("CBH_UNIT_NAME", "Meters x 10");
594 69 : if (nTemp == 4)
595 1 : poBand->SetMetadataItem("CBH_UNIT_NAME", "Feet x 10");
596 :
597 69 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 2516);
598 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
599 69 : poBand->SetMetadataItem("CBH_MIN", szTemp);
600 :
601 69 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 2520);
602 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
603 69 : poBand->SetMetadataItem("CBH_MAX", szTemp);
604 :
605 69 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 2524);
606 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
607 69 : poBand->SetMetadataItem("CBH_NUM_CLASSES", szTemp);
608 :
609 69 : *(poDS->pachHeader + 5780 + 255) = '\0';
610 69 : poBand->SetMetadataItem("CBH_FILE",
611 69 : poDS->pachHeader + 5780);
612 : }
613 : else
614 : {
615 3 : poBand->SetDescription("Coarse woody debris");
616 :
617 3 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4242);
618 3 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
619 3 : poBand->SetMetadataItem("CWD_OPTION", szTemp);
620 :
621 : // if ( nTemp == 1 )
622 : // poBand->SetMetadataItem( "CWD_UNIT_DESC", "?" );
623 : // if ( nTemp == 2 )
624 : // poBand->SetMetadataItem( "CWD_UNIT_DESC", "?" );
625 :
626 3 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3752);
627 3 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
628 3 : poBand->SetMetadataItem("CWD_MIN", szTemp);
629 :
630 3 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3756);
631 3 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
632 3 : poBand->SetMetadataItem("CWD_MAX", szTemp);
633 :
634 3 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3760);
635 3 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
636 3 : poBand->SetMetadataItem("CWD_NUM_CLASSES", szTemp);
637 :
638 3 : *(poDS->pachHeader + 6548 + 255) = '\0';
639 3 : poBand->SetMetadataItem("CWD_FILE",
640 3 : poDS->pachHeader + 6548);
641 : }
642 72 : break;
643 :
644 69 : case 8:
645 69 : poBand->SetDescription("Canopy bulk density");
646 :
647 69 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4238);
648 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
649 69 : poBand->SetMetadataItem("CBD_UNIT", szTemp);
650 :
651 69 : if (nTemp == 1)
652 3 : poBand->SetMetadataItem("CBD_UNIT_NAME", "kg/m^3");
653 69 : if (nTemp == 2)
654 3 : poBand->SetMetadataItem("CBD_UNIT_NAME", "lb/ft^3");
655 69 : if (nTemp == 3)
656 62 : poBand->SetMetadataItem("CBD_UNIT_NAME", "kg/m^3 x 100");
657 69 : if (nTemp == 4)
658 1 : poBand->SetMetadataItem("CBD_UNIT_NAME", "lb/ft^3 x 1000");
659 :
660 69 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 2928);
661 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
662 69 : poBand->SetMetadataItem("CBD_MIN", szTemp);
663 :
664 69 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 2932);
665 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
666 69 : poBand->SetMetadataItem("CBD_MAX", szTemp);
667 :
668 69 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 2936);
669 69 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
670 69 : poBand->SetMetadataItem("CBD_NUM_CLASSES", szTemp);
671 :
672 69 : *(poDS->pachHeader + 6036 + 255) = '\0';
673 69 : poBand->SetMetadataItem("CBD_FILE", poDS->pachHeader + 6036);
674 :
675 69 : break;
676 :
677 60 : case 9:
678 60 : poBand->SetDescription("Duff");
679 :
680 60 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4240);
681 60 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
682 60 : poBand->SetMetadataItem("DUFF_UNIT", szTemp);
683 :
684 60 : if (nTemp == 1)
685 59 : poBand->SetMetadataItem("DUFF_UNIT_NAME", "Mg/ha");
686 60 : if (nTemp == 2)
687 1 : poBand->SetMetadataItem("DUFF_UNIT_NAME", "t/ac");
688 :
689 60 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3340);
690 60 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
691 60 : poBand->SetMetadataItem("DUFF_MIN", szTemp);
692 :
693 60 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3344);
694 60 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
695 60 : poBand->SetMetadataItem("DUFF_MAX", szTemp);
696 :
697 60 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3348);
698 60 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
699 60 : poBand->SetMetadataItem("DUFF_NUM_CLASSES", szTemp);
700 :
701 60 : *(poDS->pachHeader + 6292 + 255) = '\0';
702 60 : poBand->SetMetadataItem("DUFF_FILE", poDS->pachHeader + 6292);
703 :
704 60 : break;
705 :
706 60 : case 10:
707 60 : poBand->SetDescription("Coarse woody debris");
708 :
709 60 : nTemp = CPL_LSBUINT16PTR(poDS->pachHeader + 4242);
710 60 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
711 60 : poBand->SetMetadataItem("CWD_OPTION", szTemp);
712 :
713 : // if ( nTemp == 1 )
714 : // poBand->SetMetadataItem( "CWD_UNIT_DESC", "?" );
715 : // if ( nTemp == 2 )
716 : // poBand->SetMetadataItem( "CWD_UNIT_DESC", "?" );
717 :
718 60 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3752);
719 60 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
720 60 : poBand->SetMetadataItem("CWD_MIN", szTemp);
721 :
722 60 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3756);
723 60 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
724 60 : poBand->SetMetadataItem("CWD_MAX", szTemp);
725 :
726 60 : nTemp = CPL_LSBSINT32PTR(poDS->pachHeader + 3760);
727 60 : snprintf(szTemp, sizeof(szTemp), "%d", nTemp);
728 60 : poBand->SetMetadataItem("CWD_NUM_CLASSES", szTemp);
729 :
730 60 : *(poDS->pachHeader + 6548 + 255) = '\0';
731 60 : poBand->SetMetadataItem("CWD_FILE", poDS->pachHeader + 6548);
732 :
733 60 : break;
734 : }
735 :
736 708 : poDS->SetBand(iBand, std::move(poBand));
737 : }
738 :
739 : /* -------------------------------------------------------------------- */
740 : /* Try to read projection file. */
741 : /* -------------------------------------------------------------------- */
742 : char *const pszDirname =
743 75 : CPLStrdup(CPLGetPathSafe(poOpenInfo->pszFilename).c_str());
744 : char *const pszBasename =
745 75 : CPLStrdup(CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str());
746 :
747 75 : poDS->osPrjFilename = CPLFormFilenameSafe(pszDirname, pszBasename, "prj");
748 : VSIStatBufL sStatBuf;
749 75 : int nRet = VSIStatL(poDS->osPrjFilename, &sStatBuf);
750 :
751 75 : if (nRet != 0 && VSIIsCaseSensitiveFS(poDS->osPrjFilename))
752 : {
753 72 : poDS->osPrjFilename =
754 144 : CPLFormFilenameSafe(pszDirname, pszBasename, "PRJ");
755 72 : nRet = VSIStatL(poDS->osPrjFilename, &sStatBuf);
756 : }
757 :
758 75 : if (nRet == 0)
759 : {
760 3 : char **papszPrj = CSLLoad(poDS->osPrjFilename);
761 :
762 3 : CPLDebug("LCP", "Loaded SRS from %s", poDS->osPrjFilename.c_str());
763 :
764 3 : poDS->m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
765 3 : if (poDS->m_oSRS.importFromESRI(papszPrj) != OGRERR_NONE)
766 : {
767 0 : poDS->m_oSRS.Clear();
768 : }
769 :
770 3 : CSLDestroy(papszPrj);
771 : }
772 :
773 75 : CPLFree(pszDirname);
774 75 : CPLFree(pszBasename);
775 :
776 : /* -------------------------------------------------------------------- */
777 : /* Initialize any PAM information. */
778 : /* -------------------------------------------------------------------- */
779 75 : poDS->SetDescription(poOpenInfo->pszFilename);
780 75 : poDS->TryLoadXML();
781 :
782 : /* -------------------------------------------------------------------- */
783 : /* Check for external overviews. */
784 : /* -------------------------------------------------------------------- */
785 150 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename,
786 75 : poOpenInfo->GetSiblingFiles());
787 :
788 75 : return poDS.release();
789 : }
790 :
791 : /************************************************************************/
792 : /* ClassifyBandData() */
793 : /* Classify a band and store 99 or fewer unique values. If there are */
794 : /* more than 99 unique values, then set pnNumClasses to -1 as a flag */
795 : /* that represents this. These are legacy values in the header, and */
796 : /* while we should never deprecate them, we could possibly not */
797 : /* calculate them by default. */
798 : /************************************************************************/
799 :
800 320 : CPLErr LCPDataset::ClassifyBandData(GDALRasterBand *poBand, GInt32 &nNumClasses,
801 : GInt32 *panClasses)
802 : {
803 320 : CPLAssert(poBand);
804 320 : CPLAssert(panClasses);
805 :
806 320 : const int nXSize = poBand->GetXSize();
807 320 : const int nYSize = poBand->GetYSize();
808 :
809 : GInt16 *panValues =
810 320 : static_cast<GInt16 *>(CPLMalloc(sizeof(GInt16) * nXSize));
811 320 : constexpr int MIN_VAL = std::numeric_limits<GInt16>::min();
812 320 : constexpr int MAX_VAL = std::numeric_limits<GInt16>::max();
813 320 : constexpr int RANGE_VAL = MAX_VAL - MIN_VAL + 1;
814 320 : GByte *pabyFlags = static_cast<GByte *>(CPLCalloc(1, RANGE_VAL));
815 :
816 : /* so that values in [-32768,32767] are mapped to [0,65535] in pabyFlags */
817 320 : constexpr int OFFSET = -MIN_VAL;
818 :
819 320 : int nFound = 0;
820 320 : bool bTooMany = false;
821 320 : CPLErr eErr = CE_None;
822 3770 : for (int iLine = 0; iLine < nYSize; iLine++)
823 : {
824 3451 : eErr = poBand->RasterIO(GF_Read, 0, iLine, nXSize, 1, panValues, nXSize,
825 : 1, GDT_Int16, 0, 0, nullptr);
826 3451 : if (eErr != CE_None)
827 0 : break;
828 37531 : for (int iPixel = 0; iPixel < nXSize; iPixel++)
829 : {
830 34081 : if (panValues[iPixel] == -9999)
831 : {
832 0 : continue;
833 : }
834 34081 : if (nFound == LCP_MAX_CLASSES)
835 : {
836 1 : CPLDebug("LCP",
837 : "Found more that %d unique values in "
838 : "band %d. Not 'classifying' the data.",
839 : LCP_MAX_CLASSES - 1, poBand->GetBand());
840 1 : nFound = -1;
841 1 : bTooMany = true;
842 1 : break;
843 : }
844 34080 : if (pabyFlags[panValues[iPixel] + OFFSET] == 0)
845 : {
846 573 : pabyFlags[panValues[iPixel] + OFFSET] = 1;
847 573 : nFound++;
848 : }
849 : }
850 3451 : if (bTooMany)
851 1 : break;
852 : }
853 320 : if (!bTooMany)
854 : {
855 : // The classes are always padded with a leading 0. This was for aligning
856 : // offsets, or making it a 1-based array instead of 0-based.
857 319 : panClasses[0] = 0;
858 20906300 : for (int j = 0, nIndex = 1; j < RANGE_VAL; j++)
859 : {
860 20906000 : if (pabyFlags[j] == 1)
861 : {
862 473 : panClasses[nIndex++] = j;
863 : }
864 : }
865 : }
866 320 : nNumClasses = nFound;
867 320 : CPLFree(pabyFlags);
868 320 : CPLFree(panValues);
869 :
870 320 : return eErr;
871 : }
872 :
873 : /************************************************************************/
874 : /* CreateCopy() */
875 : /************************************************************************/
876 :
877 68 : GDALDataset *LCPDataset::CreateCopy(const char *pszFilename,
878 : GDALDataset *poSrcDS, int bStrict,
879 : char **papszOptions,
880 : GDALProgressFunc pfnProgress,
881 : void *pProgressData)
882 :
883 : {
884 : /* -------------------------------------------------------------------- */
885 : /* Verify input options. */
886 : /* -------------------------------------------------------------------- */
887 68 : const int nBands = poSrcDS->GetRasterCount();
888 68 : if (nBands != 5 && nBands != 7 && nBands != 8 && nBands != 10)
889 : {
890 25 : CPLError(CE_Failure, CPLE_NotSupported,
891 : "LCP driver doesn't support %d bands. Must be 5, 7, 8 "
892 : "or 10 bands.",
893 : nBands);
894 25 : return nullptr;
895 : }
896 :
897 43 : GDALDataType eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
898 43 : if (eType != GDT_Int16 && bStrict)
899 : {
900 1 : CPLError(CE_Failure, CPLE_AppDefined,
901 : "LCP only supports 16-bit signed integer data types.");
902 1 : return nullptr;
903 : }
904 42 : else if (eType != GDT_Int16)
905 : {
906 0 : CPLError(CE_Warning, CPLE_AppDefined,
907 : "Setting data type to 16-bit integer.");
908 : }
909 :
910 : /* -------------------------------------------------------------------- */
911 : /* What schema do we have (ground/crown fuels) */
912 : /* -------------------------------------------------------------------- */
913 :
914 42 : const bool bHaveCrownFuels = nBands == 8 || nBands == 10;
915 42 : const bool bHaveGroundFuels = nBands == 7 || nBands == 10;
916 :
917 : // Since units are 'configurable', we should check for user
918 : // defined units. This is a bit cumbersome, but the user should
919 : // be allowed to specify none to get default units/options. Use
920 : // default units every chance we get.
921 42 : GInt16 panMetadata[LCP_MAX_BANDS] = {
922 : 0, // 0 ELEVATION_UNIT
923 : 0, // 1 SLOPE_UNIT
924 : 2, // 2 ASPECT_UNIT
925 : 0, // 3 FUEL_MODEL_OPTION
926 : 1, // 4 CANOPY_COV_UNIT
927 : 3, // 5 CANOPY_HT_UNIT
928 : 3, // 6 CBH_UNIT
929 : 3, // 7 CBD_UNIT
930 : 1, // 8 DUFF_UNIT
931 : 0, // 9 CWD_OPTION
932 : };
933 :
934 : // Check the units/options for user overrides.
935 : const char *pszTemp =
936 42 : CSLFetchNameValueDef(papszOptions, "ELEVATION_UNIT", "METERS");
937 42 : if (STARTS_WITH_CI(pszTemp, "METER"))
938 : {
939 40 : panMetadata[0] = 0;
940 : }
941 2 : else if (EQUAL(pszTemp, "FEET") || EQUAL(pszTemp, "FOOT"))
942 : {
943 1 : panMetadata[0] = 1;
944 : }
945 : else
946 : {
947 1 : CPLError(CE_Failure, CPLE_AppDefined,
948 : "Invalid value (%s) for ELEVATION_UNIT.", pszTemp);
949 1 : return nullptr;
950 : }
951 :
952 41 : pszTemp = CSLFetchNameValueDef(papszOptions, "SLOPE_UNIT", "DEGREES");
953 41 : if (EQUAL(pszTemp, "DEGREES"))
954 : {
955 39 : panMetadata[1] = 0;
956 : }
957 2 : else if (EQUAL(pszTemp, "PERCENT"))
958 : {
959 1 : panMetadata[1] = 1;
960 : }
961 : else
962 : {
963 1 : CPLError(CE_Failure, CPLE_AppDefined,
964 : "Invalid value (%s) for SLOPE_UNIT.", pszTemp);
965 1 : return nullptr;
966 : }
967 :
968 : pszTemp =
969 40 : CSLFetchNameValueDef(papszOptions, "ASPECT_UNIT", "AZIMUTH_DEGREES");
970 40 : if (EQUAL(pszTemp, "GRASS_CATEGORIES"))
971 : {
972 1 : panMetadata[2] = 0;
973 : }
974 39 : else if (EQUAL(pszTemp, "GRASS_DEGREES"))
975 : {
976 1 : panMetadata[2] = 1;
977 : }
978 38 : else if (EQUAL(pszTemp, "AZIMUTH_DEGREES"))
979 : {
980 37 : panMetadata[2] = 2;
981 : }
982 : else
983 : {
984 1 : CPLError(CE_Failure, CPLE_AppDefined,
985 : "Invalid value (%s) for ASPECT_UNIT.", pszTemp);
986 1 : return nullptr;
987 : }
988 :
989 39 : pszTemp = CSLFetchNameValueDef(papszOptions, "FUEL_MODEL_OPTION",
990 : "NO_CUSTOM_AND_NO_FILE");
991 39 : if (EQUAL(pszTemp, "NO_CUSTOM_AND_NO_FILE"))
992 : {
993 38 : panMetadata[3] = 0;
994 : }
995 1 : else if (EQUAL(pszTemp, "CUSTOM_AND_NO_FILE"))
996 : {
997 0 : panMetadata[3] = 1;
998 : }
999 1 : else if (EQUAL(pszTemp, "NO_CUSTOM_AND_FILE"))
1000 : {
1001 0 : panMetadata[3] = 2;
1002 : }
1003 1 : else if (EQUAL(pszTemp, "CUSTOM_AND_FILE"))
1004 : {
1005 0 : panMetadata[3] = 3;
1006 : }
1007 : else
1008 : {
1009 1 : CPLError(CE_Failure, CPLE_AppDefined,
1010 : "Invalid value (%s) for FUEL_MODEL_OPTION.", pszTemp);
1011 1 : return nullptr;
1012 : }
1013 :
1014 38 : pszTemp = CSLFetchNameValueDef(papszOptions, "CANOPY_COV_UNIT", "PERCENT");
1015 38 : if (EQUAL(pszTemp, "CATEGORIES"))
1016 : {
1017 1 : panMetadata[4] = 0;
1018 : }
1019 37 : else if (EQUAL(pszTemp, "PERCENT"))
1020 : {
1021 36 : panMetadata[4] = 1;
1022 : }
1023 : else
1024 : {
1025 1 : CPLError(CE_Failure, CPLE_AppDefined,
1026 : "Invalid value (%s) for CANOPY_COV_UNIT.", pszTemp);
1027 1 : return nullptr;
1028 : }
1029 :
1030 37 : if (bHaveCrownFuels)
1031 : {
1032 : pszTemp =
1033 35 : CSLFetchNameValueDef(papszOptions, "CANOPY_HT_UNIT", "METERS_X_10");
1034 35 : if (EQUAL(pszTemp, "METERS") || EQUAL(pszTemp, "METER"))
1035 : {
1036 1 : panMetadata[5] = 1;
1037 : }
1038 34 : else if (EQUAL(pszTemp, "FEET") || EQUAL(pszTemp, "FOOT"))
1039 : {
1040 1 : panMetadata[5] = 2;
1041 : }
1042 33 : else if (EQUAL(pszTemp, "METERS_X_10") || EQUAL(pszTemp, "METER_X_10"))
1043 : {
1044 31 : panMetadata[5] = 3;
1045 : }
1046 2 : else if (EQUAL(pszTemp, "FEET_X_10") || EQUAL(pszTemp, "FOOT_X_10"))
1047 : {
1048 1 : panMetadata[5] = 4;
1049 : }
1050 : else
1051 : {
1052 1 : CPLError(CE_Failure, CPLE_AppDefined,
1053 : "Invalid value (%s) for CANOPY_HT_UNIT.", pszTemp);
1054 1 : return nullptr;
1055 : }
1056 :
1057 34 : pszTemp = CSLFetchNameValueDef(papszOptions, "CBH_UNIT", "METERS_X_10");
1058 34 : if (EQUAL(pszTemp, "METERS") || EQUAL(pszTemp, "METER"))
1059 : {
1060 1 : panMetadata[6] = 1;
1061 : }
1062 33 : else if (EQUAL(pszTemp, "FEET") || EQUAL(pszTemp, "FOOT"))
1063 : {
1064 1 : panMetadata[6] = 2;
1065 : }
1066 32 : else if (EQUAL(pszTemp, "METERS_X_10") || EQUAL(pszTemp, "METER_X_10"))
1067 : {
1068 30 : panMetadata[6] = 3;
1069 : }
1070 2 : else if (EQUAL(pszTemp, "FEET_X_10") || EQUAL(pszTemp, "FOOT_X_10"))
1071 : {
1072 1 : panMetadata[6] = 4;
1073 : }
1074 : else
1075 : {
1076 1 : CPLError(CE_Failure, CPLE_AppDefined,
1077 : "Invalid value (%s) for CBH_UNIT.", pszTemp);
1078 1 : return nullptr;
1079 : }
1080 :
1081 33 : pszTemp = CSLFetchNameValueDef(papszOptions, "CBD_UNIT",
1082 : "KG_PER_CUBIC_METER_X_100");
1083 33 : if (EQUAL(pszTemp, "KG_PER_CUBIC_METER"))
1084 : {
1085 1 : panMetadata[7] = 1;
1086 : }
1087 32 : else if (EQUAL(pszTemp, "POUND_PER_CUBIC_FOOT"))
1088 : {
1089 1 : panMetadata[7] = 2;
1090 : }
1091 31 : else if (EQUAL(pszTemp, "KG_PER_CUBIC_METER_X_100"))
1092 : {
1093 29 : panMetadata[7] = 3;
1094 : }
1095 2 : else if (EQUAL(pszTemp, "POUND_PER_CUBIC_FOOT_X_1000"))
1096 : {
1097 1 : panMetadata[7] = 4;
1098 : }
1099 : else
1100 : {
1101 1 : CPLError(CE_Failure, CPLE_AppDefined,
1102 : "Invalid value (%s) for CBD_UNIT.", pszTemp);
1103 1 : return nullptr;
1104 : }
1105 : }
1106 :
1107 34 : if (bHaveGroundFuels)
1108 : {
1109 32 : pszTemp = CSLFetchNameValueDef(papszOptions, "DUFF_UNIT",
1110 : "MG_PER_HECTARE_X_10");
1111 32 : if (EQUAL(pszTemp, "MG_PER_HECTARE_X_10"))
1112 : {
1113 30 : panMetadata[8] = 1;
1114 : }
1115 2 : else if (EQUAL(pszTemp, "TONS_PER_ACRE_X_10"))
1116 : {
1117 1 : panMetadata[8] = 2;
1118 : }
1119 : else
1120 : {
1121 1 : CPLError(CE_Failure, CPLE_AppDefined,
1122 : "Invalid value (%s) for DUFF_UNIT.", pszTemp);
1123 1 : return nullptr;
1124 : }
1125 :
1126 31 : panMetadata[9] = 1;
1127 : }
1128 :
1129 : // Calculate the stats for each band. The binary file carries along
1130 : // these metadata for display purposes(?).
1131 33 : bool bCalculateStats = CPLFetchBool(papszOptions, "CALCULATE_STATS", true);
1132 :
1133 : const bool bClassifyData =
1134 33 : CPLFetchBool(papszOptions, "CLASSIFY_DATA", true);
1135 :
1136 : // We should have stats if we classify, we'll get them anyway.
1137 33 : if (bClassifyData && !bCalculateStats)
1138 : {
1139 0 : CPLError(CE_Warning, CPLE_AppDefined,
1140 : "Ignoring request to not calculate statistics, "
1141 : "because CLASSIFY_DATA was set to ON");
1142 0 : bCalculateStats = true;
1143 : }
1144 :
1145 33 : pszTemp = CSLFetchNameValueDef(papszOptions, "LINEAR_UNIT", "SET_FROM_SRS");
1146 33 : int nLinearUnits = 0;
1147 33 : bool bSetLinearUnits = false;
1148 33 : if (EQUAL(pszTemp, "SET_FROM_SRS"))
1149 : {
1150 0 : bSetLinearUnits = true;
1151 : }
1152 33 : else if (STARTS_WITH_CI(pszTemp, "METER"))
1153 : {
1154 32 : nLinearUnits = 0;
1155 : }
1156 1 : else if (EQUAL(pszTemp, "FOOT") || EQUAL(pszTemp, "FEET"))
1157 : {
1158 1 : nLinearUnits = 1;
1159 : }
1160 0 : else if (STARTS_WITH_CI(pszTemp, "KILOMETER"))
1161 : {
1162 0 : nLinearUnits = 2;
1163 : }
1164 33 : bool bCalculateLatitude = true;
1165 33 : int nLatitude = 0;
1166 33 : if (CSLFetchNameValue(papszOptions, "LATITUDE") != nullptr)
1167 : {
1168 33 : bCalculateLatitude = false;
1169 33 : nLatitude = atoi(CSLFetchNameValue(papszOptions, "LATITUDE"));
1170 33 : if (nLatitude > 90 || nLatitude < -90)
1171 : {
1172 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1173 : "Invalid value (%d) for LATITUDE.", nLatitude);
1174 0 : return nullptr;
1175 : }
1176 : }
1177 :
1178 : // If no latitude is supplied, attempt to extract the central latitude
1179 : // from the image. It must be set either manually or here, otherwise
1180 : // we fail.
1181 33 : GDALGeoTransform srcGT;
1182 33 : poSrcDS->GetGeoTransform(srcGT);
1183 33 : const OGRSpatialReference *poSrcSRS = poSrcDS->GetSpatialRef();
1184 33 : double dfLongitude = 0.0;
1185 33 : double dfLatitude = 0.0;
1186 :
1187 33 : const int nYSize = poSrcDS->GetRasterYSize();
1188 :
1189 33 : if (!bCalculateLatitude)
1190 : {
1191 33 : dfLatitude = nLatitude;
1192 : }
1193 0 : else if (poSrcSRS)
1194 : {
1195 0 : OGRSpatialReference oDstSRS;
1196 0 : oDstSRS.importFromEPSG(4269);
1197 0 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1198 : OGRCoordinateTransformation *poCT =
1199 0 : reinterpret_cast<OGRCoordinateTransformation *>(
1200 : OGRCreateCoordinateTransformation(poSrcSRS, &oDstSRS));
1201 0 : if (poCT != nullptr)
1202 : {
1203 0 : dfLatitude = srcGT[3] + srcGT[5] * nYSize / 2;
1204 : const int nErr =
1205 0 : static_cast<int>(poCT->Transform(1, &dfLongitude, &dfLatitude));
1206 0 : if (!nErr)
1207 : {
1208 0 : dfLatitude = 0.0;
1209 : // For the most part, this is an invalid LCP, but it is a
1210 : // changeable value in Flammap/Farsite, etc. We should
1211 : // probably be strict here all the time.
1212 0 : CPLError(CE_Failure, CPLE_AppDefined,
1213 : "Could not calculate latitude from spatial "
1214 : "reference and LATITUDE was not set.");
1215 0 : return nullptr;
1216 : }
1217 : }
1218 0 : OGRCoordinateTransformation::DestroyCT(poCT);
1219 : }
1220 : else
1221 : {
1222 : // See comment above on failure to transform.
1223 0 : CPLError(CE_Failure, CPLE_AppDefined,
1224 : "Could not calculate latitude from spatial reference "
1225 : "and LATITUDE was not set.");
1226 0 : return nullptr;
1227 : }
1228 : // Set the linear units if the metadata item was not already set, and we
1229 : // have an SRS.
1230 33 : if (bSetLinearUnits && poSrcSRS)
1231 : {
1232 0 : const char *pszUnit = poSrcSRS->GetAttrValue("UNIT", 0);
1233 0 : if (pszUnit == nullptr)
1234 : {
1235 0 : if (bStrict)
1236 : {
1237 0 : CPLError(CE_Failure, CPLE_AppDefined,
1238 : "Could not parse linear unit.");
1239 0 : return nullptr;
1240 : }
1241 : else
1242 : {
1243 0 : CPLError(CE_Warning, CPLE_AppDefined,
1244 : "Could not parse linear unit, using meters");
1245 0 : nLinearUnits = 0;
1246 : }
1247 : }
1248 : else
1249 : {
1250 0 : CPLDebug("LCP", "Setting linear unit to %s", pszUnit);
1251 0 : if (EQUAL(pszUnit, "meter") || EQUAL(pszUnit, "metre"))
1252 : {
1253 0 : nLinearUnits = 0;
1254 : }
1255 0 : else if (EQUAL(pszUnit, "feet") || EQUAL(pszUnit, "foot"))
1256 : {
1257 0 : nLinearUnits = 1;
1258 : }
1259 0 : else if (STARTS_WITH_CI(pszUnit, "kilomet"))
1260 : {
1261 0 : nLinearUnits = 2;
1262 : }
1263 : else
1264 : {
1265 0 : if (bStrict)
1266 0 : nLinearUnits = 0;
1267 : }
1268 0 : pszUnit = poSrcSRS->GetAttrValue("UNIT", 1);
1269 0 : if (pszUnit != nullptr)
1270 : {
1271 0 : const double dfScale = CPLAtof(pszUnit);
1272 0 : if (dfScale != 1.0)
1273 : {
1274 0 : if (bStrict)
1275 : {
1276 0 : CPLError(CE_Failure, CPLE_AppDefined,
1277 : "Unit scale is %lf (!=1.0). It is not "
1278 : "supported.",
1279 : dfScale);
1280 0 : return nullptr;
1281 : }
1282 : else
1283 : {
1284 0 : CPLError(CE_Warning, CPLE_AppDefined,
1285 : "Unit scale is %lf (!=1.0). It is not "
1286 : "supported, ignoring.",
1287 : dfScale);
1288 : }
1289 : }
1290 : }
1291 0 : }
1292 : }
1293 33 : else if (bSetLinearUnits)
1294 : {
1295 : // This can be defaulted if it is not a strict creation.
1296 0 : if (bStrict)
1297 : {
1298 0 : CPLError(CE_Failure, CPLE_AppDefined,
1299 : "Could not parse linear unit from spatial reference "
1300 : "and LINEAR_UNIT was not set.");
1301 0 : return nullptr;
1302 : }
1303 : else
1304 : {
1305 0 : CPLError(CE_Warning, CPLE_AppDefined,
1306 : "Could not parse linear unit from spatial reference "
1307 : "and LINEAR_UNIT was not set, defaulting to meters.");
1308 0 : nLinearUnits = 0;
1309 : }
1310 : }
1311 :
1312 33 : const char *pszDescription = CSLFetchNameValueDef(
1313 : papszOptions, "DESCRIPTION", "LCP file created by GDAL.");
1314 :
1315 : // Loop through and get the stats for the bands if we need to calculate
1316 : // them. This probably should be done when we copy the data over to the
1317 : // destination dataset, since we load the values into memory, but this is
1318 : // much simpler code using GDALRasterBand->GetStatistics(). We also may
1319 : // need to classify the data (number of unique values and a list of those
1320 : // values if the number of unique values is > 100. It is currently unclear
1321 : // how these data are used though, so we will implement that at some point
1322 : // if need be.
1323 33 : double *padfMin = static_cast<double *>(CPLMalloc(sizeof(double) * nBands));
1324 33 : double *padfMax = static_cast<double *>(CPLMalloc(sizeof(double) * nBands));
1325 :
1326 : // Initialize these arrays to zeros
1327 : GInt32 *panFound =
1328 33 : static_cast<GInt32 *>(VSIMalloc2(sizeof(GInt32), nBands));
1329 33 : memset(panFound, 0, sizeof(GInt32) * nBands);
1330 : GInt32 *panClasses = static_cast<GInt32 *>(
1331 33 : VSIMalloc3(sizeof(GInt32), nBands, LCP_MAX_CLASSES));
1332 33 : memset(panClasses, 0, sizeof(GInt32) * nBands * LCP_MAX_CLASSES);
1333 :
1334 33 : if (bCalculateStats)
1335 : {
1336 :
1337 353 : for (int i = 0; i < nBands; i++)
1338 : {
1339 320 : GDALRasterBand *poBand = poSrcDS->GetRasterBand(i + 1);
1340 320 : double dfDummy = 0.0;
1341 640 : CPLErr eErr = poBand->GetStatistics(
1342 320 : FALSE, TRUE, &padfMin[i], &padfMax[i], &dfDummy, &dfDummy);
1343 320 : if (eErr != CE_None)
1344 : {
1345 0 : CPLError(CE_Warning, CPLE_AppDefined,
1346 : "Failed to properly "
1347 : "calculate statistics "
1348 : "on band %d",
1349 : i);
1350 0 : padfMin[i] = 0.0;
1351 0 : padfMax[i] = 0.0;
1352 : }
1353 :
1354 : // See comment above.
1355 320 : if (bClassifyData)
1356 : {
1357 640 : eErr = ClassifyBandData(poBand, panFound[i],
1358 320 : panClasses + (i * LCP_MAX_CLASSES));
1359 320 : if (eErr != CE_None)
1360 : {
1361 0 : CPLError(CE_Warning, CPLE_AppDefined,
1362 : "Failed to classify band data on band %d.", i);
1363 : }
1364 : }
1365 : }
1366 : }
1367 :
1368 33 : VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
1369 33 : if (fp == nullptr)
1370 : {
1371 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Unable to create lcp file %s.",
1372 : pszFilename);
1373 0 : CPLFree(padfMin);
1374 0 : CPLFree(padfMax);
1375 0 : CPLFree(panFound);
1376 0 : CPLFree(panClasses);
1377 0 : return nullptr;
1378 : }
1379 :
1380 : /* -------------------------------------------------------------------- */
1381 : /* Write the header */
1382 : /* -------------------------------------------------------------------- */
1383 :
1384 33 : GInt32 nTemp = bHaveCrownFuels ? 21 : 20;
1385 33 : CPL_LSBPTR32(&nTemp);
1386 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&nTemp, 4, 1, fp));
1387 33 : nTemp = bHaveGroundFuels ? 21 : 20;
1388 33 : CPL_LSBPTR32(&nTemp);
1389 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&nTemp, 4, 1, fp));
1390 :
1391 33 : const int nXSize = poSrcDS->GetRasterXSize();
1392 33 : nTemp = static_cast<GInt32>(dfLatitude + 0.5);
1393 33 : CPL_LSBPTR32(&nTemp);
1394 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&nTemp, 4, 1, fp));
1395 33 : dfLongitude = srcGT[0] + srcGT[1] * nXSize;
1396 33 : CPL_LSBPTR64(&dfLongitude);
1397 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfLongitude, 8, 1, fp));
1398 33 : dfLongitude = srcGT[0];
1399 33 : CPL_LSBPTR64(&dfLongitude);
1400 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfLongitude, 8, 1, fp));
1401 33 : dfLatitude = srcGT[3];
1402 33 : CPL_LSBPTR64(&dfLatitude);
1403 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfLatitude, 8, 1, fp));
1404 33 : dfLatitude = srcGT[3] + srcGT[5] * nYSize;
1405 33 : CPL_LSBPTR64(&dfLatitude);
1406 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfLatitude, 8, 1, fp));
1407 :
1408 : // Swap the two classification arrays if we are writing them, and they need
1409 : // to be swapped.
1410 : #ifdef CPL_MSB
1411 : if (bClassifyData)
1412 : {
1413 : GDALSwapWords(panFound, 2, nBands, 2);
1414 : GDALSwapWords(panClasses, 2, LCP_MAX_CLASSES, 2);
1415 : }
1416 : #endif
1417 :
1418 33 : if (bCalculateStats)
1419 : {
1420 353 : for (int i = 0; i < nBands; i++)
1421 : {
1422 : // If we don't have Crown fuels, but do have Ground fuels, we
1423 : // have to 'fast forward'.
1424 320 : if (i == 5 && !bHaveCrownFuels && bHaveGroundFuels)
1425 : {
1426 1 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 3340, SEEK_SET));
1427 : }
1428 320 : nTemp = static_cast<GInt32>(padfMin[i]);
1429 320 : CPL_LSBPTR32(&nTemp);
1430 320 : CPL_IGNORE_RET_VAL(VSIFWriteL(&nTemp, 4, 1, fp));
1431 320 : nTemp = static_cast<GInt32>(padfMax[i]);
1432 320 : CPL_LSBPTR32(&nTemp);
1433 320 : CPL_IGNORE_RET_VAL(VSIFWriteL(&nTemp, 4, 1, fp));
1434 320 : if (bClassifyData)
1435 : {
1436 : // These two arrays were swapped in their entirety above.
1437 320 : CPL_IGNORE_RET_VAL(VSIFWriteL(panFound + i, 4, 1, fp));
1438 320 : CPL_IGNORE_RET_VAL(
1439 320 : VSIFWriteL(panClasses + (i * LCP_MAX_CLASSES), 4,
1440 : LCP_MAX_CLASSES, fp));
1441 : }
1442 : else
1443 : {
1444 0 : nTemp = -1;
1445 0 : CPL_LSBPTR32(&nTemp);
1446 0 : CPL_IGNORE_RET_VAL(VSIFWriteL(&nTemp, 4, 1, fp));
1447 0 : CPL_IGNORE_RET_VAL(
1448 0 : VSIFSeekL(fp, 4 * LCP_MAX_CLASSES, SEEK_CUR));
1449 : }
1450 : }
1451 : }
1452 : else
1453 : {
1454 0 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 4164, SEEK_SET));
1455 : }
1456 33 : CPLFree(padfMin);
1457 33 : CPLFree(padfMax);
1458 33 : CPLFree(panFound);
1459 33 : CPLFree(panClasses);
1460 :
1461 : // Should be at one of 3 locations, 2104, 3340, or 4164.
1462 33 : CPLAssert(VSIFTellL(fp) == 2104 || VSIFTellL(fp) == 3340 ||
1463 : VSIFTellL(fp) == 4164);
1464 33 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 4164, SEEK_SET));
1465 :
1466 : /* Image size */
1467 33 : nTemp = static_cast<GInt32>(nXSize);
1468 33 : CPL_LSBPTR32(&nTemp);
1469 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&nTemp, 4, 1, fp));
1470 33 : nTemp = static_cast<GInt32>(nYSize);
1471 33 : CPL_LSBPTR32(&nTemp);
1472 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&nTemp, 4, 1, fp));
1473 :
1474 : // X and Y boundaries.
1475 : // Max x.
1476 33 : double dfTemp = srcGT[0] + srcGT[1] * nXSize;
1477 33 : CPL_LSBPTR64(&dfTemp);
1478 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfTemp, 8, 1, fp));
1479 : // Min x.
1480 33 : dfTemp = srcGT[0];
1481 33 : CPL_LSBPTR64(&dfTemp);
1482 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfTemp, 8, 1, fp));
1483 : // Max y.
1484 33 : dfTemp = srcGT[3];
1485 33 : CPL_LSBPTR64(&dfTemp);
1486 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfTemp, 8, 1, fp));
1487 : // Min y.
1488 33 : dfTemp = srcGT[3] + srcGT[5] * nYSize;
1489 33 : CPL_LSBPTR64(&dfTemp);
1490 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfTemp, 8, 1, fp));
1491 :
1492 33 : nTemp = nLinearUnits;
1493 33 : CPL_LSBPTR32(&nTemp);
1494 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&nTemp, 4, 1, fp));
1495 :
1496 : // Resolution.
1497 : // X resolution.
1498 33 : dfTemp = srcGT[1];
1499 33 : CPL_LSBPTR64(&dfTemp);
1500 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfTemp, 8, 1, fp));
1501 : // Y resolution.
1502 33 : dfTemp = fabs(srcGT[5]);
1503 33 : CPL_LSBPTR64(&dfTemp);
1504 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfTemp, 8, 1, fp));
1505 :
1506 : #ifdef CPL_MSB
1507 : GDALSwapWords(panMetadata, 2, LCP_MAX_BANDS, 2);
1508 : #endif
1509 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(panMetadata, 2, LCP_MAX_BANDS, fp));
1510 :
1511 : // Write the source filenames.
1512 33 : char **papszFileList = poSrcDS->GetFileList();
1513 33 : if (papszFileList != nullptr)
1514 : {
1515 308 : for (int i = 0; i < nBands; i++)
1516 : {
1517 280 : if (i == 5 && !bHaveCrownFuels && bHaveGroundFuels)
1518 : {
1519 0 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 6292, SEEK_SET));
1520 : }
1521 280 : CPL_IGNORE_RET_VAL(
1522 280 : VSIFWriteL(papszFileList[0], 1,
1523 : CPLStrnlen(papszFileList[0], LCP_MAX_PATH), fp));
1524 280 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 4244 + (256 * (i + 1)), SEEK_SET));
1525 : }
1526 : }
1527 : // No file list, mem driver, etc.
1528 : else
1529 : {
1530 5 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 6804, SEEK_SET));
1531 : }
1532 33 : CSLDestroy(papszFileList);
1533 : // Should be at location 5524, 6292 or 6804.
1534 33 : CPLAssert(VSIFTellL(fp) == 5524 || VSIFTellL(fp) == 6292 ||
1535 : VSIFTellL(fp) == 6804);
1536 33 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 6804, SEEK_SET));
1537 :
1538 : // Description .
1539 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(
1540 : pszDescription, 1, CPLStrnlen(pszDescription, LCP_MAX_DESC), fp));
1541 :
1542 : // Should be at or below location 7316, all done with the header.
1543 33 : CPLAssert(VSIFTellL(fp) <= 7316);
1544 33 : CPL_IGNORE_RET_VAL(VSIFSeekL(fp, 7316, SEEK_SET));
1545 :
1546 : /* -------------------------------------------------------------------- */
1547 : /* Loop over image, copying image data. */
1548 : /* -------------------------------------------------------------------- */
1549 :
1550 33 : GInt16 *panScanline = static_cast<GInt16 *>(VSIMalloc3(2, nBands, nXSize));
1551 :
1552 33 : if (!pfnProgress(0.0, nullptr, pProgressData))
1553 : {
1554 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1555 0 : VSIFree(panScanline);
1556 0 : return nullptr;
1557 : }
1558 399 : for (int iLine = 0; iLine < nYSize; iLine++)
1559 : {
1560 3826 : for (int iBand = 0; iBand < nBands; iBand++)
1561 : {
1562 3460 : GDALRasterBand *poBand = poSrcDS->GetRasterBand(iBand + 1);
1563 6920 : CPLErr eErr = poBand->RasterIO(
1564 3460 : GF_Read, 0, iLine, nXSize, 1, panScanline + iBand, nXSize, 1,
1565 3460 : GDT_Int16, nBands * 2, static_cast<size_t>(nBands) * nXSize * 2,
1566 : nullptr);
1567 : // Not sure what to do here.
1568 3460 : if (eErr != CE_None)
1569 : {
1570 0 : CPLError(CE_Warning, CPLE_AppDefined,
1571 : "Error reported in "
1572 : "RasterIO");
1573 : }
1574 : }
1575 : #ifdef CPL_MSB
1576 : GDALSwapWordsEx(panScanline, 2, static_cast<size_t>(nBands) * nXSize,
1577 : 2);
1578 : #endif
1579 366 : CPL_IGNORE_RET_VAL(VSIFWriteL(
1580 366 : panScanline, 2, static_cast<size_t>(nBands) * nXSize, fp));
1581 :
1582 366 : if (!pfnProgress(iLine / static_cast<double>(nYSize), nullptr,
1583 : pProgressData))
1584 : {
1585 0 : VSIFree(reinterpret_cast<void *>(panScanline));
1586 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1587 0 : return nullptr;
1588 : }
1589 : }
1590 33 : VSIFree(panScanline);
1591 33 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1592 33 : if (!pfnProgress(1.0, nullptr, pProgressData))
1593 : {
1594 0 : return nullptr;
1595 : }
1596 :
1597 : // Try to write projection file. *Most* landfire data follows ESRI
1598 : // style projection files, so we use the same code as the AAIGrid driver.
1599 33 : if (poSrcSRS)
1600 : {
1601 0 : char *pszESRIProjection = nullptr;
1602 0 : const char *const apszOptions[] = {"FORMAT=WKT1_ESRI", nullptr};
1603 0 : poSrcSRS->exportToWkt(&pszESRIProjection, apszOptions);
1604 0 : if (pszESRIProjection)
1605 : {
1606 : char *const pszDirname =
1607 0 : CPLStrdup(CPLGetPathSafe(pszFilename).c_str());
1608 : char *const pszBasename =
1609 0 : CPLStrdup(CPLGetBasenameSafe(pszFilename).c_str());
1610 0 : char *pszPrjFilename = CPLStrdup(
1611 0 : CPLFormFilenameSafe(pszDirname, pszBasename, "prj").c_str());
1612 0 : fp = VSIFOpenL(pszPrjFilename, "wt");
1613 0 : if (fp != nullptr)
1614 : {
1615 0 : CPL_IGNORE_RET_VAL(VSIFWriteL(pszESRIProjection, 1,
1616 : strlen(pszESRIProjection), fp));
1617 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1618 : }
1619 : else
1620 : {
1621 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to create file %s.",
1622 : pszPrjFilename);
1623 : }
1624 0 : CPLFree(pszDirname);
1625 0 : CPLFree(pszBasename);
1626 0 : CPLFree(pszPrjFilename);
1627 : }
1628 0 : CPLFree(pszESRIProjection);
1629 : }
1630 33 : return static_cast<GDALDataset *>(GDALOpen(pszFilename, GA_ReadOnly));
1631 : }
1632 :
1633 : /************************************************************************/
1634 : /* GDALRegister_LCP() */
1635 : /************************************************************************/
1636 :
1637 2038 : void GDALRegister_LCP()
1638 :
1639 : {
1640 2038 : if (GDALGetDriverByName("LCP") != nullptr)
1641 283 : return;
1642 :
1643 1755 : GDALDriver *poDriver = new GDALDriver();
1644 :
1645 1755 : poDriver->SetDescription("LCP");
1646 1755 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1647 1755 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1648 1755 : "FARSITE v.4 Landscape File (.lcp)");
1649 1755 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "lcp");
1650 1755 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/lcp.html");
1651 :
1652 1755 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1653 :
1654 1755 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Int16");
1655 1755 : poDriver->SetMetadataItem(
1656 : GDAL_DMD_CREATIONOPTIONLIST,
1657 : "<CreationOptionList>"
1658 : " <Option name='ELEVATION_UNIT' type='string-select' "
1659 : "default='METERS' description='Elevation units'>"
1660 : " <Value>METERS</Value>"
1661 : " <Value>FEET</Value>"
1662 : " </Option>"
1663 : " <Option name='SLOPE_UNIT' type='string-select' default='DEGREES' "
1664 : "description='Slope units'>"
1665 : " <Value>DEGREES</Value>"
1666 : " <Value>PERCENT</Value>"
1667 : " </Option>"
1668 : " <Option name='ASPECT_UNIT' type='string-select' "
1669 : "default='AZIMUTH_DEGREES'>"
1670 : " <Value>GRASS_CATEGORIES</Value>"
1671 : " <Value>AZIMUTH_DEGREES</Value>"
1672 : " <Value>GRASS_DEGREES</Value>"
1673 : " </Option>"
1674 : " <Option name='FUEL_MODEL_OPTION' type='string-select' "
1675 : "default='NO_CUSTOM_AND_NO_FILE'>"
1676 : " <Value>NO_CUSTOM_AND_NO_FILE</Value>"
1677 : " <Value>CUSTOM_AND_NO_FILE</Value>"
1678 : " <Value>NO_CUSTOM_AND_FILE</Value>"
1679 : " <Value>CUSTOM_AND_FILE</Value>"
1680 : " </Option>"
1681 : " <Option name='CANOPY_COV_UNIT' type='string-select' "
1682 : "default='PERCENT'>"
1683 : " <Value>CATEGORIES</Value>"
1684 : " <Value>PERCENT</Value>"
1685 : " </Option>"
1686 : " <Option name='CANOPY_HT_UNIT' type='string-select' "
1687 : "default='METERS_X_10'>"
1688 : " <Value>METERS</Value>"
1689 : " <Value>FEET</Value>"
1690 : " <Value>METERS_X_10</Value>"
1691 : " <Value>FEET_X_10</Value>"
1692 : " </Option>"
1693 : " <Option name='CBH_UNIT' type='string-select' default='METERS_X_10'>"
1694 : " <Value>METERS</Value>"
1695 : " <Value>FEET</Value>"
1696 : " <Value>METERS_X_10</Value>"
1697 : " <Value>FEET_X_10</Value>"
1698 : " </Option>"
1699 : " <Option name='CBD_UNIT' type='string-select' "
1700 : "default='KG_PER_CUBIC_METER_X_100'>"
1701 : " <Value>KG_PER_CUBIC_METER</Value>"
1702 : " <Value>POUND_PER_CUBIC_FOOT</Value>"
1703 : " <Value>KG_PER_CUBIC_METER_X_100</Value>"
1704 : " <Value>POUND_PER_CUBIC_FOOT_X_1000</Value>"
1705 : " </Option>"
1706 : " <Option name='DUFF_UNIT' type='string-select' "
1707 : "default='MG_PER_HECTARE_X_10'>"
1708 : " <Value>MG_PER_HECTARE_X_10</Value>"
1709 : " <Value>TONS_PER_ACRE_X_10</Value>"
1710 : " </Option>"
1711 : // Kyle does not think we need to override this, but maybe?
1712 : // " <Option name='CWD_OPTION' type='boolean' default='FALSE'
1713 : // description='Override logic for setting the coarse woody presence'/>"
1714 : // */
1715 : " <Option name='CALCULATE_STATS' type='boolean' default='YES' "
1716 : "description='Write the stats to the lcp'/>"
1717 : " <Option name='CLASSIFY_DATA' type='boolean' default='YES' "
1718 : "description='Write the stats to the lcp'/>"
1719 : " <Option name='LINEAR_UNIT' type='string-select' "
1720 : "default='SET_FROM_SRS' description='Set the linear units in the lcp'>"
1721 : " <Value>SET_FROM_SRS</Value>"
1722 : " <Value>METER</Value>"
1723 : " <Value>FOOT</Value>"
1724 : " <Value>KILOMETER</Value>"
1725 : " </Option>"
1726 : " <Option name='LATITUDE' type='int' default='' description='Set the "
1727 : "latitude for the dataset, this overrides the driver trying to set it "
1728 : "programmatically in EPSG:4269'/>"
1729 : " <Option name='DESCRIPTION' type='string' default='LCP file created "
1730 : "by GDAL' description='A short description of the lcp file'/>"
1731 1755 : "</CreationOptionList>");
1732 :
1733 1755 : poDriver->pfnOpen = LCPDataset::Open;
1734 1755 : poDriver->pfnCreateCopy = LCPDataset::CreateCopy;
1735 1755 : poDriver->pfnIdentify = LCPDataset::Identify;
1736 :
1737 1755 : GetGDALDriverManager()->RegisterDriver(poDriver);
1738 : }
|
