• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

OSGeo / gdal / 15885686134

25 Jun 2025 07:44PM UTC coverage: 71.084%. Remained the same
15885686134

push

github

rouault
gdal_priv.h: fix C++11 compatibility

573814 of 807237 relevant lines covered (71.08%)

250621.56 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

88.66
/apps/gdal_footprint_lib.cpp
1
/******************************************************************************
2
 *
3
 * Project:  GDAL Utilities
4
 * Purpose:  Footprint OGR shapes into a GDAL raster.
5
 * Authors:  Even Rouault, <even dot rouault at spatialys dot com>
6
 *
7
 ******************************************************************************
8
 * Copyright (c) 2023, Even Rouault <even dot rouault at spatialys dot com>
9
 *
10
 * SPDX-License-Identifier: MIT
11
 ****************************************************************************/
12

13
#include "cpl_port.h"
14
#include "gdal_utils.h"
15
#include "gdal_utils_priv.h"
16

17
#include <cmath>
18
#include <cstdio>
19
#include <cstdlib>
20
#include <cstring>
21
#include <algorithm>
22
#include <limits>
23
#include <vector>
24

25
#include "commonutils.h"
26
#include "cpl_conv.h"
27
#include "cpl_error.h"
28
#include "cpl_progress.h"
29
#include "cpl_string.h"
30
#include "cpl_vsi.h"
31
#include "gdal.h"
32
#include "gdal_alg.h"
33
#include "gdal_priv.h"
34
#include "ogr_api.h"
35
#include "ogr_core.h"
36
#include "memdataset.h"
37
#include "ogrsf_frmts.h"
38
#include "ogr_spatialref.h"
39
#include "gdalargumentparser.h"
40

41
constexpr const char *DEFAULT_LAYER_NAME = "footprint";
42

43
/************************************************************************/
44
/*                          GDALFootprintOptions                        */
45
/************************************************************************/
46

47
struct GDALFootprintOptions
48
{
49
    /*! output format. Use the short format name. */
50
    std::string osFormat{};
51

52
    /*! the progress function to use */
53
    GDALProgressFunc pfnProgress = GDALDummyProgress;
54

55
    /*! pointer to the progress data variable */
56
    void *pProgressData = nullptr;
57

58
    bool bCreateOutput = false;
59

60
    std::string osDestLayerName{};
61

62
    /*! Layer creation options */
63
    CPLStringList aosLCO{};
64

65
    /*! Dataset creation options */
66
    CPLStringList aosDSCO{};
67

68
    /*! Overview index: 0 = first overview level */
69
    int nOvrIndex = -1;
70

71
    /** Whether output geometry should be in georeferenced coordinates, if
72
     * possible (if explicitly requested, bOutCSGeorefRequested is also set)
73
     * false = in pixel coordinates
74
     */
75
    bool bOutCSGeoref = true;
76

77
    /** Whether -t_cs georef has been explicitly set */
78
    bool bOutCSGeorefRequested = false;
79

80
    OGRSpatialReference oOutputSRS{};
81

82
    bool bSplitPolys = false;
83

84
    double dfDensifyDistance = 0;
85

86
    double dfSimplifyTolerance = 0;
87

88
    bool bConvexHull = false;
89

90
    double dfMinRingArea = 0;
91

92
    int nMaxPoints = 100;
93

94
    /*! Source bands to take into account */
95
    std::vector<int> anBands{};
96

97
    /*! Whether to combine bands unioning (true) or intersecting (false) */
98
    bool bCombineBandsUnion = true;
99

100
    /*! Field name where to write the path of the raster. Empty if not desired */
101
    std::string osLocationFieldName = "location";
102

103
    /*! Clears the osLocationFieldName var when set */
104
    bool bClearLocation = false;
105

106
    /*! Whether to force writing absolute paths in location field. */
107
    bool bAbsolutePath = false;
108

109
    std::string osSrcNoData{};
110
};
111

112
static std::unique_ptr<GDALArgumentParser> GDALFootprintAppOptionsGetParser(
79✔
113
    GDALFootprintOptions *psOptions,
114
    GDALFootprintOptionsForBinary *psOptionsForBinary)
115
{
116
    auto argParser = std::make_unique<GDALArgumentParser>(
117
        "gdal_footprint", /* bForBinary=*/psOptionsForBinary != nullptr);
79✔
118

119
    argParser->add_description(_("Compute footprint of a raster."));
79✔
120

121
    argParser->add_epilog(_("For more details, consult "
79✔
122
                            "https://gdal.org/programs/gdal_footprint.html"));
79✔
123

124
    argParser->add_argument("-b")
79✔
125
        .metavar("<band>")
158✔
126
        .scan<'i', int>()
79✔
127
        .append()
79✔
128
        .store_into(psOptions->anBands)
79✔
129
        .help(_("Band(s) of interest."));
79✔
130

131
    argParser->add_argument("-combine_bands")
79✔
132
        .choices("union", "intersection")
79✔
133
        .action([psOptions](const auto &s)
33✔
134
                { psOptions->bCombineBandsUnion = s == "union"; })
112✔
135
        .default_value("union")
79✔
136
        .help(_("Defines how the mask bands of the selected bands are combined "
137
                "to generate a single mask band, before being vectorized."));
79✔
138

139
    {
140
        auto &group = argParser->add_mutually_exclusive_group();
79✔
141

142
        group.add_argument("-ovr")
79✔
143
            .metavar("<index>")
158✔
144
            .scan<'i', int>()
79✔
145
            .store_into(psOptions->nOvrIndex)
79✔
146
            .help(_("Defines which overview level of source file must be used, "
147
                    "when overviews are available on the source raster."));
79✔
148

149
        group.add_argument("-srcnodata")
79✔
150
            .metavar("\"<value>[ <value>]...\"")
158✔
151
            .store_into(psOptions->osSrcNoData)
79✔
152
            .help(_("Set nodata value(s) for input bands."));
79✔
153
    }
154

155
    argParser->add_argument("-t_cs")
79✔
156
        .choices("pixel", "georef")
79✔
157
        .default_value("georef")
79✔
158
        .action(
159
            [psOptions](const auto &s)
40✔
160
            {
161
                const bool georefSet{s == "georef"};
20✔
162
                psOptions->bOutCSGeoref = georefSet;
20✔
163
                psOptions->bOutCSGeorefRequested = georefSet;
20✔
164
            })
79✔
165
        .help(_("Target coordinate system."));
79✔
166

167
    // Note: no store_into (requires post validation)
168
    argParser->add_argument("-t_srs")
79✔
169
        .metavar("<srs_def>")
158✔
170
        .help(_("Target CRS of the output file.."));
79✔
171

172
    argParser->add_argument("-split_polys")
79✔
173
        .flag()
79✔
174
        .store_into(psOptions->bSplitPolys)
79✔
175
        .help(_("Split multipolygons into several features each one with a "
176
                "single polygon."));
79✔
177

178
    argParser->add_argument("-convex_hull")
79✔
179
        .flag()
79✔
180
        .store_into(psOptions->bConvexHull)
79✔
181
        .help(_("Compute the convex hull of the (multi)polygons."));
79✔
182

183
    argParser->add_argument("-densify")
79✔
184
        .metavar("<value>")
158✔
185
        .scan<'g', double>()
79✔
186
        .store_into(psOptions->dfDensifyDistance)
79✔
187
        .help(_("The specified value of this option is the maximum distance "
188
                "between 2 consecutive points of the output geometry. "));
79✔
189

190
    argParser->add_argument("-simplify")
79✔
191
        .metavar("<value>")
158✔
192
        .scan<'g', double>()
79✔
193
        .store_into(psOptions->dfSimplifyTolerance)
79✔
194
        .help(_("The specified value of this option is the tolerance used to "
195
                "merge consecutive points of the output geometry."));
79✔
196

197
    argParser->add_argument("-min_ring_area")
79✔
198
        .metavar("<value>")
158✔
199
        .scan<'g', double>()
79✔
200
        .store_into(psOptions->dfMinRingArea)
79✔
201
        .help(_("The specified value of this option is the minimum area of a "
202
                "ring to be considered."));
79✔
203

204
    // Note: no store_into (requires post validation)
205
    argParser->add_argument("-max_points")
79✔
206
        .metavar("<value>|unlimited")
158✔
207
        .default_value("100")
79✔
208
        .help(_("The maximum number of points in the output geometry."));
79✔
209

210
    if (psOptionsForBinary)
79✔
211
    {
212
        argParser->add_quiet_argument(&psOptionsForBinary->bQuiet);
7✔
213
        argParser->add_open_options_argument(
214
            psOptionsForBinary->aosOpenOptions);
7✔
215
    }
216

217
    argParser->add_output_format_argument(psOptions->osFormat);
79✔
218

219
    {
220
        auto &group = argParser->add_mutually_exclusive_group();
79✔
221

222
        group.add_argument("-location_field_name")
79✔
223
            .metavar("<field_name>")
158✔
224
            .default_value("location")
79✔
225
            .store_into(psOptions->osLocationFieldName)
79✔
226
            .help(_(
227
                "Specifies the name of the field in the resulting vector "
228
                "dataset where the path of the input dataset will be stored."));
79✔
229

230
        group.add_argument("-no_location")
79✔
231
            .flag()
79✔
232
            .store_into(psOptions->bClearLocation)
79✔
233
            .help(
234
                _("Turns off the writing of the path of the input dataset as a "
235
                  "field in the output vector dataset."));
79✔
236
    }
237

238
    argParser->add_argument("-write_absolute_path")
79✔
239
        .flag()
79✔
240
        .store_into(psOptions->bAbsolutePath)
79✔
241
        .help(_("Enables writing the absolute path of the input dataset."));
79✔
242

243
    argParser->add_layer_creation_options_argument(psOptions->aosLCO);
79✔
244

245
    argParser->add_dataset_creation_options_argument(psOptions->aosDSCO);
79✔
246

247
    argParser->add_argument("-lyr_name")
79✔
248
        .metavar("<value>")
158✔
249
        .store_into(psOptions->osDestLayerName)
79✔
250
        .help(_("Name of the target layer."));
79✔
251

252
    if (psOptionsForBinary)
79✔
253
    {
254
        argParser->add_argument("-overwrite")
7✔
255
            .flag()
7✔
256
            .store_into(psOptionsForBinary->bOverwrite)
7✔
257
            .help(_("Overwrite the target layer if it exists."));
7✔
258

259
        argParser->add_argument("src_filename")
7✔
260
            .metavar("<src_filename>")
14✔
261
            .store_into(psOptionsForBinary->osSource)
7✔
262
            .help(_("Source raster file name."));
7✔
263

264
        argParser->add_argument("dst_filename")
7✔
265
            .metavar("<dst_filename>")
14✔
266
            .store_into(psOptionsForBinary->osDest)
7✔
267
            .help(_("Destination vector file name."));
7✔
268
    }
269

270
    return argParser;
79✔
271
}
272

