Coverage Report

Created: 2025-06-09 07:07

/src/gdal/frmts/ers/ersdataset.cpp
Line
Count
Source (jump to first uncovered line)
1
/******************************************************************************
2
 *
3
 * Project:  ERMapper .ers Driver
4
 * Purpose:  Implementation of .ers driver.
5
 * Author:   Frank Warmerdam, warmerdam@pobox.com
6
 *
7
 ******************************************************************************
8
 * Copyright (c) 2007, Frank Warmerdam <warmerdam@pobox.com>
9
 * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10
 *
11
 * SPDX-License-Identifier: MIT
12
 ****************************************************************************/
13
14
#include "cpl_string.h"
15
#include "ershdrnode.h"
16
#include "gdal_frmts.h"
17
#include "gdal_proxy.h"
18
#include "ogr_spatialref.h"
19
#include "rawdataset.h"
20
21
#include <limits>
22
23
/************************************************************************/
24
/* ==================================================================== */
25
/*                              ERSDataset                              */
26
/* ==================================================================== */
27
/************************************************************************/
28
29
class ERSRasterBand;
30
31
class ERSDataset final : public RawDataset
32
{
33
    friend class ERSRasterBand;
34
35
    VSILFILE *fpImage;  // Image data file.
36
    GDALDataset *poDepFile;
37
38
    int bGotTransform;
39
    double adfGeoTransform[6];
40
    OGRSpatialReference m_oSRS{};
41
42
    CPLString osRawFilename;
43
44
    int bHDRDirty;
45
    ERSHdrNode *poHeader;
46
47
    const char *Find(const char *, const char *);
48
49
    int nGCPCount;
50
    GDAL_GCP *pasGCPList;
51
    OGRSpatialReference m_oGCPSRS{};
52
53
    void ReadGCPs();
54
55
    int bHasNoDataValue;
56
    double dfNoDataValue;
57
58
    CPLString osProj, osProjForced;
59
    CPLString osDatum, osDatumForced;
60
    CPLString osUnits, osUnitsForced;
61
    void WriteProjectionInfo(const char *pszProj, const char *pszDatum,
62
                             const char *pszUnits);
63
64
    CPLStringList oERSMetadataList;
65
66
  protected:
67
    int CloseDependentDatasets() override;
68
69
    CPLErr Close() override;
70
71
  public:
72
    ERSDataset();
73
    ~ERSDataset() override;
74
75
    CPLErr FlushCache(bool bAtClosing) override;
76
    CPLErr GetGeoTransform(double *padfTransform) override;
77
    CPLErr SetGeoTransform(double *padfTransform) override;
78
    const OGRSpatialReference *GetSpatialRef() const override;
79
    CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
80
    char **GetFileList(void) override;
81
82
    int GetGCPCount() override;
83
    const OGRSpatialReference *GetGCPSpatialRef() const override;
84
    const GDAL_GCP *GetGCPs() override;
85
    CPLErr SetGCPs(int nGCPCountIn, const GDAL_GCP *pasGCPListIn,
86
                   const OGRSpatialReference *poSRS) override;
87
88
    char **GetMetadataDomainList() override;
89
    const char *GetMetadataItem(const char *pszName,
90
                                const char *pszDomain = "") override;
91
    char **GetMetadata(const char *pszDomain = "") override;
92
93
    static GDALDataset *Open(GDALOpenInfo *);
94
    static int Identify(GDALOpenInfo *);
95
    static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
96
                               int nBandsIn, GDALDataType eType,
97
                               char **papszParamList);
98
};
99
100
/************************************************************************/
101
/*                            ERSDataset()                             */
102
/************************************************************************/
103
104
ERSDataset::ERSDataset()
105
4
    : fpImage(nullptr), poDepFile(nullptr), bGotTransform(FALSE),
106
4
      bHDRDirty(FALSE), poHeader(nullptr), nGCPCount(0), pasGCPList(nullptr),
107
4
      bHasNoDataValue(FALSE), dfNoDataValue(0.0)
108
4
{
109
4
    m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
110
4
    m_oGCPSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
111
4
    adfGeoTransform[0] = 0.0;
112
4
    adfGeoTransform[1] = 1.0;
113
4
    adfGeoTransform[2] = 0.0;
114
4
    adfGeoTransform[3] = 0.0;
115
4
    adfGeoTransform[4] = 0.0;
116
4
    adfGeoTransform[5] = 1.0;
117
4
}
118
119
/************************************************************************/
120
/*                            ~ERSDataset()                            */
121
/************************************************************************/
122
123
ERSDataset::~ERSDataset()
124
125
4
{
126
4
    ERSDataset::Close();
127
4
}
128
129
/************************************************************************/
130
/*                              Close()                                 */
131
/************************************************************************/
132
133
CPLErr ERSDataset::Close()
134
4
{
135
4
    CPLErr eErr = CE_None;
136
4
    if (nOpenFlags != OPEN_FLAGS_CLOSED)
137
4
    {
138
4
        if (ERSDataset::FlushCache(true) != CE_None)
139
0
            eErr = CE_Failure;
140
141
4
        if (fpImage != nullptr)
142
3
        {
143
3
            VSIFCloseL(fpImage);
144
3
        }
145
146
4
        ERSDataset::CloseDependentDatasets();
147
148
4
        if (nGCPCount > 0)
149
0
        {
150
0
            GDALDeinitGCPs(nGCPCount, pasGCPList);
151
0
            CPLFree(pasGCPList);
152
0
        }
153
154
4
        if (poHeader != nullptr)
155
4
            delete poHeader;
156
157
4
        if (GDALPamDataset::Close() != CE_None)
158
0
            eErr = CE_Failure;
159
4
    }
160
4
    return eErr;
161
4
}
162
163
/************************************************************************/
164
/*                      CloseDependentDatasets()                        */
165
/************************************************************************/
166
167
int ERSDataset::CloseDependentDatasets()
168
4
{
169
4
    int bHasDroppedRef = GDALPamDataset::CloseDependentDatasets();
170
171
4
    if (poDepFile != nullptr)
172
0
    {
173
0
        bHasDroppedRef = TRUE;
174
175
0
        for (int iBand = 0; iBand < nBands; iBand++)
176
0
        {
177
0
            delete papoBands[iBand];
178
0
            papoBands[iBand] = nullptr;
179
0
        }
180
0
        nBands = 0;
181
182
0
        GDALClose((GDALDatasetH)poDepFile);
183
0
        poDepFile = nullptr;
184
0
    }
185
186
4
    return bHasDroppedRef;
187
4
}
188
189
/************************************************************************/
190
/*                             FlushCache()                             */
191
/************************************************************************/
192
193
CPLErr ERSDataset::FlushCache(bool bAtClosing)
194
195
4
{
196
4
    CPLErr eErr = CE_None;
197
4
    if (bHDRDirty)
198
0
    {
199
0
        VSILFILE *fpERS = VSIFOpenL(GetDescription(), "w");
200
0
        if (fpERS == nullptr)
201
0
        {
202
0
            eErr = CE_Failure;
203
0
            CPLError(CE_Failure, CPLE_OpenFailed,
204
0
                     "Unable to rewrite %s header.", GetDescription());
205
0
        }
206
0
        else
207
0
        {
208
0
            if (VSIFPrintfL(fpERS, "DatasetHeader Begin\n") <= 0)
209
0
                eErr = CE_Failure;
210
0
            poHeader->WriteSelf(fpERS, 1);
211
0
            if (VSIFPrintfL(fpERS, "DatasetHeader End\n") <= 0)
212
0
                eErr = CE_Failure;
213
0
            if (VSIFCloseL(fpERS) != 0)
214
0
                eErr = CE_Failure;
215
0
        }
216
0
    }
217
218
4
    if (RawDataset::FlushCache(bAtClosing) != CE_None)
219
0
        eErr = CE_Failure;
220
4
    return eErr;
221
4
}
222
223
/************************************************************************/
224
/*                      GetMetadataDomainList()                         */
225
/************************************************************************/
226
227
char **ERSDataset::GetMetadataDomainList()
228
0
{
229
0
    return BuildMetadataDomainList(GDALPamDataset::GetMetadataDomainList(),
230
0
                                   TRUE, "ERS", nullptr);
231
0
}
232
233
/************************************************************************/
234
/*                           GetMetadataItem()                          */
235
/************************************************************************/
236
237
const char *ERSDataset::GetMetadataItem(const char *pszName,
238
                                        const char *pszDomain)
