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