273
/************************************************************************/
274
/*                       GDALFootprintMaskBand                          */
275
/************************************************************************/
276

277
class GDALFootprintMaskBand final : public GDALRasterBand
278
{
279
    GDALRasterBand *m_poSrcBand = nullptr;
280

281
    CPL_DISALLOW_COPY_ASSIGN(GDALFootprintMaskBand)
282

283
  public:
284
    explicit GDALFootprintMaskBand(GDALRasterBand *poSrcBand)
53✔
285
        : m_poSrcBand(poSrcBand)
53✔
286
    {
287
        nRasterXSize = m_poSrcBand->GetXSize();
53✔
288
        nRasterYSize = m_poSrcBand->GetYSize();
53✔
289
        eDataType = GDT_Byte;
53✔
290
        m_poSrcBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
53✔
291
    }
53✔
292

293
  protected:
294
    CPLErr IReadBlock(int nBlockXOff, int nBlockYOff, void *pData) override;
295

296
    CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
297
                     int nYSize, void *pData, int nBufXSize, int nBufYSize,
298
                     GDALDataType eBufType, GSpacing nPixelSpace,
299
                     GSpacing nLineSpace,
300
                     GDALRasterIOExtraArg *psExtraArg) override;
301
};
302

303
CPLErr GDALFootprintMaskBand::IReadBlock(int nBlockXOff, int nBlockYOff,
×
304
                                         void *pData)
305
{
306
    int nWindowXSize;
307
    int nWindowYSize;
308
    m_poSrcBand->GetActualBlockSize(nBlockXOff, nBlockYOff, &nWindowXSize,
×
309
                                    &nWindowYSize);
310
    GDALRasterIOExtraArg sExtraArg;
311
    INIT_RASTERIO_EXTRA_ARG(sExtraArg);
×
312
    return IRasterIO(GF_Read, nBlockXOff * nBlockXSize,
×
313
                     nBlockYOff * nBlockYSize, nWindowXSize, nWindowYSize,
×
314
                     pData, nWindowXSize, nWindowYSize, GDT_Byte, 1,
315
                     nBlockXSize, &sExtraArg);
×
316
}
317

318
CPLErr GDALFootprintMaskBand::IRasterIO(
2,184✔
319
    GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
320
    void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
321
    GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
322
{
323
    if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&
2,184✔
324
        eBufType == GDT_Byte && nPixelSpace == 1)
1,092✔
325
    {
326
        // Request when band seen as the mask band for GDALPolygonize()
327

328
        if (m_poSrcBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize, pData,
1,092✔
329
                                  nBufXSize, nBufYSize, eBufType, nPixelSpace,
330
                                  nLineSpace, psExtraArg) != CE_None)
1,092✔
331
        {
332
            return CE_Failure;
×
333
        }
334
        GByte *pabyData = static_cast<GByte *>(pData);
1,092✔
335
        for (int iY = 0; iY < nYSize; ++iY)
2,184✔
336
        {
337
            for (int iX = 0; iX < nXSize; ++iX)
21,430✔
338
            {
339
                if (pabyData[iX])
20,338✔
340
                    pabyData[iX] = 1;
20,108✔
341
            }
342
            pabyData += nLineSpace;
1,092✔
343
        }
344

345
        return CE_None;
1,092✔
346
    }
347

348
    if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&
1,092✔
349
        eBufType == GDT_Int64 &&
1,092✔
350
        nPixelSpace == static_cast<int>(sizeof(int64_t)) &&
1,092✔
351
        (nLineSpace % nPixelSpace) == 0)
1,092✔
352
    {
353
        // Request when band seen as the value band for GDALPolygonize()
354

355
        if (m_poSrcBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize, pData,
1,092✔
356
                                  nBufXSize, nBufYSize, eBufType, nPixelSpace,
357
                                  nLineSpace, psExtraArg) != CE_None)
1,092✔
358
        {
359
            return CE_Failure;
×
360
        }
361
        int64_t *panData = static_cast<int64_t *>(pData);
1,092✔
362
        for (int iY = 0; iY < nYSize; ++iY)
2,184✔
363
        {
364
            for (int iX = 0; iX < nXSize; ++iX)
21,430✔
365
            {
366
                if (panData[iX])
20,338✔
367
                    panData[iX] = 1;
20,108✔
368
            }
369
            panData += (nLineSpace / nPixelSpace);
1,092✔
370
        }
371

372
        return CE_None;
1,092✔
373
    }
374

375
    return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
×
376
                                     pData, nBufXSize, nBufYSize, eBufType,
377
                                     nPixelSpace, nLineSpace, psExtraArg);
×
378
}
379

380
/************************************************************************/
381
/*                   GDALFootprintCombinedMaskBand                      */
382
/************************************************************************/
383