239
0
{
240
0
    if (pszDomain != nullptr && EQUAL(pszDomain, "ERS") && pszName != nullptr)
241
0
    {
242
0
        if (EQUAL(pszName, "PROJ"))
243
0
            return osProj.size() ? osProj.c_str() : nullptr;
244
0
        if (EQUAL(pszName, "DATUM"))
245
0
            return osDatum.size() ? osDatum.c_str() : nullptr;
246
0
        if (EQUAL(pszName, "UNITS"))
247
0
            return osUnits.size() ? osUnits.c_str() : nullptr;
248
0
    }
249
0
    return GDALPamDataset::GetMetadataItem(pszName, pszDomain);
250
0
}
251
252
/************************************************************************/
253
/*                            GetMetadata()                             */
254
/************************************************************************/
255
256
char **ERSDataset::GetMetadata(const char *pszDomain)
257
258
0
{
259
0
    if (pszDomain != nullptr && EQUAL(pszDomain, "ERS"))
260
0
    {
261
0
        oERSMetadataList.Clear();
262
0
        if (!osProj.empty())
263
0
            oERSMetadataList.AddString(
264
0
                CPLSPrintf("%s=%s", "PROJ", osProj.c_str()));
265
0
        if (!osDatum.empty())
266
0
            oERSMetadataList.AddString(
267
0
                CPLSPrintf("%s=%s", "DATUM", osDatum.c_str()));
268
0
        if (!osUnits.empty())
269
0
            oERSMetadataList.AddString(
270
0
                CPLSPrintf("%s=%s", "UNITS", osUnits.c_str()));
271
0
        return oERSMetadataList.List();
272
0
    }
273
274
0
    return GDALPamDataset::GetMetadata(pszDomain);
275
0
}
276
277
/************************************************************************/
278
/*                            GetGCPCount()                             */
279
/************************************************************************/
280
281
int ERSDataset::GetGCPCount()
282
283
0
{
284
0
    return nGCPCount;
285
0
}
286
287
/************************************************************************/
288
/*                          GetGCPSpatialRef()                          */
289
/************************************************************************/
290
291
const OGRSpatialReference *ERSDataset::GetGCPSpatialRef() const
292
293
0
{
294
0
    return m_oGCPSRS.IsEmpty() ? nullptr : &m_oGCPSRS;
295
0
}
296
297
/************************************************************************/
298
/*                               GetGCPs()                              */
299
/************************************************************************/
300
301
const GDAL_GCP *ERSDataset::GetGCPs()
302
303
0
{
304
0
    return pasGCPList;
305
0
}
306
307
/************************************************************************/
308
/*                              SetGCPs()                               */
309
/************************************************************************/
310
311
CPLErr ERSDataset::SetGCPs(int nGCPCountIn, const GDAL_GCP *pasGCPListIn,
312
                           const OGRSpatialReference *poSRS)
313
314
0
{
315
    /* -------------------------------------------------------------------- */
316
    /*      Clean old gcps.                                                 */
317
    /* -------------------------------------------------------------------- */
318
0
    m_oGCPSRS.Clear();
319
320
0
    if (nGCPCount > 0)
321
0
    {
322
0
        GDALDeinitGCPs(nGCPCount, pasGCPList);
323
0
        CPLFree(pasGCPList);
324
325
0
        pasGCPList = nullptr;
326
0
        nGCPCount = 0;
327
0
    }
328
329
    /* -------------------------------------------------------------------- */
330
    /*      Copy new ones.                                                  */
331
    /* -------------------------------------------------------------------- */
332
0
    nGCPCount = nGCPCountIn;
333
0
    pasGCPList = GDALDuplicateGCPs(nGCPCount, pasGCPListIn);
334
0
    if (poSRS)
335
0
        m_oGCPSRS = *poSRS;
336
337
    /* -------------------------------------------------------------------- */
338
    /*      Setup the header contents corresponding to these GCPs.          */
339
    /* -------------------------------------------------------------------- */
340
0
    bHDRDirty = TRUE;
341
342
0
    poHeader->Set("RasterInfo.WarpControl.WarpType", "Polynomial");
343
0
    if (nGCPCount > 6)
344
0
        poHeader->Set("RasterInfo.WarpControl.WarpOrder", "2");
345
0
    else
346
0
        poHeader->Set("RasterInfo.WarpControl.WarpOrder", "1");
347
0
    poHeader->Set("RasterInfo.WarpControl.WarpSampling", "Nearest");
348
349
    /* -------------------------------------------------------------------- */
350
    /*      Translate the projection.                                       */
351
    /* -------------------------------------------------------------------- */
352
0
    char szERSProj[32], szERSDatum[32], szERSUnits[32];
353
354
0
    m_oGCPSRS.exportToERM(szERSProj, szERSDatum, szERSUnits);
355
356
    /* Write the above computed values, unless they have been overridden by */
357
    /* the creation options PROJ, DATUM or UNITS */
358
359
0
    poHeader->Set("RasterInfo.WarpControl.CoordinateSpace.Datum",
360
0
                  CPLString().Printf("\"%s\"", (osDatum.size())
361
0
                                                   ? osDatum.c_str()
362
0
                                                   : szERSDatum));
363
0
    poHeader->Set("RasterInfo.WarpControl.CoordinateSpace.Projection",
364
0
                  CPLString().Printf("\"%s\"", (osProj.size()) ? osProj.c_str()
365
0
                                                               : szERSProj));
366
0
    poHeader->Set("RasterInfo.WarpControl.CoordinateSpace.CoordinateType",
367
0
                  CPLString().Printf("EN"));
368
0
    poHeader->Set("RasterInfo.WarpControl.CoordinateSpace.Units",
369
0
                  CPLString().Printf("\"%s\"", (osUnits.size())
370
0
                                                   ? osUnits.c_str()
371
0
                                                   : szERSUnits));
372
0
    poHeader->Set("RasterInfo.WarpControl.CoordinateSpace.Rotation", "0:0:0.0");
373
374
    /* -------------------------------------------------------------------- */
375
    /*      Translate the GCPs.                                             */
376
    /* -------------------------------------------------------------------- */
377
0
    CPLString osControlPoints = "{\n";
378
379
0
    for (int iGCP = 0; iGCP < nGCPCount; iGCP++)
380
0
    {
381
0
        CPLString osLine;
382
383
0
        CPLString osId = pasGCPList[iGCP].pszId;
384
0
        if (osId.empty())
385
0
            osId.Printf("%d", iGCP + 1);
386
387
0
        osLine.Printf(
388
0
            "\t\t\t\t\"%s\"\tYes\tYes\t%.6f\t%.6f\t%.15g\t%.15g\t%.15g\n",
389
0
            osId.c_str(), pasGCPList[iGCP].dfGCPPixel,
390
0
            pasGCPList[iGCP].dfGCPLine, pasGCPList[iGCP].dfGCPX,
391
0
            pasGCPList[iGCP].dfGCPY, pasGCPList[iGCP].dfGCPZ);
392
0
        osControlPoints += osLine;
393
0
    }
394
0
    osControlPoints += "\t\t}";
395
396
0
    poHeader->Set("RasterInfo.WarpControl.ControlPoints", osControlPoints);
397
398
0
    return CE_None;
399
0
}
400
401
/************************************************************************/
402
/*                          GetSpatialRef()                             */
403
/************************************************************************/
404
405
const OGRSpatialReference *ERSDataset::GetSpatialRef() const
406
407
0
{
408
    // try xml first
409
0
    const auto poSRS = GDALPamDataset::GetSpatialRef();
410
0
    if (poSRS)
411
0
        return poSRS;
412
413
0
    return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
414
0
}
415
416
/************************************************************************/
417
/*                           SetSpatialRef()                            */
418
/************************************************************************/
419
420
CPLErr ERSDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
421
422
0
{
423
0
    if (poSRS == nullptr && m_oSRS.IsEmpty())
424
0
        return CE_None;
425
0
    if (poSRS != nullptr && poSRS->IsSame(&m_oSRS))
426
0
        return CE_None;
427
428
0
    m_oSRS.Clear();
429
0
    if (poSRS)
430
0
        m_oSRS = *poSRS;
431
432
0
    char szERSProj[32], szERSDatum[32], szERSUnits[32];
433
434
0
    m_oSRS.exportToERM(szERSProj, szERSDatum, szERSUnits);
435
436
    /* Write the above computed values, unless they have been overridden by */
437
    /* the creation options PROJ, DATUM or UNITS */
438
0
    if (!osProjForced.empty())
439
0
        osProj = osProjForced;
440
0
    else
441
0
        osProj = szERSProj;
442
0
    if (!osDatumForced.empty())
443
0
        osDatum = osDatumForced;
444
0
    else
445
0
        osDatum = szERSDatum;
446
0
    if (!osUnitsForced.empty())
447
0
        osUnits = osUnitsForced;
448
0
    else
449
0
        osUnits = szERSUnits;
450
451
0
    WriteProjectionInfo(osProj, osDatum, osUnits);
452
453
0
    return CE_None;
454
0
}
455
456
/************************************************************************/
457
/*                         WriteProjectionInfo()                        */
458
/************************************************************************/
459
460
void ERSDataset::WriteProjectionInfo(const char *pszProj, const char *pszDatum,
461
                                     const char *pszUnits)
