Coverage Report

Created: 2025-06-09 08:44

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