384
class GDALFootprintCombinedMaskBand final : public GDALRasterBand
385
{
386
    std::vector<GDALRasterBand *> m_apoSrcBands{};
387

388
    /** Whether to combine bands unioning (true) or intersecting (false) */
389
    bool m_bUnion = false;
390

391
  public:
392
    explicit GDALFootprintCombinedMaskBand(
7✔
393
        const std::vector<GDALRasterBand *> &apoSrcBands, bool bUnion)
394
        : m_apoSrcBands(apoSrcBands), m_bUnion(bUnion)
7✔
395
    {
396
        nRasterXSize = m_apoSrcBands[0]->GetXSize();
7✔
397
        nRasterYSize = m_apoSrcBands[0]->GetYSize();
7✔
398
        eDataType = GDT_Byte;
7✔
399
        m_apoSrcBands[0]->GetBlockSize(&nBlockXSize, &nBlockYSize);
7✔
400
    }
7✔
401

402
  protected:
403
    CPLErr IReadBlock(int nBlockXOff, int nBlockYOff, void *pData) override;
404

405
    CPLErr IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize,
406
                     int nYSize, void *pData, int nBufXSize, int nBufYSize,
407
                     GDALDataType eBufType, GSpacing nPixelSpace,
408
                     GSpacing nLineSpace,
409
                     GDALRasterIOExtraArg *psExtraArg) override;
410
};
411

412
CPLErr GDALFootprintCombinedMaskBand::IReadBlock(int nBlockXOff, int nBlockYOff,
×
413
                                                 void *pData)
414
{
415
    int nWindowXSize;
416
    int nWindowYSize;
417
    m_apoSrcBands[0]->GetActualBlockSize(nBlockXOff, nBlockYOff, &nWindowXSize,
×
418
                                         &nWindowYSize);
419
    GDALRasterIOExtraArg sExtraArg;
420
    INIT_RASTERIO_EXTRA_ARG(sExtraArg);
×
421
    return IRasterIO(GF_Read, nBlockXOff * nBlockXSize,
×
422
                     nBlockYOff * nBlockYSize, nWindowXSize, nWindowYSize,
×
423
                     pData, nWindowXSize, nWindowYSize, GDT_Byte, 1,
424
                     nBlockXSize, &sExtraArg);
×
425
}
426

427
CPLErr GDALFootprintCombinedMaskBand::IRasterIO(
28✔
428
    GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
429
    void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
430
    GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
431
{
432
    if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&
28✔
433
        eBufType == GDT_Byte && nPixelSpace == 1)
14✔
434
    {
435
        // Request when band seen as the mask band for GDALPolygonize()
436
        {
437
            GByte *pabyData = static_cast<GByte *>(pData);
14✔
438
            for (int iY = 0; iY < nYSize; ++iY)
28✔
439
            {
440
                memset(pabyData, m_bUnion ? 0 : 1, nXSize);
14✔
441
                pabyData += nLineSpace;
14✔
442
            }
443
        }
444

445
        std::vector<GByte> abyTmp(static_cast<size_t>(nXSize) * nYSize);
28✔
446
        for (auto poBand : m_apoSrcBands)
42✔
447
        {
448
            if (poBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,
28✔
449
                                 abyTmp.data(), nBufXSize, nBufYSize, GDT_Byte,
28✔
450
                                 1, nXSize, psExtraArg) != CE_None)
28✔
451
            {
452
                return CE_Failure;
×
453
            }
454
            GByte *pabyData = static_cast<GByte *>(pData);
28✔
455
            size_t iTmp = 0;
28✔
456
            for (int iY = 0; iY < nYSize; ++iY)
56✔
457
            {
458
                if (m_bUnion)
28✔
459
                {
460
                    for (int iX = 0; iX < nXSize; ++iX, ++iTmp)
60✔
461
                    {
462
                        if (abyTmp[iTmp])
44✔
463
                            pabyData[iX] = 1;
20✔
464
                    }
465
                }
466
                else
467
                {
468
                    for (int iX = 0; iX < nXSize; ++iX, ++iTmp)
48✔
469
                    {
470
                        if (abyTmp[iTmp] == 0)
36✔
471
                            pabyData[iX] = 0;
18✔
472
                    }
473
                }
474
                pabyData += nLineSpace;
28✔
475
            }
476
        }
477

478
        return CE_None;
14✔
479
    }
480

481
    if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&
14✔
482
        eBufType == GDT_Int64 &&
14✔
483
        nPixelSpace == static_cast<int>(sizeof(int64_t)) &&
14✔
484
        (nLineSpace % nPixelSpace) == 0)
14✔
485
    {
486
        // Request when band seen as the value band for GDALPolygonize()
487
        {
488
            int64_t *panData = static_cast<int64_t *>(pData);
14✔
489
            for (int iY = 0; iY < nYSize; ++iY)
28✔
490
            {
491
                if (m_bUnion)
14✔
492
                {
493
                    memset(panData, 0, nXSize * sizeof(int64_t));
8✔
494
                }
495
                else
496
                {
497
                    int64_t nOne = 1;
6✔
498
                    GDALCopyWords(&nOne, GDT_Int64, 0, panData, GDT_Int64,
6✔
499
                                  sizeof(int64_t), nXSize);
500
                }
501
                panData += (nLineSpace / nPixelSpace);
14✔
502
            }
503
        }
504

505
        std::vector<GByte> abyTmp(static_cast<size_t>(nXSize) * nYSize);
28✔
506
        for (auto poBand : m_apoSrcBands)
42✔
507
        {
508
            if (poBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,
28✔
509
                                 abyTmp.data(), nBufXSize, nBufYSize, GDT_Byte,
28✔
510
                                 1, nXSize, psExtraArg) != CE_None)
28✔
511
            {
512
                return CE_Failure;
×
513
            }
514
            size_t iTmp = 0;
28✔
515
            int64_t *panData = static_cast<int64_t *>(pData);
28✔
516
            for (int iY = 0; iY < nYSize; ++iY)
56✔
517
            {
518
                if (m_bUnion)
28✔
519
                {
520
                    for (int iX = 0; iX < nXSize; ++iX, ++iTmp)
60✔
521
                    {
522
                        if (abyTmp[iTmp])
44✔
523
                            panData[iX] = 1;
20✔
524
                    }
525
                }
526
                else
527
                {
528
                    for (int iX = 0; iX < nXSize; ++iX, ++iTmp)
48✔
529
                    {
530
                        if (abyTmp[iTmp] == 0)
36✔
531
                            panData[iX] = 0;
18✔
532
                    }
533
                }
534
                panData += (nLineSpace / nPixelSpace);
28✔
535
            }
536
        }
537
        return CE_None;
14✔
538
    }
539

540
    return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
×
541
                                     pData, nBufXSize, nBufYSize, eBufType,
542
                                     nPixelSpace, nLineSpace, psExtraArg);
×
543
}
544

545
/************************************************************************/
546
/*                    GetOutputLayerAndUpdateDstDS()                    */
547
/************************************************************************/
548

549
static OGRLayer *
550
GetOutputLayerAndUpdateDstDS(const char *pszDest, GDALDatasetH &hDstDS,
74✔
551
                             GDALDataset *poSrcDS,
552
                             const GDALFootprintOptions *psOptions)
553
{
554

555
    if (pszDest == nullptr)
74✔
556
        pszDest = GDALGetDescription(hDstDS);
3✔
557

558
    /* -------------------------------------------------------------------- */
559
    /*      Create output dataset if needed                                 */
560
    /* -------------------------------------------------------------------- */
561
    const bool bCreateOutput = psOptions->bCreateOutput || hDstDS == nullptr;
74✔
562

563
    GDALDriverH hDriver = nullptr;
74✔
564
    if (bCreateOutput)
74✔
565
    {
566
        std::string osFormat(psOptions->osFormat);
68✔
567
        if (osFormat.empty())
68✔
568
        {
569
            const auto aoDrivers = GetOutputDriversFor(pszDest, GDAL_OF_VECTOR);
8✔
570
            if (aoDrivers.empty())
8✔
571
            {
572
                CPLError(CE_Failure, CPLE_AppDefined,
1✔
573
                         "Cannot guess driver for %s", pszDest);
574
                return nullptr;
1✔
575
            }
576
            else
577
            {
578
                if (aoDrivers.size() > 1)
7✔
579
                {
580
                    CPLError(CE_Warning, CPLE_AppDefined,
2✔
581
                             "Several drivers matching %s extension. Using %s",
582
                             CPLGetExtensionSafe(pszDest).c_str(),
2✔
583
                             aoDrivers[0].c_str());
1✔
584
                }
585
                osFormat = aoDrivers[0];
7✔
586
            }
587
        }
588

589
        /* ----------------------------------------------------------------- */
590
        /*      Find the output driver. */
591
        /* ----------------------------------------------------------------- */
592
        hDriver = GDALGetDriverByName(osFormat.c_str());
67✔
593
        char **papszDriverMD =
594
            hDriver ? GDALGetMetadata(hDriver, nullptr) : nullptr;
67✔
595
        if (hDriver == nullptr ||
133✔
596
            !CPLTestBool(CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_VECTOR,
66✔
597
                                              "FALSE")) ||
133✔
598
            !CPLTestBool(
65✔
599
                CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_CREATE, "FALSE")))