462
0
{
463
0
    bHDRDirty = TRUE;
464
0
    poHeader->Set("CoordinateSpace.Datum",
465
0
                  CPLString().Printf("\"%s\"", pszDatum));
466
0
    poHeader->Set("CoordinateSpace.Projection",
467
0
                  CPLString().Printf("\"%s\"", pszProj));
468
0
    poHeader->Set("CoordinateSpace.CoordinateType", CPLString().Printf("EN"));
469
0
    poHeader->Set("CoordinateSpace.Units",
470
0
                  CPLString().Printf("\"%s\"", pszUnits));
471
0
    poHeader->Set("CoordinateSpace.Rotation", "0:0:0.0");
472
473
    /* -------------------------------------------------------------------- */
474
    /*      It seems that CoordinateSpace needs to come before              */
475
    /*      RasterInfo.  Try moving it up manually.                         */
476
    /* -------------------------------------------------------------------- */
477
0
    int iRasterInfo = -1;
478
0
    int iCoordSpace = -1;
479
480
0
    for (int i = 0; i < poHeader->nItemCount; i++)
481
0
    {
482
0
        if (EQUAL(poHeader->papszItemName[i], "RasterInfo"))
483
0
            iRasterInfo = i;
484
485
0
        if (EQUAL(poHeader->papszItemName[i], "CoordinateSpace"))
486
0
        {
487
0
            iCoordSpace = i;
488
0
            break;
489
0
        }
490
0
    }
491
492
0
    if (iCoordSpace > iRasterInfo && iRasterInfo != -1)
493
0
    {
494
0
        for (int i = iCoordSpace; i > 0 && i != iRasterInfo; i--)
495
0
        {
496
497
0
            ERSHdrNode *poTemp = poHeader->papoItemChild[i];
498
0
            poHeader->papoItemChild[i] = poHeader->papoItemChild[i - 1];
499
0
            poHeader->papoItemChild[i - 1] = poTemp;
500
501
0
            char *pszTemp = poHeader->papszItemName[i];
502
0
            poHeader->papszItemName[i] = poHeader->papszItemName[i - 1];
503
0
            poHeader->papszItemName[i - 1] = pszTemp;
504
505
0
            pszTemp = poHeader->papszItemValue[i];
506
0
            poHeader->papszItemValue[i] = poHeader->papszItemValue[i - 1];
507
0
            poHeader->papszItemValue[i - 1] = pszTemp;
508
0
        }
509
0
    }
510
0
}
511
512
/************************************************************************/
513
/*                          GetGeoTransform()                           */
514
/************************************************************************/
515
516
CPLErr ERSDataset::GetGeoTransform(double *padfTransform)
517
518
0
{
519
0
    if (bGotTransform)
520
0
    {
521
0
        memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
522
0
        return CE_None;
523
0
    }
524
525
0
    return GDALPamDataset::GetGeoTransform(padfTransform);
526
0
}
527
528
/************************************************************************/
529
/*                          SetGeoTransform()                           */
530
/************************************************************************/
531
532
CPLErr ERSDataset::SetGeoTransform(double *padfTransform)
533
534
0
{
535
0
    if (memcmp(padfTransform, adfGeoTransform, sizeof(double) * 6) == 0)
536
0
        return CE_None;
537
538
0
    if (adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0)
539
0
    {
540
0
        CPLError(CE_Failure, CPLE_AppDefined,
541
0
                 "Rotated and skewed geotransforms not currently supported for "
542
0
                 "ERS driver.");
543
0
        return CE_Failure;
544
0
    }
545
546
0
    bGotTransform = TRUE;
547
0
    memcpy(adfGeoTransform, padfTransform, sizeof(double) * 6);
548
549
0
    bHDRDirty = TRUE;
550
551
0
    poHeader->Set("RasterInfo.CellInfo.Xdimension",
552
0
                  CPLString().Printf("%.15g", fabs(adfGeoTransform[1])));
553
0
    poHeader->Set("RasterInfo.CellInfo.Ydimension",
554
0
                  CPLString().Printf("%.15g", fabs(adfGeoTransform[5])));
555
0
    poHeader->Set("RasterInfo.RegistrationCoord.Eastings",
556
0
                  CPLString().Printf("%.15g", adfGeoTransform[0]));
557
0
    poHeader->Set("RasterInfo.RegistrationCoord.Northings",
558
0
                  CPLString().Printf("%.15g", adfGeoTransform[3]));
559
560
0
    if (CPLAtof(poHeader->Find("RasterInfo.RegistrationCellX", "0")) != 0.0 ||
561
0
        CPLAtof(poHeader->Find("RasterInfo.RegistrationCellY", "0")) != 0.0)
562
0
    {
563
        // Reset RegistrationCellX/Y to 0 if the header gets rewritten (#5493)
564
0
        poHeader->Set("RasterInfo.RegistrationCellX", "0");
565
0
        poHeader->Set("RasterInfo.RegistrationCellY", "0");
566
0
    }
567
568
0
    return CE_None;
569
0
}
570
571
/************************************************************************/
572
/*                             ERSDMS2Dec()                             */
573
/*                                                                      */
574
/*      Convert ERS DMS format to decimal degrees.   Input is like      */
575
/*      "-180:00:00".                                                   */
576
/************************************************************************/
577
578
static double ERSDMS2Dec(const char *pszDMS)
579
580
0
{
581
0
    char **papszTokens = CSLTokenizeStringComplex(pszDMS, ":", FALSE, FALSE);
582
583
0
    if (CSLCount(papszTokens) != 3)
584
0
    {
585
0
        CSLDestroy(papszTokens);
586
0
        return CPLAtof(pszDMS);
587
0
    }
588
589
0
    double dfResult = fabs(CPLAtof(papszTokens[0])) +
590
0
                      CPLAtof(papszTokens[1]) / 60.0 +
591
0
                      CPLAtof(papszTokens[2]) / 3600.0;
592
593
0
    if (CPLAtof(papszTokens[0]) < 0)
594
0
        dfResult *= -1;
595
596
0
    CSLDestroy(papszTokens);
597
0
    return dfResult;
598
0
}
599
600
/************************************************************************/
601
/*                            GetFileList()                             */
602
/************************************************************************/
603
604
char **ERSDataset::GetFileList()
605
606
0
{
607
0
    static thread_local int nRecLevel = 0;
608
0
    if (nRecLevel > 0)
609
0
        return nullptr;
610
611
    // Main data file, etc.
612
0
    char **papszFileList = GDALPamDataset::GetFileList();
613
614
    // Add raw data file if we have one.
615
0
    if (!osRawFilename.empty())
616
0
        papszFileList = CSLAddString(papszFileList, osRawFilename);
617
618
    // If we have a dependent file, merge its list of files in.
619
0
    if (poDepFile)
620
0
    {
621
0
        nRecLevel++;
622
0
        char **papszDepFiles = poDepFile->GetFileList();
623
0
        nRecLevel--;
624
0
        papszFileList = CSLInsertStrings(papszFileList, -1, papszDepFiles);
625
0
        CSLDestroy(papszDepFiles);
626
0
    }
627
628
0
    return papszFileList;
629
0
}
630
631
/************************************************************************/
632
/*                              ReadGCPs()                              */
633
/*                                                                      */
634
/*      Read the GCPs from the header.                                  */
635
/************************************************************************/
636
637
void ERSDataset::ReadGCPs()
638
639
0
{
640
0
    const char *pszCP =
641
0
        poHeader->Find("RasterInfo.WarpControl.ControlPoints", nullptr);
642
643
0
    if (pszCP == nullptr)
644
0
        return;
645
646
    /* -------------------------------------------------------------------- */
647
    /*      Parse the control points.  They will look something like:       */
648
    /*                                                                      */
649
    /*   "1035" Yes No 2344.650885 3546.419458 483270.73 3620906.21 3.105   */
650
    /* -------------------------------------------------------------------- */
651
0
    char **papszTokens = CSLTokenizeStringComplex(pszCP, "{ \t}", TRUE, FALSE);
652
0
    int nItemCount = CSLCount(papszTokens);
653
654
    /* -------------------------------------------------------------------- */
655
    /*      Work out if we have elevation values or not.                    */
656
    /* -------------------------------------------------------------------- */
657
0
    int nItemsPerLine;
658
659
0
    if (nItemCount == 7)
660
0
        nItemsPerLine = 7;
661
0
    else if (nItemCount == 8)
662
0
        nItemsPerLine = 8;
663
0
    else if (nItemCount < 14)
664
0
    {
665
0
        CPLDebug("ERS", "Invalid item count for ControlPoints");
666
0
        CSLDestroy(papszTokens);
667
0
        return;
668
0
    }
669
0
    else if (EQUAL(papszTokens[8], "Yes") || EQUAL(papszTokens[8], "No"))
670
0
        nItemsPerLine = 7;
671
0
    else if (EQUAL(papszTokens[9], "Yes") || EQUAL(papszTokens[9], "No"))
672
0
        nItemsPerLine = 8;
673
0
    else
674
0
    {
675
0
        CPLDebug("ERS", "Invalid format for ControlPoints");
676
0
        CSLDestroy(papszTokens);
677
0
        return;
678
0
    }
679
680
    /* -------------------------------------------------------------------- */
681
    /*      Setup GCPs.                                                     */
682
    /* -------------------------------------------------------------------- */
683
0
    CPLAssert(nGCPCount == 0);
684
685
0
    nGCPCount = nItemCount / nItemsPerLine;
686
0
    pasGCPList = (GDAL_GCP *)CPLCalloc(nGCPCount, sizeof(GDAL_GCP));
687
0
    GDALInitGCPs(nGCPCount, pasGCPList);
688
689
0
    for (int iGCP = 0; iGCP < nGCPCount; iGCP++)
690
0
    {
691
0
        GDAL_GCP *psGCP = pasGCPList + iGCP;
692
693
0
        CPLFree(psGCP->pszId);
694
0
        psGCP->pszId = CPLStrdup(papszTokens[iGCP * nItemsPerLine + 0]);
695
0
        psGCP->dfGCPPixel = CPLAtof(papszTokens[iGCP * nItemsPerLine + 3]);
696
0
        psGCP->dfGCPLine = CPLAtof(papszTokens[iGCP * nItemsPerLine + 4]);
697
0
        psGCP->dfGCPX = CPLAtof(papszTokens[iGCP * nItemsPerLine + 5]);
698
0
        psGCP->dfGCPY = CPLAtof(papszTokens[iGCP * nItemsPerLine + 6]);
699
0
        if (nItemsPerLine == 8)
700
0
            psGCP->dfGCPZ = CPLAtof(papszTokens[iGCP * nItemsPerLine + 7]);
701
0
    }
702
703
0
    CSLDestroy(papszTokens);
704
705
    /* -------------------------------------------------------------------- */
706
    /*      Parse the GCP projection.                                       */
707
    /* -------------------------------------------------------------------- */
708
0
    osProj =
709
0
        poHeader->Find("RasterInfo.WarpControl.CoordinateSpace.Projection", "");
710
0
    osDatum =
711
0
        poHeader->Find("RasterInfo.WarpControl.CoordinateSpace.Datum", "");
712
0
    osUnits =
713
0
        poHeader->Find("RasterInfo.WarpControl.CoordinateSpace.Units", "");
714
715
0
    m_oGCPSRS.importFromERM(!osProj.empty() ? osProj.c_str() : "RAW",
716
0
                            !osDatum.empty() ? osDatum.c_str() : "WGS84",
717
0
                            !osUnits.empty() ? osUnits.c_str() : "METERS");
718
0
}
719
720
/************************************************************************/
721
/* ==================================================================== */
722
/*                             ERSRasterBand                            */
723
/* ==================================================================== */
724
/************************************************************************/
725
726
class ERSRasterBand final : public RawRasterBand
727
{
728
  public:
729
    ERSRasterBand(GDALDataset *poDS, int nBand, VSILFILE *fpRaw,
730
                  vsi_l_offset nImgOffset, int nPixelOffset, int nLineOffset,
731
                  GDALDataType eDataType, int bNativeOrder);
732
733
    double GetNoDataValue(int *pbSuccess = nullptr) override;
734
    CPLErr SetNoDataValue(double) override;
735
};
736
737
/************************************************************************/
738
/*                           ERSRasterBand()                            */
739
/************************************************************************/
740
741
ERSRasterBand::ERSRasterBand(GDALDataset *poDSIn, int nBandIn,
742
                             VSILFILE *fpRawIn, vsi_l_offset nImgOffsetIn,
743
                             int nPixelOffsetIn, int nLineOffsetIn,
744
                             GDALDataType eDataTypeIn, int bNativeOrderIn)
