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 51708 : int LCPDataset::Identify(GDALOpenInfo *poOpenInfo)
161 :
162 : {
163 : /* -------------------------------------------------------------------- */
164 : /* Verify that this is a FARSITE v.4 LCP file */
165 : /* -------------------------------------------------------------------- */
166 51708 : if (poOpenInfo->nHeaderBytes < 50)
167 48121 : return FALSE;
168 :
169 : /* check if first three fields have valid data */
170 3587 : if ((CPL_LSBSINT32PTR(poOpenInfo->pabyHeader) != 20 &&
171 3573 : 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 3419 : return FALSE;
178 : }
179 :
180 : /* -------------------------------------------------------------------- */
181 : /* Check file extension */
182 : /* -------------------------------------------------------------------- */
183 : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
184 168 : const char *pszFileExtension = CPLGetExtension(poOpenInfo->pszFilename);
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 268 : for (int i = 0; i <= nTemp; i++)
461 : {
462 193 : const int nTemp2 = CPL_LSBSINT32PTR(poDS->pachHeader +
463 : (1292 + (i * 4)));
464 193 : 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 75 : char *const pszDirname = CPLStrdup(CPLGetPath(poOpenInfo->pszFilename));
744 : char *const pszBasename =
745 75 : CPLStrdup(CPLGetBasename(poOpenInfo->pszFilename));
746 :
747 75 : poDS->osPrjFilename = CPLFormFilename(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 = CPLFormFilename(pszDirname, pszBasename, "PRJ");
754 72 : nRet = VSIStatL(poDS->osPrjFilename, &sStatBuf);
755 : }
756 :
757 75 : if (nRet == 0)
758 : {
759 3 : char **papszPrj = CSLLoad(poDS->osPrjFilename);
760 :
761 3 : CPLDebug("LCP", "Loaded SRS from %s", poDS->osPrjFilename.c_str());
762 :
763 3 : poDS->m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
764 3 : if (poDS->m_oSRS.importFromESRI(papszPrj) != OGRERR_NONE)
765 : {
766 0 : poDS->m_oSRS.Clear();
767 : }
768 :
769 3 : CSLDestroy(papszPrj);
770 : }
771 :
772 75 : CPLFree(pszDirname);
773 75 : CPLFree(pszBasename);
774 :
775 : /* -------------------------------------------------------------------- */
776 : /* Initialize any PAM information. */
777 : /* -------------------------------------------------------------------- */
778 75 : poDS->SetDescription(poOpenInfo->pszFilename);
779 75 : poDS->TryLoadXML();
780 :
781 : /* -------------------------------------------------------------------- */
782 : /* Check for external overviews. */
783 : /* -------------------------------------------------------------------- */
784 150 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename,
785 75 : poOpenInfo->GetSiblingFiles());
786 :
787 75 : return poDS.release();
788 : }
789 :
790 : /************************************************************************/
791 : /* ClassifyBandData() */
792 : /* Classify a band and store 99 or fewer unique values. If there are */
793 : /* more than 99 unique values, then set pnNumClasses to -1 as a flag */
794 : /* that represents this. These are legacy values in the header, and */
795 : /* while we should never deprecate them, we could possibly not */
796 : /* calculate them by default. */
797 : /************************************************************************/
798 :
799 320 : CPLErr LCPDataset::ClassifyBandData(GDALRasterBand *poBand, GInt32 &nNumClasses,
800 : GInt32 *panClasses)
801 : {
802 320 : CPLAssert(poBand);
803 320 : CPLAssert(panClasses);
804 :
805 320 : const int nXSize = poBand->GetXSize();
806 320 : const int nYSize = poBand->GetYSize();
807 :
808 : GInt16 *panValues =
809 320 : static_cast<GInt16 *>(CPLMalloc(sizeof(GInt16) * nXSize));
810 320 : constexpr int MIN_VAL = std::numeric_limits<GInt16>::min();
811 320 : constexpr int MAX_VAL = std::numeric_limits<GInt16>::max();
812 320 : constexpr int RANGE_VAL = MAX_VAL - MIN_VAL + 1;
813 320 : GByte *pabyFlags = static_cast<GByte *>(CPLCalloc(1, RANGE_VAL));
814 :
815 : /* so that values in [-32768,32767] are mapped to [0,65535] in pabyFlags */
816 320 : constexpr int OFFSET = -MIN_VAL;
817 :
818 320 : int nFound = 0;
819 320 : bool bTooMany = false;
820 320 : CPLErr eErr = CE_None;
821 3770 : for (int iLine = 0; iLine < nYSize; iLine++)
822 : {
823 3451 : eErr = poBand->RasterIO(GF_Read, 0, iLine, nXSize, 1, panValues, nXSize,
824 : 1, GDT_Int16, 0, 0, nullptr);
825 3451 : if (eErr != CE_None)
826 0 : break;
827 37531 : for (int iPixel = 0; iPixel < nXSize; iPixel++)
828 : {
829 34081 : if (panValues[iPixel] == -9999)
830 : {
831 0 : continue;
832 : }
833 34081 : if (nFound == LCP_MAX_CLASSES)
834 : {
835 1 : CPLDebug("LCP",
836 : "Found more that %d unique values in "
837 : "band %d. Not 'classifying' the data.",
838 : LCP_MAX_CLASSES - 1, poBand->GetBand());
839 1 : nFound = -1;
840 1 : bTooMany = true;
841 1 : break;
842 : }
843 34080 : if (pabyFlags[panValues[iPixel] + OFFSET] == 0)
844 : {
845 577 : pabyFlags[panValues[iPixel] + OFFSET] = 1;
846 577 : nFound++;
847 : }
848 : }
849 3451 : if (bTooMany)
850 1 : break;
851 : }
852 320 : if (!bTooMany)
853 : {
854 : // The classes are always padded with a leading 0. This was for aligning
855 : // offsets, or making it a 1-based array instead of 0-based.
856 319 : panClasses[0] = 0;
857 20906300 : for (int j = 0, nIndex = 1; j < RANGE_VAL; j++)
858 : {
859 20906000 : if (pabyFlags[j] == 1)
860 : {
861 477 : panClasses[nIndex++] = j;
862 : }
863 : }
864 : }
865 320 : nNumClasses = nFound;
866 320 : CPLFree(pabyFlags);
867 320 : CPLFree(panValues);
868 :
869 320 : return eErr;
870 : }
871 :
872 : /************************************************************************/
873 : /* CreateCopy() */
874 : /************************************************************************/
875 :
876 68 : GDALDataset *LCPDataset::CreateCopy(const char *pszFilename,
877 : GDALDataset *poSrcDS, int bStrict,
878 : char **papszOptions,
879 : GDALProgressFunc pfnProgress,
880 : void *pProgressData)
881 :
882 : {
883 : /* -------------------------------------------------------------------- */
884 : /* Verify input options. */
885 : /* -------------------------------------------------------------------- */
886 68 : const int nBands = poSrcDS->GetRasterCount();
887 68 : if (nBands != 5 && nBands != 7 && nBands != 8 && nBands != 10)
888 : {
889 25 : CPLError(CE_Failure, CPLE_NotSupported,
890 : "LCP driver doesn't support %d bands. Must be 5, 7, 8 "
891 : "or 10 bands.",
892 : nBands);
893 25 : return nullptr;
894 : }
895 :
896 43 : GDALDataType eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
897 43 : if (eType != GDT_Int16 && bStrict)
898 : {
899 1 : CPLError(CE_Failure, CPLE_AppDefined,
900 : "LCP only supports 16-bit signed integer data types.");
901 1 : return nullptr;
902 : }
903 42 : else if (eType != GDT_Int16)
904 : {
905 0 : CPLError(CE_Warning, CPLE_AppDefined,
906 : "Setting data type to 16-bit integer.");
907 : }
908 :
909 : /* -------------------------------------------------------------------- */
910 : /* What schema do we have (ground/crown fuels) */
911 : /* -------------------------------------------------------------------- */
912 :
913 42 : const bool bHaveCrownFuels = nBands == 8 || nBands == 10;
914 42 : const bool bHaveGroundFuels = nBands == 7 || nBands == 10;
915 :
916 : // Since units are 'configurable', we should check for user
917 : // defined units. This is a bit cumbersome, but the user should
918 : // be allowed to specify none to get default units/options. Use
919 : // default units every chance we get.
920 42 : GInt16 panMetadata[LCP_MAX_BANDS] = {
921 : 0, // 0 ELEVATION_UNIT
922 : 0, // 1 SLOPE_UNIT
923 : 2, // 2 ASPECT_UNIT
924 : 0, // 3 FUEL_MODEL_OPTION
925 : 1, // 4 CANOPY_COV_UNIT
926 : 3, // 5 CANOPY_HT_UNIT
927 : 3, // 6 CBH_UNIT
928 : 3, // 7 CBD_UNIT
929 : 1, // 8 DUFF_UNIT
930 : 0, // 9 CWD_OPTION
931 : };
932 :
933 : // Check the units/options for user overrides.
934 : const char *pszTemp =
935 42 : CSLFetchNameValueDef(papszOptions, "ELEVATION_UNIT", "METERS");
936 42 : if (STARTS_WITH_CI(pszTemp, "METER"))
937 : {
938 40 : panMetadata[0] = 0;
939 : }
940 2 : else if (EQUAL(pszTemp, "FEET") || EQUAL(pszTemp, "FOOT"))
941 : {
942 1 : panMetadata[0] = 1;
943 : }
944 : else
945 : {
946 1 : CPLError(CE_Failure, CPLE_AppDefined,
947 : "Invalid value (%s) for ELEVATION_UNIT.", pszTemp);
948 1 : return nullptr;
949 : }
950 :
951 41 : pszTemp = CSLFetchNameValueDef(papszOptions, "SLOPE_UNIT", "DEGREES");
952 41 : if (EQUAL(pszTemp, "DEGREES"))
953 : {
954 39 : panMetadata[1] = 0;
955 : }
956 2 : else if (EQUAL(pszTemp, "PERCENT"))
957 : {
958 1 : panMetadata[1] = 1;
959 : }
960 : else
961 : {
962 1 : CPLError(CE_Failure, CPLE_AppDefined,
963 : "Invalid value (%s) for SLOPE_UNIT.", pszTemp);
964 1 : return nullptr;
965 : }
966 :
967 : pszTemp =
968 40 : CSLFetchNameValueDef(papszOptions, "ASPECT_UNIT", "AZIMUTH_DEGREES");
969 40 : if (EQUAL(pszTemp, "GRASS_CATEGORIES"))
970 : {
971 1 : panMetadata[2] = 0;
972 : }
973 39 : else if (EQUAL(pszTemp, "GRASS_DEGREES"))
974 : {
975 1 : panMetadata[2] = 1;
976 : }
977 38 : else if (EQUAL(pszTemp, "AZIMUTH_DEGREES"))
978 : {
979 37 : panMetadata[2] = 2;
980 : }
981 : else
982 : {
983 1 : CPLError(CE_Failure, CPLE_AppDefined,
984 : "Invalid value (%s) for ASPECT_UNIT.", pszTemp);
985 1 : return nullptr;
986 : }
987 :
988 39 : pszTemp = CSLFetchNameValueDef(papszOptions, "FUEL_MODEL_OPTION",
989 : "NO_CUSTOM_AND_NO_FILE");
990 39 : if (EQUAL(pszTemp, "NO_CUSTOM_AND_NO_FILE"))
991 : {
992 38 : panMetadata[3] = 0;
993 : }
994 1 : else if (EQUAL(pszTemp, "CUSTOM_AND_NO_FILE"))
995 : {
996 0 : panMetadata[3] = 1;
997 : }
998 1 : else if (EQUAL(pszTemp, "NO_CUSTOM_AND_FILE"))
999 : {
1000 0 : panMetadata[3] = 2;
1001 : }
1002 1 : else if (EQUAL(pszTemp, "CUSTOM_AND_FILE"))
1003 : {
1004 0 : panMetadata[3] = 3;
1005 : }
1006 : else
1007 : {
1008 1 : CPLError(CE_Failure, CPLE_AppDefined,
1009 : "Invalid value (%s) for FUEL_MODEL_OPTION.", pszTemp);
1010 1 : return nullptr;
1011 : }
1012 :
1013 38 : pszTemp = CSLFetchNameValueDef(papszOptions, "CANOPY_COV_UNIT", "PERCENT");
1014 38 : if (EQUAL(pszTemp, "CATEGORIES"))
1015 : {
1016 1 : panMetadata[4] = 0;
1017 : }
1018 37 : else if (EQUAL(pszTemp, "PERCENT"))
1019 : {
1020 36 : panMetadata[4] = 1;
1021 : }
1022 : else
1023 : {
1024 1 : CPLError(CE_Failure, CPLE_AppDefined,
1025 : "Invalid value (%s) for CANOPY_COV_UNIT.", pszTemp);
1026 1 : return nullptr;
1027 : }
1028 :
1029 37 : if (bHaveCrownFuels)
1030 : {
1031 : pszTemp =
1032 35 : CSLFetchNameValueDef(papszOptions, "CANOPY_HT_UNIT", "METERS_X_10");
1033 35 : if (EQUAL(pszTemp, "METERS") || EQUAL(pszTemp, "METER"))
1034 : {
1035 1 : panMetadata[5] = 1;
1036 : }
1037 34 : else if (EQUAL(pszTemp, "FEET") || EQUAL(pszTemp, "FOOT"))
1038 : {
1039 1 : panMetadata[5] = 2;
1040 : }
1041 33 : else if (EQUAL(pszTemp, "METERS_X_10") || EQUAL(pszTemp, "METER_X_10"))
1042 : {
1043 31 : panMetadata[5] = 3;
1044 : }
1045 2 : else if (EQUAL(pszTemp, "FEET_X_10") || EQUAL(pszTemp, "FOOT_X_10"))
1046 : {
1047 1 : panMetadata[5] = 4;
1048 : }
1049 : else
1050 : {
1051 1 : CPLError(CE_Failure, CPLE_AppDefined,
1052 : "Invalid value (%s) for CANOPY_HT_UNIT.", pszTemp);
1053 1 : return nullptr;
1054 : }
1055 :
1056 34 : pszTemp = CSLFetchNameValueDef(papszOptions, "CBH_UNIT", "METERS_X_10");
1057 34 : if (EQUAL(pszTemp, "METERS") || EQUAL(pszTemp, "METER"))
1058 : {
1059 1 : panMetadata[6] = 1;
1060 : }
1061 33 : else if (EQUAL(pszTemp, "FEET") || EQUAL(pszTemp, "FOOT"))
1062 : {
1063 1 : panMetadata[6] = 2;
1064 : }
1065 32 : else if (EQUAL(pszTemp, "METERS_X_10") || EQUAL(pszTemp, "METER_X_10"))
1066 : {
1067 30 : panMetadata[6] = 3;
1068 : }
1069 2 : else if (EQUAL(pszTemp, "FEET_X_10") || EQUAL(pszTemp, "FOOT_X_10"))
1070 : {
1071 1 : panMetadata[6] = 4;
1072 : }
1073 : else
1074 : {
1075 1 : CPLError(CE_Failure, CPLE_AppDefined,
1076 : "Invalid value (%s) for CBH_UNIT.", pszTemp);
1077 1 : return nullptr;
1078 : }
1079 :
1080 33 : pszTemp = CSLFetchNameValueDef(papszOptions, "CBD_UNIT",
1081 : "KG_PER_CUBIC_METER_X_100");
1082 33 : if (EQUAL(pszTemp, "KG_PER_CUBIC_METER"))
1083 : {
1084 1 : panMetadata[7] = 1;
1085 : }
1086 32 : else if (EQUAL(pszTemp, "POUND_PER_CUBIC_FOOT"))
1087 : {
1088 1 : panMetadata[7] = 2;
1089 : }
1090 31 : else if (EQUAL(pszTemp, "KG_PER_CUBIC_METER_X_100"))
1091 : {
1092 29 : panMetadata[7] = 3;
1093 : }
1094 2 : else if (EQUAL(pszTemp, "POUND_PER_CUBIC_FOOT_X_1000"))
1095 : {
1096 1 : panMetadata[7] = 4;
1097 : }
1098 : else
1099 : {
1100 1 : CPLError(CE_Failure, CPLE_AppDefined,
1101 : "Invalid value (%s) for CBD_UNIT.", pszTemp);
1102 1 : return nullptr;
1103 : }
1104 : }
1105 :
1106 34 : if (bHaveGroundFuels)
1107 : {
1108 32 : pszTemp = CSLFetchNameValueDef(papszOptions, "DUFF_UNIT",
1109 : "MG_PER_HECTARE_X_10");
1110 32 : if (EQUAL(pszTemp, "MG_PER_HECTARE_X_10"))
1111 : {
1112 30 : panMetadata[8] = 1;
1113 : }
1114 2 : else if (EQUAL(pszTemp, "TONS_PER_ACRE_X_10"))
1115 : {
1116 1 : panMetadata[8] = 2;
1117 : }
1118 : else
1119 : {
1120 1 : CPLError(CE_Failure, CPLE_AppDefined,
1121 : "Invalid value (%s) for DUFF_UNIT.", pszTemp);
1122 1 : return nullptr;
1123 : }
1124 :
1125 31 : panMetadata[9] = 1;
1126 : }
1127 :
1128 : // Calculate the stats for each band. The binary file carries along
1129 : // these metadata for display purposes(?).
1130 33 : bool bCalculateStats = CPLFetchBool(papszOptions, "CALCULATE_STATS", true);
1131 :
1132 : const bool bClassifyData =
1133 33 : CPLFetchBool(papszOptions, "CLASSIFY_DATA", true);
1134 :
1135 : // We should have stats if we classify, we'll get them anyway.
1136 33 : if (bClassifyData && !bCalculateStats)
1137 : {
1138 0 : CPLError(CE_Warning, CPLE_AppDefined,
1139 : "Ignoring request to not calculate statistics, "
1140 : "because CLASSIFY_DATA was set to ON");
1141 0 : bCalculateStats = true;
1142 : }
1143 :
1144 33 : pszTemp = CSLFetchNameValueDef(papszOptions, "LINEAR_UNIT", "SET_FROM_SRS");
1145 33 : int nLinearUnits = 0;
1146 33 : bool bSetLinearUnits = false;
1147 33 : if (EQUAL(pszTemp, "SET_FROM_SRS"))
1148 : {
1149 0 : bSetLinearUnits = true;
1150 : }
1151 33 : else if (STARTS_WITH_CI(pszTemp, "METER"))
1152 : {
1153 32 : nLinearUnits = 0;
1154 : }
1155 1 : else if (EQUAL(pszTemp, "FOOT") || EQUAL(pszTemp, "FEET"))
1156 : {
1157 1 : nLinearUnits = 1;
1158 : }
1159 0 : else if (STARTS_WITH_CI(pszTemp, "KILOMETER"))
1160 : {
1161 0 : nLinearUnits = 2;
1162 : }
1163 33 : bool bCalculateLatitude = true;
1164 33 : int nLatitude = 0;
1165 33 : if (CSLFetchNameValue(papszOptions, "LATITUDE") != nullptr)
1166 : {
1167 33 : bCalculateLatitude = false;
1168 33 : nLatitude = atoi(CSLFetchNameValue(papszOptions, "LATITUDE"));
1169 33 : if (nLatitude > 90 || nLatitude < -90)
1170 : {
1171 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1172 : "Invalid value (%d) for LATITUDE.", nLatitude);
1173 0 : return nullptr;
1174 : }
1175 : }
1176 :
1177 : // If no latitude is supplied, attempt to extract the central latitude
1178 : // from the image. It must be set either manually or here, otherwise
1179 : // we fail.
1180 33 : double adfSrcGeoTransform[6] = {0.0};
1181 33 : poSrcDS->GetGeoTransform(adfSrcGeoTransform);
1182 33 : const OGRSpatialReference *poSrcSRS = poSrcDS->GetSpatialRef();
1183 33 : double dfLongitude = 0.0;
1184 33 : double dfLatitude = 0.0;
1185 :
1186 33 : const int nYSize = poSrcDS->GetRasterYSize();
1187 :
1188 33 : if (!bCalculateLatitude)
1189 : {
1190 33 : dfLatitude = nLatitude;
1191 : }
1192 0 : else if (poSrcSRS)
1193 : {
1194 0 : OGRSpatialReference oDstSRS;
1195 0 : oDstSRS.importFromEPSG(4269);
1196 0 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1197 : OGRCoordinateTransformation *poCT =
1198 0 : reinterpret_cast<OGRCoordinateTransformation *>(
1199 : OGRCreateCoordinateTransformation(poSrcSRS, &oDstSRS));
1200 0 : if (poCT != nullptr)
1201 : {
1202 0 : dfLatitude =
1203 0 : adfSrcGeoTransform[3] + adfSrcGeoTransform[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 = adfSrcGeoTransform[0] + adfSrcGeoTransform[1] * nXSize;
1396 33 : CPL_LSBPTR64(&dfLongitude);
1397 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfLongitude, 8, 1, fp));
1398 33 : dfLongitude = adfSrcGeoTransform[0];
1399 33 : CPL_LSBPTR64(&dfLongitude);
1400 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfLongitude, 8, 1, fp));
1401 33 : dfLatitude = adfSrcGeoTransform[3];
1402 33 : CPL_LSBPTR64(&dfLatitude);
1403 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfLatitude, 8, 1, fp));
1404 33 : dfLatitude = adfSrcGeoTransform[3] + adfSrcGeoTransform[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 = adfSrcGeoTransform[0] + adfSrcGeoTransform[1] * nXSize;
1477 33 : CPL_LSBPTR64(&dfTemp);
1478 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfTemp, 8, 1, fp));
1479 : // Min x.
1480 33 : dfTemp = adfSrcGeoTransform[0];
1481 33 : CPL_LSBPTR64(&dfTemp);
1482 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfTemp, 8, 1, fp));
1483 : // Max y.
1484 33 : dfTemp = adfSrcGeoTransform[3];
1485 33 : CPL_LSBPTR64(&dfTemp);
1486 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfTemp, 8, 1, fp));
1487 : // Min y.
1488 33 : dfTemp = adfSrcGeoTransform[3] + adfSrcGeoTransform[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 = adfSrcGeoTransform[1];
1499 33 : CPL_LSBPTR64(&dfTemp);
1500 33 : CPL_IGNORE_RET_VAL(VSIFWriteL(&dfTemp, 8, 1, fp));
1501 : // Y resolution.
1502 33 : dfTemp = fabs(adfSrcGeoTransform[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 0 : char *const pszDirname = CPLStrdup(CPLGetPath(pszFilename));
1607 0 : char *const pszBasename = CPLStrdup(CPLGetBasename(pszFilename));
1608 : char *pszPrjFilename =
1609 0 : CPLStrdup(CPLFormFilename(pszDirname, pszBasename, "prj"));
1610 0 : fp = VSIFOpenL(pszPrjFilename, "wt");
1611 0 : if (fp != nullptr)
1612 : {
1613 0 : CPL_IGNORE_RET_VAL(VSIFWriteL(pszESRIProjection, 1,
1614 : strlen(pszESRIProjection), fp));
1615 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1616 : }
1617 : else
1618 : {
1619 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to create file %s.",
1620 : pszPrjFilename);
1621 : }
1622 0 : CPLFree(pszDirname);
1623 0 : CPLFree(pszBasename);
1624 0 : CPLFree(pszPrjFilename);
1625 : }
1626 0 : CPLFree(pszESRIProjection);
1627 : }
1628 33 : return static_cast<GDALDataset *>(GDALOpen(pszFilename, GA_ReadOnly));
1629 : }
1630 :
1631 : /************************************************************************/
1632 : /* GDALRegister_LCP() */
1633 : /************************************************************************/
1634 :
1635 1595 : void GDALRegister_LCP()
1636 :
1637 : {
1638 1595 : if (GDALGetDriverByName("LCP") != nullptr)
1639 302 : return;
1640 :
1641 1293 : GDALDriver *poDriver = new GDALDriver();
1642 :
1643 1293 : poDriver->SetDescription("LCP");
1644 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1645 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1646 1293 : "FARSITE v.4 Landscape File (.lcp)");
1647 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "lcp");
1648 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/lcp.html");
1649 :
1650 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1651 :
1652 1293 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Int16");
1653 1293 : poDriver->SetMetadataItem(
1654 : GDAL_DMD_CREATIONOPTIONLIST,
1655 : "<CreationOptionList>"
1656 : " <Option name='ELEVATION_UNIT' type='string-select' "
1657 : "default='METERS' description='Elevation units'>"
1658 : " <Value>METERS</Value>"
1659 : " <Value>FEET</Value>"
1660 : " </Option>"
1661 : " <Option name='SLOPE_UNIT' type='string-select' default='DEGREES' "
1662 : "description='Slope units'>"
1663 : " <Value>DEGREES</Value>"
1664 : " <Value>PERCENT</Value>"
1665 : " </Option>"
1666 : " <Option name='ASPECT_UNIT' type='string-select' "
1667 : "default='AZIMUTH_DEGREES'>"
1668 : " <Value>GRASS_CATEGORIES</Value>"
1669 : " <Value>AZIMUTH_DEGREES</Value>"
1670 : " <Value>GRASS_DEGREES</Value>"
1671 : " </Option>"
1672 : " <Option name='FUEL_MODEL_OPTION' type='string-select' "
1673 : "default='NO_CUSTOM_AND_NO_FILE'>"
1674 : " <Value>NO_CUSTOM_AND_NO_FILE</Value>"
1675 : " <Value>CUSTOM_AND_NO_FILE</Value>"
1676 : " <Value>NO_CUSTOM_AND_FILE</Value>"
1677 : " <Value>CUSTOM_AND_FILE</Value>"
1678 : " </Option>"
1679 : " <Option name='CANOPY_COV_UNIT' type='string-select' "
1680 : "default='PERCENT'>"
1681 : " <Value>CATEGORIES</Value>"
1682 : " <Value>PERCENT</Value>"
1683 : " </Option>"
1684 : " <Option name='CANOPY_HT_UNIT' type='string-select' "
1685 : "default='METERS_X_10'>"
1686 : " <Value>METERS</Value>"
1687 : " <Value>FEET</Value>"
1688 : " <Value>METERS_X_10</Value>"
1689 : " <Value>FEET_X_10</Value>"
1690 : " </Option>"
1691 : " <Option name='CBH_UNIT' type='string-select' default='METERS_X_10'>"
1692 : " <Value>METERS</Value>"
1693 : " <Value>FEET</Value>"
1694 : " <Value>METERS_X_10</Value>"
1695 : " <Value>FEET_X_10</Value>"
1696 : " </Option>"
1697 : " <Option name='CBD_UNIT' type='string-select' "
1698 : "default='KG_PER_CUBIC_METER_X_100'>"
1699 : " <Value>KG_PER_CUBIC_METER</Value>"
1700 : " <Value>POUND_PER_CUBIC_FOOT</Value>"
1701 : " <Value>KG_PER_CUBIC_METER_X_100</Value>"
1702 : " <Value>POUND_PER_CUBIC_FOOT_X_1000</Value>"
1703 : " </Option>"
1704 : " <Option name='DUFF_UNIT' type='string-select' "
1705 : "default='MG_PER_HECTARE_X_10'>"
1706 : " <Value>MG_PER_HECTARE_X_10</Value>"
1707 : " <Value>TONS_PER_ACRE_X_10</Value>"
1708 : " </Option>"
1709 : // Kyle does not think we need to override this, but maybe?
1710 : // " <Option name='CWD_OPTION' type='boolean' default='FALSE'
1711 : // description='Override logic for setting the coarse woody presence'/>"
1712 : // */
1713 : " <Option name='CALCULATE_STATS' type='boolean' default='YES' "
1714 : "description='Write the stats to the lcp'/>"
1715 : " <Option name='CLASSIFY_DATA' type='boolean' default='YES' "
1716 : "description='Write the stats to the lcp'/>"
1717 : " <Option name='LINEAR_UNIT' type='string-select' "
1718 : "default='SET_FROM_SRS' description='Set the linear units in the lcp'>"
1719 : " <Value>SET_FROM_SRS</Value>"
1720 : " <Value>METER</Value>"
1721 : " <Value>FOOT</Value>"
1722 : " <Value>KILOMETER</Value>"
1723 : " </Option>"
1724 : " <Option name='LATITUDE' type='int' default='' description='Set the "
1725 : "latitude for the dataset, this overrides the driver trying to set it "
1726 : "programmatically in EPSG:4269'/>"
1727 : " <Option name='DESCRIPTION' type='string' default='LCP file created "
1728 : "by GDAL' description='A short description of the lcp file'/>"
1729 1293 : "</CreationOptionList>");
1730 :
1731 1293 : poDriver->pfnOpen = LCPDataset::Open;
1732 1293 : poDriver->pfnCreateCopy = LCPDataset::CreateCopy;
1733 1293 : poDriver->pfnIdentify = LCPDataset::Identify;
1734 :
1735 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1736 : }
|