600
        {
601
            CPLError(CE_Failure, CPLE_NotSupported,
2✔
602
                     "Output driver `%s' not recognised or does not support "
603
                     "direct output file creation.",
604
                     osFormat.c_str());
605
            return nullptr;
2✔
606
        }
607

608
        hDstDS = GDALCreate(hDriver, pszDest, 0, 0, 0, GDT_Unknown,
65✔
609
                            psOptions->aosDSCO.List());
610
        if (!hDstDS)
65✔
611
        {
612
            return nullptr;
1✔
613
        }
614
    }
615

616
    /* -------------------------------------------------------------------- */
617
    /*      Open or create target layer.                                    */
618
    /* -------------------------------------------------------------------- */
619
    auto poDstDS = GDALDataset::FromHandle(hDstDS);
70✔
620
    OGRLayer *poLayer = nullptr;
70✔
621

622
    if (!bCreateOutput)
70✔
623
    {
624
        if (poDstDS->GetLayerCount() == 1 && poDstDS->GetDriver() &&
10✔
625
            EQUAL(poDstDS->GetDriver()->GetDescription(), "ESRI Shapefile"))
4✔
626
        {
627
            poLayer = poDstDS->GetLayer(0);
3✔
628
        }
629
        else if (!psOptions->osDestLayerName.empty())
3✔
630
        {
631
            poLayer =
632
                poDstDS->GetLayerByName(psOptions->osDestLayerName.c_str());
3✔
633
            if (!poLayer)
3✔
634
            {
635
                CPLError(CE_Failure, CPLE_AppDefined, "Cannot find layer %s",
1✔
636
                         psOptions->osDestLayerName.c_str());
637
                return nullptr;
1✔
638
            }
639
        }
640
        else
641
        {
642
            poLayer = poDstDS->GetLayerByName(DEFAULT_LAYER_NAME);
×
643
        }
644
    }
645

646
    if (!poLayer)
69✔
647
    {
648
        std::string osDestLayerName = psOptions->osDestLayerName;
64✔
649
        if (osDestLayerName.empty())
64✔
650
        {
651
            if (poDstDS->GetDriver() &&
118✔
652
                EQUAL(poDstDS->GetDriver()->GetDescription(), "ESRI Shapefile"))
59✔
653
            {
654
                osDestLayerName = CPLGetBasenameSafe(pszDest);
5✔
655
            }
656
            else
657
            {
658
                osDestLayerName = DEFAULT_LAYER_NAME;
54✔
659
            }
660
        }
661

662
        std::unique_ptr<OGRSpatialReference, OGRSpatialReferenceReleaser> poSRS;
×
663
        if (psOptions->bOutCSGeoref)
64✔
664
        {
665
            if (!psOptions->oOutputSRS.IsEmpty())
48✔
666
            {
667
                poSRS.reset(psOptions->oOutputSRS.Clone());
3✔
668
            }
669
            else if (auto poSrcSRS = poSrcDS->GetSpatialRef())
45✔
670
            {
671
                poSRS.reset(poSrcSRS->Clone());
17✔
672
            }
673
        }
674

675
        poLayer = poDstDS->CreateLayer(
192✔
676
            osDestLayerName.c_str(), poSRS.get(),
64✔
677
            psOptions->bSplitPolys ? wkbPolygon : wkbMultiPolygon,
64✔
678
            const_cast<char **>(psOptions->aosLCO.List()));
679

680
        if (!psOptions->osLocationFieldName.empty())
64✔
681
        {
682
            OGRFieldDefn oFieldDefn(psOptions->osLocationFieldName.c_str(),
683
                                    OFTString);
62✔
684
            if (poLayer->CreateField(&oFieldDefn) != OGRERR_NONE)
62✔
685
                return nullptr;
×
686
        }
687
    }
688

689
    return poLayer;
69✔
690
}
691

692
/************************************************************************/
693
/*                 GeoTransformCoordinateTransformation                 */
694
/************************************************************************/
695

696
class GeoTransformCoordinateTransformation final
697
    : public OGRCoordinateTransformation
698
{
699
    const GDALGeoTransform &m_gt;
700

701
  public:
702
    explicit GeoTransformCoordinateTransformation(const GDALGeoTransform &gt)
28✔
703
        : m_gt(gt)
28✔
704
    {
705
    }
28✔
706

707
    const OGRSpatialReference *GetSourceCS() const override;
708

709
    const OGRSpatialReference *GetTargetCS() const override
84✔
710
    {
711
        return nullptr;
84✔
712
    }
713

714
    OGRCoordinateTransformation *Clone() const override
×
715
    {
716
        return new GeoTransformCoordinateTransformation(m_gt);
×
717
    }
718

719
    OGRCoordinateTransformation *GetInverse() const override
×
720
    {
721
        CPLError(CE_Failure, CPLE_NotSupported,
×
722
                 "GeoTransformCoordinateTransformation::GetInverse() not "
723
                 "implemented");
724
        return nullptr;
×
725
    }
726

727
    int Transform(size_t nCount, double *x, double *y, double * /* z */,
28✔
728
                  double * /* t */, int *pabSuccess) override
729
    {
730
        for (size_t i = 0; i < nCount; ++i)
168✔
731
        {
732
            const double X = m_gt[0] + x[i] * m_gt[1] + y[i] * m_gt[2];
140✔
733
            const double Y = m_gt[3] + x[i] * m_gt[4] + y[i] * m_gt[5];
140✔
734
            x[i] = X;
140✔
735
            y[i] = Y;
140✔
736
            if (pabSuccess)
140✔
737
                pabSuccess[i] = TRUE;
140✔
738
        }
739
        return TRUE;
28✔
740
    }
741
};
742

743
const OGRSpatialReference *
744
GeoTransformCoordinateTransformation::GetSourceCS() const
×
745
{
746
    return nullptr;
×
747
}
748

749
/************************************************************************/
750
/*                             CountPoints()                            */
751
/************************************************************************/
752

753
static size_t CountPoints(const OGRGeometry *poGeom)
142✔
754
{
755
    if (poGeom->getGeometryType() == wkbMultiPolygon)
142✔
756
    {
757
        size_t n = 0;
49✔
758
        for (auto *poPoly : poGeom->toMultiPolygon())
99✔
759
        {
760
            n += CountPoints(poPoly);
50✔
761
        }
762
        return n;
49✔
763
    }
764
    else if (poGeom->getGeometryType() == wkbPolygon)
93✔
765
    {
766
        size_t n = 0;
93✔
767
        for (auto *poRing : poGeom->toPolygon())
187✔
768
        {
769
            n += poRing->getNumPoints() - 1;
94✔
770
        }
771
        return n;
93✔
772
    }
773
    return 0;
×
774
}
775

776
/************************************************************************/
777
/*                   GetMinDistanceBetweenTwoPoints()                   */
778
/************************************************************************/
779