745
0
    : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
746
0
                    nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
747
0
                    RawRasterBand::OwnFP::NO)
748
0
{
749
0
}
750
751
/************************************************************************/
752
/*                           GetNoDataValue()                           */
753
/************************************************************************/
754
755
double ERSRasterBand::GetNoDataValue(int *pbSuccess)
756
0
{
757
0
    ERSDataset *poGDS = cpl::down_cast<ERSDataset *>(poDS);
758
0
    if (poGDS->bHasNoDataValue)
759
0
    {
760
0
        if (pbSuccess)
761
0
            *pbSuccess = TRUE;
762
0
        return poGDS->dfNoDataValue;
763
0
    }
764
765
0
    return RawRasterBand::GetNoDataValue(pbSuccess);
766
0
}
767
768
/************************************************************************/
769
/*                           SetNoDataValue()                           */
770
/************************************************************************/
771
772
CPLErr ERSRasterBand::SetNoDataValue(double dfNoDataValue)
773
0
{
774
0
    ERSDataset *poGDS = cpl::down_cast<ERSDataset *>(poDS);
775
0
    if (!poGDS->bHasNoDataValue || poGDS->dfNoDataValue != dfNoDataValue)
776
0
    {
777
0
        poGDS->bHasNoDataValue = TRUE;
778
0
        poGDS->dfNoDataValue = dfNoDataValue;
779
780
0
        poGDS->bHDRDirty = TRUE;
781
0
        poGDS->poHeader->Set("RasterInfo.NullCellValue",
782
0
                             CPLString().Printf("%.16g", dfNoDataValue));
783
0
    }
784
0
    return CE_None;
785
0
}
786
787
/************************************************************************/
788
/*                              Identify()                              */
789
/************************************************************************/
790
791
int ERSDataset::Identify(GDALOpenInfo *poOpenInfo)
792
793
16.5k
{
794
    /* -------------------------------------------------------------------- */
795
    /*      We assume the user selects the .ers file.                       */
796
    /* -------------------------------------------------------------------- */
797
16.5k
    CPLString osHeader((const char *)poOpenInfo->pabyHeader,
798
16.5k
                       poOpenInfo->nHeaderBytes);
799
800
16.5k
    if (osHeader.ifind("Algorithm Begin") != std::string::npos)
801
0
    {
802
0
        CPLError(CE_Failure, CPLE_OpenFailed,
803
0
                 "%s appears to be an algorithm ERS file, which is not "
804
0
                 "currently supported.",
805
0
                 poOpenInfo->pszFilename);
806
0
        return FALSE;
807
0
    }
808
809
16.5k
    if (osHeader.ifind("DatasetHeader ") != std::string::npos)
810
694
        return TRUE;
811
812
15.8k
    return FALSE;
813
16.5k
}
814
815
/************************************************************************/
816
/*                         ERSProxyRasterBand                           */
817
/************************************************************************/
818
819
namespace
820
{
821
822
static int &GetRecLevel()
823
347
{
824
347
    static thread_local int nRecLevel = 0;
825
347
    return nRecLevel;
826
347
}
827
828
class ERSProxyRasterBand final : public GDALProxyRasterBand
829
{
830
  public:
831
    explicit ERSProxyRasterBand(GDALRasterBand *poUnderlyingBand)
832
0
        : m_poUnderlyingBand(poUnderlyingBand)
833
0
    {
834
0
        poUnderlyingBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
835
0
        eDataType = poUnderlyingBand->GetRasterDataType();
836
0
    }
837
838
    int GetOverviewCount() override;
839
840
  protected:
841
    GDALRasterBand *RefUnderlyingRasterBand(bool /*bForceOpen*/) const override
842
0
    {
843
0
        return m_poUnderlyingBand;
844
0
    }
845
846
  private:
847
    GDALRasterBand *m_poUnderlyingBand;
848
};
849
850
int ERSProxyRasterBand::GetOverviewCount()
851
0
{
852
0
    int &nRecLevel = GetRecLevel();
853
0
    nRecLevel++;
854
0
    int nRet = GDALProxyRasterBand::GetOverviewCount();
855
0
    nRecLevel--;
856
0
    return nRet;
857
0
}
858
859
}  // namespace
860
861
/************************************************************************/
862
/*                                Open()                                */
863
/************************************************************************/
864
865
GDALDataset *ERSDataset::Open(GDALOpenInfo *poOpenInfo)
866
867
347
{
868
347
    if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
869
0
        return nullptr;
870
871
347
    int &nRecLevel = GetRecLevel();
872
    // cppcheck-suppress knownConditionTrueFalse
873
347
    if (nRecLevel)
874
0
    {
875
0
        CPLError(CE_Failure, CPLE_AppDefined,
876
0
                 "Attempt at recursively opening ERS dataset");
877
0
        return nullptr;
878
0
    }
879
880
    /* -------------------------------------------------------------------- */
881
    /*      Ingest the file as a tree of header nodes.                      */
882
    /* -------------------------------------------------------------------- */
883
347
    ERSHdrNode *poHeader = new ERSHdrNode();
884
885
347
    if (!poHeader->ParseHeader(poOpenInfo->fpL))
886
325
    {
887
325
        delete poHeader;
888
325
        VSIFCloseL(poOpenInfo->fpL);
889
325
        poOpenInfo->fpL = nullptr;
890
325
        return nullptr;
891
325
    }
892
893
22
    VSIFCloseL(poOpenInfo->fpL);
894
22
    poOpenInfo->fpL = nullptr;
895
896
    /* -------------------------------------------------------------------- */
897
    /*      Do we have the minimum required information from this header?   */
898
    /* -------------------------------------------------------------------- */
899
22
    if (poHeader->Find("RasterInfo.NrOfLines") == nullptr ||
900
22
        poHeader->Find("RasterInfo.NrOfCellsPerLine") == nullptr ||
901
22
        poHeader->Find("RasterInfo.NrOfBands") == nullptr)
902
18
    {
903
18
        if (poHeader->FindNode("Algorithm") != nullptr)
904
0
        {
905
0
            CPLError(CE_Failure, CPLE_OpenFailed,
906
0
                     "%s appears to be an algorithm ERS file, which is not "
907
0
                     "currently supported.",
908
0
                     poOpenInfo->pszFilename);
909
0
        }
910
18
        delete poHeader;
911
18
        return nullptr;
912
18
    }
913
914
    /* -------------------------------------------------------------------- */
915
    /*      Create a corresponding GDALDataset.                             */
916
    /* -------------------------------------------------------------------- */
917
4
    auto poDS = std::make_unique<ERSDataset>();
918
4
    poDS->poHeader = poHeader;
919
4
    poDS->eAccess = poOpenInfo->eAccess;
920
921
    /* -------------------------------------------------------------------- */
922
    /*      Capture some information from the file that is of interest.     */
923
    /* -------------------------------------------------------------------- */
924
4
    int nBands = atoi(poHeader->Find("RasterInfo.NrOfBands"));
925
4
    poDS->nRasterXSize = atoi(poHeader->Find("RasterInfo.NrOfCellsPerLine"));
926
4
    poDS->nRasterYSize = atoi(poHeader->Find("RasterInfo.NrOfLines"));
927
928
4
    if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
929
4
        !GDALCheckBandCount(nBands, FALSE))
