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