780
static double GetMinDistanceBetweenTwoPoints(const OGRGeometry *poGeom)
6✔
781
{
782
    if (poGeom->getGeometryType() == wkbMultiPolygon)
6✔
783
    {
784
        double v = std::numeric_limits<double>::max();
2✔
785
        for (auto *poPoly : poGeom->toMultiPolygon())
4✔
786
        {
787
            v = std::min(v, GetMinDistanceBetweenTwoPoints(poPoly));
2✔
788
        }
789
        return v;
2✔
790
    }
791
    else if (poGeom->getGeometryType() == wkbPolygon)
4✔
792
    {
793
        double v = std::numeric_limits<double>::max();
2✔
794
        for (auto *poRing : poGeom->toPolygon())
4✔
795
        {
796
            v = std::min(v, GetMinDistanceBetweenTwoPoints(poRing));
2✔
797
        }
798
        return v;
2✔
799
    }
800
    else if (poGeom->getGeometryType() == wkbLineString)
2✔
801
    {
802
        double v = std::numeric_limits<double>::max();
2✔
803
        const auto poLS = poGeom->toLineString();
2✔
804
        const int nNumPoints = poLS->getNumPoints();
2✔
805
        for (int i = 0; i < nNumPoints - 1; ++i)
44✔
806
        {
807
            const double x1 = poLS->getX(i);
42✔
808
            const double y1 = poLS->getY(i);
42✔
809
            const double x2 = poLS->getX(i + 1);
42✔
810
            const double y2 = poLS->getY(i + 1);
42✔
811
            const double d = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);
42✔
812
            if (d > 0)
42✔
813
                v = std::min(v, d);
42✔
814
        }
815
        return sqrt(v);
2✔
816
    }
817
    return 0;
×
818
}
819

820
/************************************************************************/
821
/*                       GDALFootprintProcess()                         */
822
/************************************************************************/
823

824
static bool GDALFootprintProcess(GDALDataset *poSrcDS, OGRLayer *poDstLayer,
69✔
825
                                 const GDALFootprintOptions *psOptions)