930
0
    {
931
0
        return nullptr;
932
0
    }
933
934
    /* -------------------------------------------------------------------- */
935
    /*     Get the HeaderOffset if it exists in the header                  */
936
    /* -------------------------------------------------------------------- */
937
4
    GIntBig nHeaderOffset = 0;
938
4
    const char *pszHeaderOffset = poHeader->Find("HeaderOffset");
939
4
    if (pszHeaderOffset != nullptr)
940
0
    {
941
0
        nHeaderOffset = CPLAtoGIntBig(pszHeaderOffset);
942
0
        if (nHeaderOffset < 0)
943
0
        {
944
0
            CPLError(CE_Failure, CPLE_AppDefined,
945
0
                     "Illegal value for HeaderOffset: %s", pszHeaderOffset);
946
0
            return nullptr;
947
0
        }
948
0
    }
949
950
    /* -------------------------------------------------------------------- */
951
    /*      Establish the data type.                                        */
952
    /* -------------------------------------------------------------------- */
953
4
    CPLString osCellType =
954
4
        poHeader->Find("RasterInfo.CellType", "Unsigned8BitInteger");
955
4
    GDALDataType eType;
956
4
    if (EQUAL(osCellType, "Unsigned8BitInteger"))
957
4
        eType = GDT_Byte;
958
0
    else if (EQUAL(osCellType, "Signed8BitInteger"))
959
0
        eType = GDT_Int8;
960
0
    else if (EQUAL(osCellType, "Unsigned16BitInteger"))
961
0
        eType = GDT_UInt16;
962
0
    else if (EQUAL(osCellType, "Signed16BitInteger"))
963
0
        eType = GDT_Int16;
964
0
    else if (EQUAL(osCellType, "Unsigned32BitInteger"))
965
0
        eType = GDT_UInt32;
966
0
    else if (EQUAL(osCellType, "Signed32BitInteger"))
967
0
        eType = GDT_Int32;
968
0
    else if (EQUAL(osCellType, "IEEE4ByteReal"))
969
0
        eType = GDT_Float32;
970
0
    else if (EQUAL(osCellType, "IEEE8ByteReal"))
971
0
        eType = GDT_Float64;
972
0
    else
973
0
    {
974
0
        CPLDebug("ERS", "Unknown CellType '%s'", osCellType.c_str());
975
0
        eType = GDT_Byte;
976
0
    }
977
978
    /* -------------------------------------------------------------------- */
979
    /*      Pick up the word order.                                         */
980
    /* -------------------------------------------------------------------- */
981
4
    const int bNative =
982
4
#ifdef CPL_LSB
983
4
        EQUAL(poHeader->Find("ByteOrder", "LSBFirst"), "LSBFirst")
984
#else
985
        EQUAL(poHeader->Find("ByteOrder", "MSBFirst"), "MSBFirst")
986
#endif
987
4
        ;
988
    /* -------------------------------------------------------------------- */
989
    /*      Figure out the name of the target file.                         */
990
    /* -------------------------------------------------------------------- */
991
4
    CPLString osPath = CPLGetPathSafe(poOpenInfo->pszFilename);
992
4
    CPLString osDataFile = poHeader->Find("DataFile", "");
993
994
4
    if (osDataFile.length() == 0)  // just strip off extension.
995
4
    {
996
4
        osDataFile = CPLGetFilename(poOpenInfo->pszFilename);
997
4
        osDataFile = osDataFile.substr(0, osDataFile.find_last_of('.'));
998
4
    }
999
1000
4
    CPLString osDataFilePath = CPLFormFilenameSafe(osPath, osDataFile, nullptr);
1001
1002
    /* -------------------------------------------------------------------- */
1003
    /*      DataSetType = Translated files are links to things like ecw     */
1004
    /*      files.                                                          */
1005
    /* -------------------------------------------------------------------- */
1006
4
    if (EQUAL(poHeader->Find("DataSetType", ""), "Translated"))
1007
0
    {
1008
0
        nRecLevel++;
1009
0
        poDS->poDepFile = GDALDataset::FromHandle(
1010
0
            GDALOpen(osDataFilePath, poOpenInfo->eAccess));
1011
0
        nRecLevel--;
1012
1013
0
        if (poDS->poDepFile != nullptr &&
1014
0
            poDS->poDepFile->GetRasterXSize() == poDS->GetRasterXSize() &&
1015
0
            poDS->poDepFile->GetRasterYSize() == poDS->GetRasterYSize() &&
1016
0
            poDS->poDepFile->GetRasterCount() >= nBands)
1017
0
        {
1018
0
            for (int iBand = 0; iBand < nBands; iBand++)
1019
0
            {
1020
                // Assume pixel interleaved.
1021
0
                poDS->SetBand(iBand + 1,
1022
0
                              new ERSProxyRasterBand(
1023
0
                                  poDS->poDepFile->GetRasterBand(iBand + 1)));
1024
0
            }
1025
0
        }
1026
0
        else
1027
0
        {
1028
0
            delete poDS->poDepFile;
1029
0
            poDS->poDepFile = nullptr;
1030
0
        }
1031
0
    }
1032
1033
    /* ==================================================================== */
1034
    /*      While ERStorage indicates a raw file.                           */
1035
    /* ==================================================================== */
1036
4
    else if (EQUAL(poHeader->Find("DataSetType", ""), "ERStorage"))