826
{
827
    std::unique_ptr<OGRCoordinateTransformation> poCT_SRS;
69✔
828
    const OGRSpatialReference *poDstSRS = poDstLayer->GetSpatialRef();
69✔
829
    if (!psOptions->oOutputSRS.IsEmpty())
69✔
830
        poDstSRS = &(psOptions->oOutputSRS);
3✔
831
    if (poDstSRS)
69✔
832
    {
833
        auto poSrcSRS = poSrcDS->GetSpatialRef();
22✔
834
        if (!poSrcSRS)
22✔
835
        {
836
            CPLError(CE_Failure, CPLE_AppDefined,
1✔
837
                     "Output layer has CRS, but input is not georeferenced");
838
            return false;
1✔
839
        }
840
        poCT_SRS.reset(OGRCreateCoordinateTransformation(poSrcSRS, poDstSRS));
21✔
841
        if (!poCT_SRS)
21✔
842
            return false;
×
843
    }
844

845
    std::vector<int> anBands = psOptions->anBands;
136✔
846
    const int nBandCount = poSrcDS->GetRasterCount();
68✔
847
    if (anBands.empty())
68✔
848
    {
849
        for (int i = 1; i <= nBandCount; ++i)
139✔
850
            anBands.push_back(i);
76✔
851
    }
852

853
    std::vector<GDALRasterBand *> apoSrcMaskBands;
136✔
854
    const CPLStringList aosSrcNoData(
855
        CSLTokenizeString2(psOptions->osSrcNoData.c_str(), " ", 0));
136✔
856
    std::vector<double> adfSrcNoData;
136✔
857
    if (!psOptions->osSrcNoData.empty())
68✔
858
    {
859
        if (aosSrcNoData.size() != 1 &&
6✔
860
            static_cast<size_t>(aosSrcNoData.size()) != anBands.size())
2✔
861
        {
862
            CPLError(CE_Failure, CPLE_AppDefined,
1✔
863
                     "Number of values in -srcnodata should be 1 or the number "
864
                     "of bands");
865
            return false;
1✔
866
        }
867
        for (int i = 0; i < aosSrcNoData.size(); ++i)
7✔
868
        {
869
            adfSrcNoData.emplace_back(CPLAtof(aosSrcNoData[i]));
4✔
870
        }
871
    }
872
    bool bGlobalMask = true;
67✔
873
    std::vector<std::unique_ptr<GDALRasterBand>> apoTmpNoDataMaskBands;
134✔
874
    for (size_t i = 0; i < anBands.size(); ++i)
142✔
875
    {
876
        const int nBand = anBands[i];
80✔
877
        if (nBand <= 0 || nBand > nBandCount)
80✔
878
        {
879
            CPLError(CE_Failure, CPLE_AppDefined, "Invalid band number: %d",
1✔
880
                     nBand);
881
            return false;
5✔
882
        }
883
        auto poBand = poSrcDS->GetRasterBand(nBand);
79✔
884
        if (!adfSrcNoData.empty())
79✔
885
        {
886
            bGlobalMask = false;
4✔
887
            apoTmpNoDataMaskBands.emplace_back(
888
                std::make_unique<GDALNoDataMaskBand>(
12✔
889
                    poBand, adfSrcNoData.size() == 1 ? adfSrcNoData[0]
6✔
890
                                                     : adfSrcNoData[i]));
6✔
891
            apoSrcMaskBands.push_back(apoTmpNoDataMaskBands.back().get());
4✔
892
        }
893
        else
894
        {
895
            GDALRasterBand *poMaskBand;
896
            const int nMaskFlags = poBand->GetMaskFlags();
75✔
897
            if (poBand->GetColorInterpretation() == GCI_AlphaBand)
75✔
898
            {
899
                poMaskBand = poBand;
2✔
900
            }
901
            else
902
            {
903
                if ((nMaskFlags & GMF_PER_DATASET) == 0)
73✔
904
                {
905
                    bGlobalMask = false;
65✔
906
                }
907
                poMaskBand = poBand->GetMaskBand();
73✔
908
            }
909
            if (psOptions->nOvrIndex >= 0)
75✔
910
            {
911
                if (nMaskFlags == GMF_NODATA)
11✔
912
                {
913
                    // If the mask band is based on nodata, we don't need
914
                    // to check the overviews of the mask band, but we
915
                    // can take the mask band of the overviews
916
                    auto poOvrBand = poBand->GetOverview(psOptions->nOvrIndex);
5✔
917
                    if (!poOvrBand)
5✔
918
                    {
919
                        if (poBand->GetOverviewCount() == 0)
2✔
920
                        {
921
                            CPLError(
1✔
922
                                CE_Failure, CPLE_AppDefined,
923
                                "Overview index %d invalid for this dataset. "
924
                                "Bands of this dataset have no "
925
                                "precomputed overviews",
926
                                psOptions->nOvrIndex);
1✔
927
                        }
928
                        else
929
                        {
930
                            CPLError(
1✔
931
                                CE_Failure, CPLE_AppDefined,
932
                                "Overview index %d invalid for this dataset. "
933
                                "Value should be in [0,%d] range",
934
                                psOptions->nOvrIndex,
1✔
935
                                poBand->GetOverviewCount() - 1);
1✔
936
                        }
937
                        return false;
4✔
938
                    }
939
                    if (poOvrBand->GetMaskFlags() != GMF_NODATA)
3✔
940
                    {
941
                        CPLError(CE_Failure, CPLE_AppDefined,
×
942
                                 "poOvrBand->GetMaskFlags() != GMF_NODATA");
943
                        return false;
×
944
                    }
945
                    poMaskBand = poOvrBand->GetMaskBand();
3✔
946
                }
947
                else
948
                {
949
                    poMaskBand = poMaskBand->GetOverview(psOptions->nOvrIndex);
6✔
950
                    if (!poMaskBand)
6✔
951
                    {
952
                        if (poBand->GetMaskBand()->GetOverviewCount() == 0)
2✔
953
                        {
954
                            CPLError(
1✔
955
                                CE_Failure, CPLE_AppDefined,
956
                                "Overview index %d invalid for this dataset. "
957
                                "Mask bands of this dataset have no "
958
                                "precomputed overviews",
959
                                psOptions->nOvrIndex);
1✔
960
                        }
961
                        else
962
                        {
963
                            CPLError(
1✔
964
                                CE_Failure, CPLE_AppDefined,
965
                                "Overview index %d invalid for this dataset. "
966
                                "Value should be in [0,%d] range",
967
                                psOptions->nOvrIndex,
1✔
968
                                poBand->GetMaskBand()->GetOverviewCount() - 1);
1✔
969
                        }
970
                        return false;
2✔
971
                    }
972
                }
973
            }
974
            apoSrcMaskBands.push_back(poMaskBand);
71✔
975
        }
976
    }
977

978
    std::unique_ptr<OGRCoordinateTransformation> poCT_GT;
62✔
979
    GDALGeoTransform gt;
62✔
980
    if (psOptions->bOutCSGeoref && poSrcDS->GetGeoTransform(gt) == CE_None)
62✔
981
    {
982
        auto poMaskBand = apoSrcMaskBands[0];
25✔
983
        gt.Rescale(double(poSrcDS->GetRasterXSize()) / poMaskBand->GetXSize(),
25✔
984
                   double(poSrcDS->GetRasterYSize()) / poMaskBand->GetYSize());
25✔
985
        poCT_GT = std::make_unique<GeoTransformCoordinateTransformation>(gt);
25✔
986
    }
987
    else if (psOptions->bOutCSGeorefRequested)
37✔
988
    {
989
        CPLError(CE_Failure, CPLE_AppDefined,
2✔
990
                 "Georeferenced coordinates requested, but "
991
                 "input dataset has no geotransform.");
992
        return false;
2✔
993
    }
994
    else if (psOptions->nOvrIndex >= 0)
35✔
995
    {
996
        // Transform from overview pixel coordinates to full resolution
997
        // pixel coordinates
998
        auto poMaskBand = apoSrcMaskBands[0];
3✔
999
        gt[1] = double(poSrcDS->GetRasterXSize()) / poMaskBand->GetXSize();
3✔
1000
        gt[2] = 0;
3✔
1001
        gt[4] = 0;
3✔
1002
        gt[5] = double(poSrcDS->GetRasterYSize()) / poMaskBand->GetYSize();
3✔
1003
        poCT_GT = std::make_unique<GeoTransformCoordinateTransformation>(gt);
3✔
1004
    }
1005

1006
    std::unique_ptr<GDALRasterBand> poMaskForRasterize;
60✔
1007
    if (bGlobalMask || anBands.size() == 1)
60✔
1008
    {
1009
        poMaskForRasterize =
1010
            std::make_unique<GDALFootprintMaskBand>(apoSrcMaskBands[0]);
53✔
1011
    }
1012
    else
1013
    {
1014
        poMaskForRasterize = std::make_unique<GDALFootprintCombinedMaskBand>(
7✔
1015
            apoSrcMaskBands, psOptions->bCombineBandsUnion);
7✔
1016
    }
1017

1018
    auto hBand = GDALRasterBand::ToHandle(poMaskForRasterize.get());
60✔
1019
    auto poMemLayer = std::make_unique<OGRMemLayer>("", nullptr, wkbUnknown);
120✔
1020
    const CPLErr eErr =
1021
        GDALPolygonize(hBand, hBand, OGRLayer::ToHandle(poMemLayer.get()),
60✔
1022
                       /* iPixValField = */ -1,
1023
                       /* papszOptions = */ nullptr, psOptions->pfnProgress,
60✔
1024
                       psOptions->pProgressData);
60✔
1025
    if (eErr != CE_None)
60✔
1026
    {
1027
        return false;
×
1028
    }
1029

1030
    if (!psOptions->bSplitPolys)
60✔
1031
    {
1032
        auto poMP = std::make_unique<OGRMultiPolygon>();
116✔
1033
        for (auto &&poFeature : poMemLayer.get())
116✔
1034
        {
1035
            auto poGeom =
1036
                std::unique_ptr<OGRGeometry>(poFeature->StealGeometry());
116✔
1037
            CPLAssert(poGeom);
58✔
1038
            if (poGeom->getGeometryType() == wkbPolygon)
58✔
1039
            {
1040
                poMP->addGeometry(std::move(poGeom));
58✔
1041
            }
1042
        }
1043
        poMemLayer = std::make_unique<OGRMemLayer>("", nullptr, wkbUnknown);
58✔
1044
        auto poFeature =
1045
            std::make_unique<OGRFeature>(poMemLayer->GetLayerDefn());
58✔
1046
        poFeature->SetGeometryDirectly(poMP.release());
58✔
1047
        CPL_IGNORE_RET_VAL(poMemLayer->CreateFeature(poFeature.get()));
58✔
1048
    }
1049

1050
    for (auto &&poFeature : poMemLayer.get())
122✔
1051
    {
1052
        auto poGeom = std::unique_ptr<OGRGeometry>(poFeature->StealGeometry());
62✔
1053
        CPLAssert(poGeom);
62✔
1054
        if (poGeom->IsEmpty())
62✔
1055
            continue;
3✔
1056

1057
        auto poDstFeature =
1058
            std::make_unique<OGRFeature>(poDstLayer->GetLayerDefn());
59✔
1059
        poDstFeature->SetFrom(poFeature.get());
59✔
1060

1061
        if (poCT_GT)
59✔
1062
        {
1063
            if (poGeom->transform(poCT_GT.get()) != OGRERR_NONE)
28✔
1064
                return false;
×
1065
        }
1066

1067
        if (psOptions->dfDensifyDistance > 0)
59✔
1068
        {
1069
            OGREnvelope sEnvelope;
2✔
1070
            poGeom->getEnvelope(&sEnvelope);
2✔
1071
            // Some sanity check to avoid insane memory allocations
1072
            if (sEnvelope.MaxX - sEnvelope.MinX >
2✔
1073
                    1e6 * psOptions->dfDensifyDistance ||
2✔
1074
                sEnvelope.MaxY - sEnvelope.MinY >
2✔
1075
                    1e6 * psOptions->dfDensifyDistance)
2✔
1076
            {
1077
                CPLError(CE_Failure, CPLE_AppDefined,
×
1078
                         "Densification distance too small compared to "
1079
                         "geometry extent");
1080
                return false;
×
1081
            }
1082
            poGeom->segmentize(psOptions->dfDensifyDistance);
2✔
1083
        }
1084

1085
        if (poCT_SRS)
59✔
1086
        {
1087
            if (poGeom->transform(poCT_SRS.get()) != OGRERR_NONE)
21✔
1088
                return false;
×
1089
        }
1090

1091
        if (psOptions->dfMinRingArea != 0)
59✔
1092
        {
1093
            if (poGeom->getGeometryType() == wkbMultiPolygon)
4✔
1094
            {
1095
                auto poMP = std::make_unique<OGRMultiPolygon>();
8✔
1096
                for (auto *poPoly : poGeom->toMultiPolygon())
8✔
1097
                {
1098
                    auto poNewPoly = std::make_unique<OGRPolygon>();
8✔
1099
                    for (auto *poRing : poPoly)
8✔
1100
                    {
1101
                        if (poRing->get_Area() >= psOptions->dfMinRingArea)
4✔
1102
                        {
1103
                            poNewPoly->addRing(poRing);
2✔
1104
                        }
1105
                    }
1106
                    if (!poNewPoly->IsEmpty())
4✔
1107
                        poMP->addGeometry(std::move(poNewPoly));
2✔
1108
                }
1109
                poGeom = std::move(poMP);
4✔
1110
            }
1111
            else if (poGeom->getGeometryType() == wkbPolygon)
×
1112
            {
1113
                auto poNewPoly = std::make_unique<OGRPolygon>();
×
1114
                for (auto *poRing : poGeom->toPolygon())
×
1115
                {
1116
                    if (poRing->get_Area() >= psOptions->dfMinRingArea)
×
1117
                    {
1118
                        poNewPoly->addRing(poRing);
×
1119
                    }
1120
                }
1121
                poGeom = std::move(poNewPoly);
×
1122
            }
1123
            if (poGeom->IsEmpty())
4✔
1124
                continue;
2✔
1125
        }
1126

1127
        if (psOptions->bConvexHull)
57✔
1128
        {
1129
            poGeom.reset(poGeom->ConvexHull());
2✔
1130
            if (!poGeom || poGeom->IsEmpty())
2✔
1131
                continue;
×
1132
        }
1133

1134
        if (psOptions->dfSimplifyTolerance != 0)
57✔
1135
        {
1136
            poGeom.reset(poGeom->Simplify(psOptions->dfSimplifyTolerance));
2✔
1137
            if (!poGeom || poGeom->IsEmpty())
2✔
1138
                continue;
×
1139
        }
1140

1141
        if (psOptions->nMaxPoints > 0 &&
114✔
1142
            CountPoints(poGeom.get()) >
57✔
1143
                static_cast<size_t>(psOptions->nMaxPoints))
57✔
1144
        {
1145
            OGREnvelope sEnvelope;
2✔
1146
            poGeom->getEnvelope(&sEnvelope);
2✔
1147
            double tolMin = GetMinDistanceBetweenTwoPoints(poGeom.get());
2✔
1148
            double tolMax = std::max(sEnvelope.MaxY - sEnvelope.MinY,
4✔
1149
                                     sEnvelope.MaxX - sEnvelope.MinX);
2✔
1150
            for (int i = 0; i < 20; ++i)
42✔
1151
            {
1152
                const double tol = (tolMin + tolMax) / 2;
40✔
1153
                std::unique_ptr<OGRGeometry> poSimplifiedGeom(
1154
                    poGeom->Simplify(tol));
40✔
1155
                if (!poSimplifiedGeom || poSimplifiedGeom->IsEmpty())
40✔
1156
                {
1157
                    tolMax = tol;
5✔
1158
                    continue;
5✔
1159
                }
1160
                const auto nPoints = CountPoints(poSimplifiedGeom.get());
35✔
1161
                if (nPoints == static_cast<size_t>(psOptions->nMaxPoints))
35✔
1162
                {
1163
                    tolMax = tol;
×
1164
                    break;
×
1165
                }
1166
                else if (nPoints < static_cast<size_t>(psOptions->nMaxPoints))
35✔
1167
                {
1168
                    tolMax = tol;
35✔
1169
                }
1170
                else
1171
                {
1172
                    tolMin = tol;
×
1173
                }
1174
            }
1175
            poGeom.reset(poGeom->Simplify(tolMax));
2✔
1176
            if (!poGeom || poGeom->IsEmpty())
2✔
1177
                continue;
×
1178
        }
1179

1180
        if (!psOptions->bSplitPolys && poGeom->getGeometryType() == wkbPolygon)
57✔
1181
            poGeom.reset(
6✔
1182
                OGRGeometryFactory::forceToMultiPolygon(poGeom.release()));
1183

1184
        poDstFeature->SetGeometryDirectly(poGeom.release());
57✔
1185

1186
        if (!psOptions->osLocationFieldName.empty())
57✔
1187
        {
1188

1189
            std::string osFilename = poSrcDS->GetDescription();
110✔
1190
            // Make sure it is a file before building absolute path name.
1191
            VSIStatBufL sStatBuf;
1192
            if (psOptions->bAbsolutePath &&
112✔
1193
                CPLIsFilenameRelative(osFilename.c_str()) &&
57✔
1194
                VSIStatL(osFilename.c_str(), &sStatBuf) == 0)
2✔
1195
            {
1196
                char *pszCurDir = CPLGetCurrentDir();
2✔
1197
                if (pszCurDir)
2✔
1198
                {
1199
                    osFilename = CPLProjectRelativeFilenameSafe(
4✔
1200
                        pszCurDir, osFilename.c_str());
2✔
1201
                    CPLFree(pszCurDir);
2✔
1202
                }
1203
            }
1204
            poDstFeature->SetField(psOptions->osLocationFieldName.c_str(),
55✔
1205
                                   osFilename.c_str());
1206
        }
1207

1208
        if (poDstLayer->CreateFeature(poDstFeature.get()) != OGRERR_NONE)
57✔
1209
        {
1210
            return false;
×
1211
        }
1212
    }
1213

1214
    return true;
60✔
1215
}
1216

1217
/************************************************************************/
1218
/*                  GDALFootprintAppGetParserUsage()                    */
1219
/************************************************************************/
1220

1221
std::string GDALFootprintAppGetParserUsage()
×
1222
{
1223
    try
1224
    {
1225
        GDALFootprintOptions sOptions;
×
1226
        GDALFootprintOptionsForBinary sOptionsForBinary;
×
1227
        auto argParser =
1228
            GDALFootprintAppOptionsGetParser(&sOptions, &sOptionsForBinary);
×
1229
        return argParser->usage();
×
1230
    }
1231
    catch (const std::exception &err)
×
1232
    {
1233
        CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
×
1234
                 err.what());
×
1235
        return std::string();
×
1236
    }
1237
}
1238

1239
/************************************************************************/
1240
/*                             GDALFootprint()                          */
1241
/************************************************************************/
1242

1243
/* clang-format off */
1244
/**
1245
 * Computes the footprint of a raster.
1246
 *
1247
 * This is the equivalent of the
1248
 * <a href="/programs/gdal_footprint.html">gdal_footprint</a> utility.
1249
 *
1250
 * GDALFootprintOptions* must be allocated and freed with
1251
 * GDALFootprintOptionsNew() and GDALFootprintOptionsFree() respectively.
1252
 * pszDest and hDstDS cannot be used at the same time.
1253
 *
1254
 * @param pszDest the vector destination dataset path or NULL.
1255
 * @param hDstDS the vector destination dataset or NULL.
1256
 * @param hSrcDataset the raster source dataset handle.
1257
 * @param psOptionsIn the options struct returned by GDALFootprintOptionsNew()
1258
 * or NULL.
1259
 * @param pbUsageError pointer to a integer output variable to store if any
1260
 * usage error has occurred or NULL.
1261
 * @return the output dataset (new dataset that must be closed using
1262
 * GDALClose(), or hDstDS is not NULL) or NULL in case of error.
1263
 *
1264
 * @since GDAL 3.8
1265
 */
1266
/* clang-format on */
1267

1268
GDALDatasetH GDALFootprint(const char *pszDest, GDALDatasetH hDstDS,
78✔
1269
                           GDALDatasetH hSrcDataset,
1270
                           const GDALFootprintOptions *psOptionsIn,
1271
                           int *pbUsageError)