1037
4
    {
1038
        // Open data file.
1039
4
        if (poOpenInfo->eAccess == GA_Update)
1040
0
            poDS->fpImage = VSIFOpenL(osDataFilePath, "r+");
1041
4
        else
1042
4
            poDS->fpImage = VSIFOpenL(osDataFilePath, "r");
1043
1044
4
        poDS->osRawFilename = std::move(osDataFilePath);
1045
1046
4
        if (poDS->fpImage != nullptr && nBands > 0)
1047
3
        {
1048
3
            int iWordSize = GDALGetDataTypeSizeBytes(eType);
1049
1050
3
            const auto knIntMax = std::numeric_limits<int>::max();
1051
3
            if (nBands > knIntMax / iWordSize ||
1052
3
                poDS->nRasterXSize > knIntMax / (nBands * iWordSize))
1053
0
            {
1054
0
                CPLError(CE_Failure, CPLE_AppDefined,
1055
0
                         "int overflow: too large nBands and/or nRasterXSize");
1056
0
                return nullptr;
1057
0
            }
1058
1059
3
            if (!RAWDatasetCheckMemoryUsage(
1060
3
                    poDS->nRasterXSize, poDS->nRasterYSize, nBands, iWordSize,
1061
3
                    iWordSize, iWordSize * nBands * poDS->nRasterXSize,
1062
3
                    nHeaderOffset,
1063
3
                    static_cast<vsi_l_offset>(iWordSize) * poDS->nRasterXSize,
1064
3
                    poDS->fpImage))
1065
3
            {
1066
3
                return nullptr;
1067
3
            }
1068
0
            if (nHeaderOffset > std::numeric_limits<GIntBig>::max() -
1069
0
                                    static_cast<GIntBig>(nBands - 1) *
1070
0
                                        iWordSize * poDS->nRasterXSize)
1071
0
            {
1072
0
                CPLError(CE_Failure, CPLE_AppDefined,
1073
0
                         "int overflow: too large nHeaderOffset");
1074
0
                return nullptr;
1075
0
            }
1076
1077
0
            for (int iBand = 0; iBand < nBands; iBand++)
1078
0
            {
1079
                // Assume pixel interleaved.
1080
0
                auto poBand = std::make_unique<ERSRasterBand>(
1081
0
                    poDS.get(), iBand + 1, poDS->fpImage,
1082
0
                    nHeaderOffset + static_cast<vsi_l_offset>(iWordSize) *
1083
0
                                        iBand * poDS->nRasterXSize,
1084
0
                    iWordSize, iWordSize * nBands * poDS->nRasterXSize, eType,
1085
0
                    bNative);
1086
0
                if (!poBand->IsValid())
1087
0
                    return nullptr;
1088
0
                poDS->SetBand(iBand + 1, std::move(poBand));
1089
0
            }
1090
0
        }
1091
4
    }
1092
1093
    /* -------------------------------------------------------------------- */
1094
    /*      Otherwise we have an error!                                     */
1095
    /* -------------------------------------------------------------------- */
1096
1
    if (poDS->nBands == 0)
1097
1
    {
1098
1
        return nullptr;
1099
1
    }
1100
1101
    /* -------------------------------------------------------------------- */
1102
    /*      Look for band descriptions.                                     */
1103
    /* -------------------------------------------------------------------- */
1104
0
    ERSHdrNode *poRI = poHeader->FindNode("RasterInfo");
1105
1106
0
    for (int iChild = 0, iBand = 0;
1107
0
         poRI != nullptr && iChild < poRI->nItemCount && iBand < poDS->nBands;
1108
0
         iChild++)
1109
0
    {
1110
0
        if (poRI->papoItemChild[iChild] != nullptr &&
1111
0
            EQUAL(poRI->papszItemName[iChild], "BandId"))
1112
0
        {
1113
0
            const char *pszValue =
1114
0
                poRI->papoItemChild[iChild]->Find("Value", nullptr);
1115
1116
0
            iBand++;
1117
0
            if (pszValue)
1118
0
            {
1119
0
                CPLPushErrorHandler(CPLQuietErrorHandler);
1120
0
                poDS->GetRasterBand(iBand)->SetDescription(pszValue);
1121
0
                CPLPopErrorHandler();
1122
0
            }
1123
1124
0
            pszValue = poRI->papoItemChild[iChild]->Find("Units", nullptr);
1125
0
            if (pszValue)
1126
0
            {
1127
0
                CPLPushErrorHandler(CPLQuietErrorHandler);
1128
0
                poDS->GetRasterBand(iBand)->SetUnitType(pszValue);
1129
0
                CPLPopErrorHandler();
1130
0
            }
1131
0
        }
1132
0
    }
1133
1134
    /* -------------------------------------------------------------------- */
1135
    /*      Look for projection.                                            */
1136
    /* -------------------------------------------------------------------- */
1137
0
    poDS->osProj = poHeader->Find("CoordinateSpace.Projection", "");
1138
0
    poDS->osDatum = poHeader->Find("CoordinateSpace.Datum", "");
1139
0
    poDS->osUnits = poHeader->Find("CoordinateSpace.Units", "");
1140
1141
0
    poDS->m_oSRS.importFromERM(
1142
0
        !poDS->osProj.empty() ? poDS->osProj.c_str() : "RAW",
1143
0
        !poDS->osDatum.empty() ? poDS->osDatum.c_str() : "WGS84",
1144
0
        !poDS->osUnits.empty() ? poDS->osUnits.c_str() : "METERS");
1145
1146
    /* -------------------------------------------------------------------- */
1147
    /*      Look for the geotransform.                                      */
1148
    /* -------------------------------------------------------------------- */
1149
0
    if (poHeader->Find("RasterInfo.RegistrationCoord.Eastings", nullptr))
1150
0
    {
1151
0
        poDS->bGotTransform = TRUE;
1152
0
        poDS->adfGeoTransform[0] = CPLAtof(
1153
0
            poHeader->Find("RasterInfo.RegistrationCoord.Eastings", ""));
1154
0
        poDS->adfGeoTransform[1] =
1155
0
            CPLAtof(poHeader->Find("RasterInfo.CellInfo.Xdimension", "1.0"));
1156
0
        poDS->adfGeoTransform[2] = 0.0;
1157
0
        poDS->adfGeoTransform[3] = CPLAtof(
1158
0
            poHeader->Find("RasterInfo.RegistrationCoord.Northings", ""));
1159
0
        poDS->adfGeoTransform[4] = 0.0;
1160
0
        poDS->adfGeoTransform[5] =
1161
0
            -CPLAtof(poHeader->Find("RasterInfo.CellInfo.Ydimension", "1.0"));
1162
0
    }
1163
0
    else if (poHeader->Find("RasterInfo.RegistrationCoord.Latitude", nullptr) &&
1164
0
             poHeader->Find("RasterInfo.CellInfo.Xdimension", nullptr))
1165
0
    {
1166
0
        poDS->bGotTransform = TRUE;
1167
0
        poDS->adfGeoTransform[0] = ERSDMS2Dec(
1168
0
            poHeader->Find("RasterInfo.RegistrationCoord.Longitude", ""));
1169
0
        poDS->adfGeoTransform[1] =
1170
0
            CPLAtof(poHeader->Find("RasterInfo.CellInfo.Xdimension", ""));
1171
0
        poDS->adfGeoTransform[2] = 0.0;
1172
0
        poDS->adfGeoTransform[3] = ERSDMS2Dec(
1173
0
            poHeader->Find("RasterInfo.RegistrationCoord.Latitude", ""));
1174
0
        poDS->adfGeoTransform[4] = 0.0;
1175
0
        poDS->adfGeoTransform[5] =
1176
0
            -CPLAtof(poHeader->Find("RasterInfo.CellInfo.Ydimension", ""));
1177
0
    }
1178
1179
    /* -------------------------------------------------------------------- */
1180
    /*      Adjust if we have a registration cell.                          */
1181
    /* -------------------------------------------------------------------- */
1182
1183
    /* http://geospatial.intergraph.com/Libraries/Tech_Docs/ERDAS_ER_Mapper_Customization_Guide.sflb.ashx
1184
     */
1185
    /* Page 27 : */
1186
    /* RegistrationCellX and RegistrationCellY : The image X and Y
1187
       coordinates of the cell which corresponds to the Registration
1188
       Coordinate. Note that the RegistrationCellX and
1189
       RegistrationCellY can be fractional values. If
1190
       RegistrationCellX and RegistrationCellY are not specified,
1191
       they are assumed to be (0,0), which is the top left corner of the
1192
       image.
1193
       */
1194
0
    double dfCellX =
1195
0
        CPLAtof(poHeader->Find("RasterInfo.RegistrationCellX", "0"));
1196
0
    double dfCellY =
1197
0
        CPLAtof(poHeader->Find("RasterInfo.RegistrationCellY", "0"));
1198
1199
0
    if (poDS->bGotTransform)
1200
0
    {
1201
0
        poDS->adfGeoTransform[0] -= dfCellX * poDS->adfGeoTransform[1] +
1202
0
                                    dfCellY * poDS->adfGeoTransform[2];
1203
0
        poDS->adfGeoTransform[3] -= dfCellX * poDS->adfGeoTransform[4] +
1204
0
                                    dfCellY * poDS->adfGeoTransform[5];
1205
0
    }
1206
1207
    /* -------------------------------------------------------------------- */
1208
    /*      Check for null values.                                          */
1209
    /* -------------------------------------------------------------------- */
1210
0
    if (poHeader->Find("RasterInfo.NullCellValue", nullptr))
1211
0
    {
1212
0
        poDS->bHasNoDataValue = TRUE;
1213
0
        poDS->dfNoDataValue =
1214
0
            CPLAtofM(poHeader->Find("RasterInfo.NullCellValue"));
1215
1216
0
        if (poDS->poDepFile != nullptr)
1217
0
        {
1218
0
            CPLPushErrorHandler(CPLQuietErrorHandler);
1219
1220
0
            for (int iBand = 1; iBand <= poDS->nBands; iBand++)
1221
0
                poDS->GetRasterBand(iBand)->SetNoDataValue(poDS->dfNoDataValue);
1222
1223
0
            CPLPopErrorHandler();
1224
0
        }
1225
0
    }
1226
1227
    /* -------------------------------------------------------------------- */
1228
    /*      Do we have an "All" region?                                     */
1229
    /* -------------------------------------------------------------------- */
1230
0
    ERSHdrNode *poAll = nullptr;
1231
1232
0
    for (int iChild = 0; poRI != nullptr && iChild < poRI->nItemCount; iChild++)
1233
0
    {
1234
0
        if (poRI->papoItemChild[iChild] != nullptr &&
1235
0
            EQUAL(poRI->papszItemName[iChild], "RegionInfo"))
1236
0
        {
1237
0
            if (EQUAL(poRI->papoItemChild[iChild]->Find("RegionName", ""),
1238
0
                      "All"))
1239
0
                poAll = poRI->papoItemChild[iChild];
1240
0
        }
1241
0
    }
1242
1243
    /* -------------------------------------------------------------------- */
1244
    /*      Do we have statistics?                                          */
1245
    /* -------------------------------------------------------------------- */
1246
0
    if (poAll && poAll->FindNode("Stats"))
1247
0
    {
1248
0
        CPLPushErrorHandler(CPLQuietErrorHandler);
1249
1250
0
        for (int iBand = 1; iBand <= poDS->nBands; iBand++)
1251
0
        {
1252
0
            const char *pszValue =
1253
0
                poAll->FindElem("Stats.MinimumValue", iBand - 1);
1254
1255
0
            if (pszValue)
1256
0
                poDS->GetRasterBand(iBand)->SetMetadataItem(
1257
0
                    "STATISTICS_MINIMUM", pszValue);
1258
1259
0
            pszValue = poAll->FindElem("Stats.MaximumValue", iBand - 1);
1260
1261
0
            if (pszValue)
1262
0
                poDS->GetRasterBand(iBand)->SetMetadataItem(
1263
0
                    "STATISTICS_MAXIMUM", pszValue);
1264
1265
0
            pszValue = poAll->FindElem("Stats.MeanValue", iBand - 1);
1266
1267
0
            if (pszValue)
1268
0
                poDS->GetRasterBand(iBand)->SetMetadataItem("STATISTICS_MEAN",
1269
0
                                                            pszValue);
1270
1271
0
            pszValue = poAll->FindElem("Stats.MedianValue", iBand - 1);
1272
1273
0
            if (pszValue)
1274
0
                poDS->GetRasterBand(iBand)->SetMetadataItem("STATISTICS_MEDIAN",
1275
0
                                                            pszValue);
1276
0
        }
1277
1278
0
        CPLPopErrorHandler();
1279
0
    }
1280
1281
    /* -------------------------------------------------------------------- */
1282
    /*      Do we have GCPs.                                                */
1283
    /* -------------------------------------------------------------------- */
1284
0
    if (poHeader->FindNode("RasterInfo.WarpControl"))
1285
0
        poDS->ReadGCPs();
1286
1287
    /* -------------------------------------------------------------------- */
1288
    /*      Initialize any PAM information.                                 */
1289
    /* -------------------------------------------------------------------- */
1290
0
    poDS->SetDescription(poOpenInfo->pszFilename);
1291
0
    poDS->TryLoadXML();
1292
1293
    // if no SR in xml, try aux
1294
0
    const OGRSpatialReference *poSRS = poDS->GDALPamDataset::GetSpatialRef();
1295
0
    if (poSRS == nullptr)
1296
0
    {
1297
        // try aux
1298
0
        auto poAuxDS = std::unique_ptr<GDALDataset>(GDALFindAssociatedAuxFile(
1299
0
            poOpenInfo->pszFilename, GA_ReadOnly, poDS.get()));
1300
0
        if (poAuxDS)
1301
0
        {
1302
0
            poSRS = poAuxDS->GetSpatialRef();
1303
0
            if (poSRS)
1304
0
            {
1305
0
                poDS->m_oSRS = *poSRS;
1306
0
            }
1307
0
        }
1308
0
    }
1309
    /* -------------------------------------------------------------------- */
1310
    /*      Check for overviews.                                            */
1311
    /* -------------------------------------------------------------------- */
1312
0
    poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
1313
1314
0
    return poDS.release();
1315
1
}
1316
1317
/************************************************************************/
1318
/*                               Create()                               */
1319
/************************************************************************/
1320
1321
GDALDataset *ERSDataset::Create(const char *pszFilename, int nXSize, int nYSize,
1322
                                int nBandsIn, GDALDataType eType,
1323
                                char **papszOptions)
1324
1325
0
{
1326
    /* -------------------------------------------------------------------- */
1327
    /*      Verify settings.                                                */
1328
    /* -------------------------------------------------------------------- */
1329
0
    if (nBandsIn <= 0)
1330
0
    {
1331
0
        CPLError(CE_Failure, CPLE_NotSupported,
1332
0
                 "ERS driver does not support %d bands.\n", nBandsIn);
1333
0
        return nullptr;
1334
0
    }
1335
1336
0
    if (eType != GDT_Byte && eType != GDT_Int8 && eType != GDT_Int16 &&
1337
0
        eType != GDT_UInt16 && eType != GDT_Int32 && eType != GDT_UInt32 &&
1338
0
        eType != GDT_Float32 && eType != GDT_Float64)
1339
0
    {
1340
0
        CPLError(
1341
0
            CE_Failure, CPLE_AppDefined,
1342
0
            "The ERS driver does not supporting creating files of types %s.",
1343
0
            GDALGetDataTypeName(eType));
1344
0
        return nullptr;
1345
0
    }
1346
1347
    /* -------------------------------------------------------------------- */
1348
    /*      Work out the name we want to use for the .ers and binary        */
1349
    /*      data files.                                                     */
1350
    /* -------------------------------------------------------------------- */
1351
0
    CPLString osBinFile, osErsFile;
1352
1353
0
    if (EQUAL(CPLGetExtensionSafe(pszFilename).c_str(), "ers"))
1354
0
    {
1355
0
        osErsFile = pszFilename;
1356
0
        osBinFile = osErsFile.substr(0, osErsFile.length() - 4);
1357
0
    }
1358
0
    else
1359
0
    {
1360
0
        osBinFile = pszFilename;
1361
0
        osErsFile = osBinFile + ".ers";
1362
0
    }
1363
1364
    /* -------------------------------------------------------------------- */
1365
    /*      Work out some values we will write.                             */
1366
    /* -------------------------------------------------------------------- */
1367
0
    const char *pszCellType = "Unsigned8BitInteger";
1368
0
    CPL_IGNORE_RET_VAL(pszCellType);  // Make CSA happy
1369
1370
0
    if (eType == GDT_Byte)
1371
0
        pszCellType = "Unsigned8BitInteger";
1372
0
    else if (eType == GDT_Int8)
1373
0
        pszCellType = "Signed8BitInteger";
1374
0
    else if (eType == GDT_Int16)
1375
0
        pszCellType = "Signed16BitInteger";
1376
0
    else if (eType == GDT_UInt16)
1377
0
        pszCellType = "Unsigned16BitInteger";
1378
0
    else if (eType == GDT_Int32)
1379
0
        pszCellType = "Signed32BitInteger";
1380
0
    else if (eType == GDT_UInt32)
1381
0
        pszCellType = "Unsigned32BitInteger";
1382
0
    else if (eType == GDT_Float32)
1383
0
        pszCellType = "IEEE4ByteReal";
1384
0
    else if (eType == GDT_Float64)
1385
0
        pszCellType = "IEEE8ByteReal";
1386
0
    else
1387
0
    {
1388
0
        CPLAssert(false);
1389
0
    }
1390
1391
    /* -------------------------------------------------------------------- */
1392
    /*      Handling for signed eight bit data.                             */
1393
    /* -------------------------------------------------------------------- */
1394
0
    const char *pszPixelType = CSLFetchNameValue(papszOptions, "PIXELTYPE");
1395
0
    if (pszPixelType && EQUAL(pszPixelType, "SIGNEDBYTE") && eType == GDT_Byte)
1396
0
        pszCellType = "Signed8BitInteger";
1397
1398
    /* -------------------------------------------------------------------- */
1399
    /*      Write binary file.                                              */
1400
    /* -------------------------------------------------------------------- */
1401
0
    VSILFILE *fpBin = VSIFOpenL(osBinFile, "w");
1402
1403
0
    if (fpBin == nullptr)
1404
0
    {
1405
0
        CPLError(CE_Failure, CPLE_FileIO, "Failed to create %s:\n%s",
1406
0
                 osBinFile.c_str(), VSIStrerror(errno));
1407
0
        return nullptr;
1408
0
    }
1409
1410
0
    GUIntBig nSize =
1411
0
        nXSize * (GUIntBig)nYSize * nBandsIn * (GDALGetDataTypeSize(eType) / 8);
1412
0
    GByte byZero = 0;
1413
0
    if (VSIFSeekL(fpBin, nSize - 1, SEEK_SET) != 0 ||
1414
0
        VSIFWriteL(&byZero, 1, 1, fpBin) != 1)
1415
0
    {
1416
0
        CPLError(CE_Failure, CPLE_FileIO, "Failed to write %s:\n%s",
1417
0
                 osBinFile.c_str(), VSIStrerror(errno));
1418
0
        VSIFCloseL(fpBin);
1419
0
        return nullptr;
1420
0
    }
1421
0
    VSIFCloseL(fpBin);
1422
1423
    /* -------------------------------------------------------------------- */
1424
    /*      Try writing header file.                                        */
1425
    /* -------------------------------------------------------------------- */
1426
0
    VSILFILE *fpERS = VSIFOpenL(osErsFile, "w");
1427
1428
0
    if (fpERS == nullptr)
1429
0
    {
1430
0
        CPLError(CE_Failure, CPLE_FileIO, "Failed to create %s:\n%s",
1431
0
                 osErsFile.c_str(), VSIStrerror(errno));
1432
0
        return nullptr;
1433
0
    }
1434
1435
0
    VSIFPrintfL(fpERS, "DatasetHeader Begin\n");
1436
0
    VSIFPrintfL(fpERS, "\tVersion\t\t = \"6.0\"\n");
1437
0
    VSIFPrintfL(fpERS, "\tName\t\t= \"%s\"\n", CPLGetFilename(osErsFile));
1438
1439
    // Last updated requires timezone info which we don't necessarily get
1440
    // get from VSICTime() so perhaps it is better to omit this.
1441
    //    VSIFPrintfL( fpERS, "\tLastUpdated\t= %s",
1442
    //                 VSICTime( VSITime( NULL ) ) );
1443
1444
0
    VSIFPrintfL(fpERS, "\tDataSetType\t= ERStorage\n");
1445
0
    VSIFPrintfL(fpERS, "\tDataType\t= Raster\n");
1446
0
    VSIFPrintfL(fpERS, "\tByteOrder\t= LSBFirst\n");
1447
0
    VSIFPrintfL(fpERS, "\tRasterInfo Begin\n");
1448
0
    VSIFPrintfL(fpERS, "\t\tCellType\t= %s\n", pszCellType);
1449
0
    VSIFPrintfL(fpERS, "\t\tNrOfLines\t= %d\n", nYSize);
1450
0
    VSIFPrintfL(fpERS, "\t\tNrOfCellsPerLine\t= %d\n", nXSize);
1451
0
    VSIFPrintfL(fpERS, "\t\tNrOfBands\t= %d\n", nBandsIn);
1452
0
    VSIFPrintfL(fpERS, "\tRasterInfo End\n");
1453
0
    if (VSIFPrintfL(fpERS, "DatasetHeader End\n") < 17)
1454
0
    {
1455
0
        CPLError(CE_Failure, CPLE_FileIO, "Failed to write %s:\n%s",
1456
0
                 osErsFile.c_str(), VSIStrerror(errno));
1457
0
        return nullptr;
1458
0
    }
1459
1460
0
    VSIFCloseL(fpERS);
1461
1462
    /* -------------------------------------------------------------------- */
1463
    /*      Reopen.                                                         */
1464
    /* -------------------------------------------------------------------- */
1465
0
    GDALOpenInfo oOpenInfo(osErsFile, GA_Update);
1466
0
    ERSDataset *poDS = cpl::down_cast<ERSDataset *>(Open(&oOpenInfo));
1467
0
    if (poDS == nullptr)
1468
0
        return nullptr;
1469
1470
    /* -------------------------------------------------------------------- */
1471
    /*      Fetch DATUM, PROJ and UNITS creation option                     */
1472
    /* -------------------------------------------------------------------- */
1473
0
    const char *pszDatum = CSLFetchNameValue(papszOptions, "DATUM");
1474
0
    if (pszDatum)
1475
0
    {
1476
0
        poDS->osDatumForced = pszDatum;
1477
0
        poDS->osDatum = pszDatum;
1478
0
    }
1479
0
    const char *pszProj = CSLFetchNameValue(papszOptions, "PROJ");
1480
0
    if (pszProj)
1481
0
    {
1482
0
        poDS->osProjForced = pszProj;
1483
0
        poDS->osProj = pszProj;
1484
0
    }
1485
0
    const char *pszUnits = CSLFetchNameValue(papszOptions, "UNITS");
1486
0
    if (pszUnits)
1487
0
    {
1488
0
        poDS->osUnitsForced = pszUnits;
1489
0
        poDS->osUnits = pszUnits;
1490
0
    }
1491
1492
0
    if (pszDatum || pszProj || pszUnits)
1493
0
    {
1494
0
        poDS->WriteProjectionInfo(pszProj ? pszProj : "RAW",
1495
0
                                  pszDatum ? pszDatum : "RAW",
1496
0
                                  pszUnits ? pszUnits : "METERS");
1497
0
    }
1498
1499
0
    return poDS;
1500
0
}
1501
1502
/************************************************************************/
1503
/*                         GDALRegister_ERS()                           */
1504
/************************************************************************/
1505
1506
void GDALRegister_ERS()
1507
1508
2
{
1509
2
    if (GDALGetDriverByName("ERS") != nullptr)
1510
0
        return;
1511
1512
2
    GDALDriver *poDriver = new GDALDriver();
1513
1514
2
    poDriver->SetDescription("ERS");
1515
2
    poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1516
2
    poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "ERMapper .ers Labelled");
1517
2
    poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ers.html");
1518
2
    poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "ers");
1519
2
    poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
1520
2
                              "Byte Int8 Int16 UInt16 Int32 UInt32 "
1521
2
                              "Float32 Float64");
1522
1523
2
    poDriver->SetMetadataItem(
1524
2
        GDAL_DMD_CREATIONOPTIONLIST,
1525
2
        "<CreationOptionList>"
1526
2
        "   <Option name='PIXELTYPE' type='string' description='(deprecated, "
1527
2
        "use Int8 datatype) By setting this to SIGNEDBYTE, a new Byte file can "
1528
2
        "be forced to be written as signed byte'/>"
1529
2
        "   <Option name='PROJ' type='string' description='ERS Projection "
1530
2
        "Name'/>"
1531
2
        "   <Option name='DATUM' type='string' description='ERS Datum Name' />"
1532
2
        "   <Option name='UNITS' type='string-select' description='ERS "
1533
2
        "Projection Units'>"
1534
2
        "       <Value>METERS</Value>"
1535
2
        "       <Value>FEET</Value>"
1536
2
        "   </Option>"
1537
2
        "</CreationOptionList>");
1538
1539
2
    poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1540
1541
2
    poDriver->pfnOpen = ERSDataset::Open;
1542
2
    poDriver->pfnIdentify = ERSDataset::Identify;
1543
2
    poDriver->pfnCreate = ERSDataset::Create;
1544
1545
2
    GetGDALDriverManager()->RegisterDriver(poDriver);
1546
2
}