1272
{
1273
    if (pszDest == nullptr && hDstDS == nullptr)
78✔
1274
    {
1275
        CPLError(CE_Failure, CPLE_AppDefined,
1✔
1276
                 "pszDest == NULL && hDstDS == NULL");
1277

1278
        if (pbUsageError)
1✔
1279
            *pbUsageError = TRUE;
×
1280
        return nullptr;
1✔
1281
    }
1282
    if (hSrcDataset == nullptr)
77✔
1283
    {
1284
        CPLError(CE_Failure, CPLE_AppDefined, "hSrcDataset== NULL");
1✔
1285

1286
        if (pbUsageError)
1✔
1287
            *pbUsageError = TRUE;
×
1288
        return nullptr;
1✔
1289
    }
1290
    if (hDstDS != nullptr && psOptionsIn && psOptionsIn->bCreateOutput)
76✔
1291
    {
1292
        CPLError(CE_Failure, CPLE_AppDefined,
1✔
1293
                 "hDstDS != NULL but options that imply creating a new dataset "
1294
                 "have been set.");
1295

1296
        if (pbUsageError)
1✔
1297
            *pbUsageError = TRUE;
×
1298
        return nullptr;
1✔
1299
    }
1300

1301
    GDALFootprintOptions *psOptionsToFree = nullptr;
75✔
1302
    const GDALFootprintOptions *psOptions = psOptionsIn;
75✔
1303
    if (psOptions == nullptr)
75✔
1304
    {
1305
        psOptionsToFree = GDALFootprintOptionsNew(nullptr, nullptr);
1✔
1306
        psOptions = psOptionsToFree;
1✔
1307
    }
1308

1309
    const bool bCloseOutDSOnError = hDstDS == nullptr;
75✔
1310

1311
    auto poSrcDS = GDALDataset::FromHandle(hSrcDataset);
75✔
1312
    if (poSrcDS->GetRasterCount() == 0)
75✔
1313
    {
1314
        CPLError(CE_Failure, CPLE_AppDefined,
1✔
1315
                 "Input dataset has no raster band.%s",
1316
                 poSrcDS->GetMetadata("SUBDATASETS")
1✔
1317
                     ? " You need to specify one subdataset."
1318
                     : "");
1319
        GDALFootprintOptionsFree(psOptionsToFree);
1✔
1320
        if (bCloseOutDSOnError)
1✔
1321
            GDALClose(hDstDS);
×
1322
        return nullptr;
1✔
1323
    }
1324

1325
    auto poLayer =
1326
        GetOutputLayerAndUpdateDstDS(pszDest, hDstDS, poSrcDS, psOptions);
74✔
1327
    if (!poLayer)
74✔
1328
    {
1329
        GDALFootprintOptionsFree(psOptionsToFree);
5✔
1330
        if (hDstDS && bCloseOutDSOnError)
5✔
1331
            GDALClose(hDstDS);
×
1332
        return nullptr;
5✔
1333
    }
1334

1335
    if (!GDALFootprintProcess(poSrcDS, poLayer, psOptions))
69✔
1336
    {
1337
        GDALFootprintOptionsFree(psOptionsToFree);
9✔
1338
        if (bCloseOutDSOnError)
9✔
1339
            GDALClose(hDstDS);
8✔
1340
        return nullptr;
9✔
1341
    }
1342

1343
    GDALFootprintOptionsFree(psOptionsToFree);
60✔
1344

1345
    return hDstDS;
60✔
1346
}
1347

1348
/************************************************************************/
1349
/*                           GDALFootprintOptionsNew()                  */
1350
/************************************************************************/
1351

1352
/**
1353
 * Allocates a GDALFootprintOptions struct.
1354
 *
1355
 * @param papszArgv NULL terminated list of options (potentially including
1356
 * filename and open options too), or NULL. The accepted options are the ones of
1357
 * the <a href="/programs/gdal_rasterize.html">gdal_rasterize</a> utility.
1358
 * @param psOptionsForBinary (output) may be NULL (and should generally be
1359
 * NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with
1360
 * GDALFootprintOptionsForBinaryNew() prior to this function. Will be filled
1361
 * with potentially present filename, open options,...
1362
 * @return pointer to the allocated GDALFootprintOptions struct. Must be freed
1363
 * with GDALFootprintOptionsFree().
1364
 *
1365
 * @since GDAL 3.8
1366
 */
1367

1368
GDALFootprintOptions *
1369
GDALFootprintOptionsNew(char **papszArgv,
79✔
1370
                        GDALFootprintOptionsForBinary *psOptionsForBinary)
1371
{
1372
    auto psOptions = std::make_unique<GDALFootprintOptions>();
158✔
1373

1374
    /* -------------------------------------------------------------------- */
1375
    /*      Parse arguments.                                                */
1376
    /* -------------------------------------------------------------------- */
1377

1378
    CPLStringList aosArgv;
158✔
1379

1380
    if (papszArgv)
79✔
1381
    {
1382
        const int nArgc = CSLCount(papszArgv);
78✔
1383
        for (int i = 0; i < nArgc; i++)
682✔
1384
        {
1385
            aosArgv.AddString(papszArgv[i]);
604✔
1386
        }
1387
    }
1388

1389
    try
1390
    {
1391
        auto argParser = GDALFootprintAppOptionsGetParser(psOptions.get(),
1392
                                                          psOptionsForBinary);
80✔
1393

1394
        argParser->parse_args_without_binary_name(aosArgv.List());
79✔
1395

1396
        if (argParser->is_used("-t_srs"))
78✔
1397
        {
1398

1399
            const std::string osVal(argParser->get<std::string>("-t_srs"));
4✔
1400
            if (psOptions->oOutputSRS.SetFromUserInput(osVal.c_str()) !=
4✔
1401
                OGRERR_NONE)
1402
            {
1403
                CPLError(CE_Failure, CPLE_AppDefined,
×
1404
                         "Failed to process SRS definition: %s", osVal.c_str());
1405
                return nullptr;
×
1406
            }
1407
            psOptions->oOutputSRS.SetAxisMappingStrategy(
4✔
1408
                OAMS_TRADITIONAL_GIS_ORDER);
1409
        }
1410

1411
        if (argParser->is_used("-max_points"))
78✔
1412
        {
1413
            const std::string maxPoints{
1414
                argParser->get<std::string>("-max_points")};
32✔
1415
            if (maxPoints == "unlimited")
32✔
1416
            {
1417
                psOptions->nMaxPoints = 0;
×
1418
            }
1419
            else
1420
            {
1421
                psOptions->nMaxPoints = atoi(maxPoints.c_str());
32✔
1422
                if (psOptions->nMaxPoints > 0 && psOptions->nMaxPoints < 3)
32✔
1423
                {
1424
                    CPLError(CE_Failure, CPLE_NotSupported,
×
1425
                             "Invalid value for -max_points");
1426
                    return nullptr;
×
1427
                }
1428
            }
1429
        }
1430

1431
        psOptions->bCreateOutput = !psOptions->osFormat.empty();
78✔
1432
    }
1433
    catch (const std::exception &err)
1✔
1434
    {
1435
        CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
1✔
1436
                 err.what());
1✔
1437
        return nullptr;
1✔
1438
    }
1439

1440
    if (!psOptions->bOutCSGeoref && !psOptions->oOutputSRS.IsEmpty())
78✔
1441
    {
1442
        CPLError(CE_Failure, CPLE_AppDefined,
1✔
1443
                 "-t_cs pixel and -t_srs are mutually exclusive.");
1444
        return nullptr;
1✔
1445
    }
1446

1447
    if (psOptions->bClearLocation)
77✔
1448
    {
1449
        psOptions->osLocationFieldName.clear();
2✔
1450
    }
1451

1452
    if (psOptionsForBinary)
77✔
1453
    {
1454
        psOptionsForBinary->bCreateOutput = psOptions->bCreateOutput;
7✔
1455
        psOptionsForBinary->osFormat = psOptions->osFormat;
7✔
1456
        psOptionsForBinary->osDestLayerName = psOptions->osDestLayerName;
7✔
1457
    }
1458

1459
    return psOptions.release();
77✔
1460
}
1461

1462
/************************************************************************/
1463
/*                       GDALFootprintOptionsFree()                     */
1464
/************************************************************************/
1465

1466
/**
1467
 * Frees the GDALFootprintOptions struct.
1468
 *
1469
 * @param psOptions the options struct for GDALFootprint().
1470
 *
1471
 * @since GDAL 3.8
1472
 */
1473

1474
void GDALFootprintOptionsFree(GDALFootprintOptions *psOptions)
150✔
1475
{
1476
    delete psOptions;
150✔
1477
}
150✔
1478

1479
/************************************************************************/
1480
/*                  GDALFootprintOptionsSetProgress()                   */
1481
/************************************************************************/
1482

1483
/**
1484
 * Set a progress function.
1485
 *
1486
 * @param psOptions the options struct for GDALFootprint().
1487
 * @param pfnProgress the progress callback.
1488
 * @param pProgressData the user data for the progress callback.
1489
 *
1490
 * @since GDAL 3.8
1491
 */
1492

1493
void GDALFootprintOptionsSetProgress(GDALFootprintOptions *psOptions,
37✔
1494
                                     GDALProgressFunc pfnProgress,
1495
                                     void *pProgressData)
1496
{
1497
    psOptions->pfnProgress = pfnProgress ? pfnProgress : GDALDummyProgress;
37✔
1498
    psOptions->pProgressData = pProgressData;
37✔
1499
}
37✔
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc