LCOV - code coverage report
Current view: top level - apps - gdal_footprint_lib.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 563 635 88.7 %
Date: 2025-07-09 17:50:03 Functions: 18 24 75.0 %

          Line data    Source code
       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          79 : static std::unique_ptr<GDALArgumentParser> GDALFootprintAppOptionsGetParser(
     113             :     GDALFootprintOptions *psOptions,
     114             :     GDALFootprintOptionsForBinary *psOptionsForBinary)
     115             : {
     116             :     auto argParser = std::make_unique<GDALArgumentParser>(
     117          79 :         "gdal_footprint", /* bForBinary=*/psOptionsForBinary != nullptr);
     118             : 
     119          79 :     argParser->add_description(_("Compute footprint of a raster."));
     120             : 
     121          79 :     argParser->add_epilog(_("For more details, consult "
     122          79 :                             "https://gdal.org/programs/gdal_footprint.html"));
     123             : 
     124          79 :     argParser->add_argument("-b")
     125         158 :         .metavar("<band>")
     126          79 :         .scan<'i', int>()
     127          79 :         .append()
     128          79 :         .store_into(psOptions->anBands)
     129          79 :         .help(_("Band(s) of interest."));
     130             : 
     131          79 :     argParser->add_argument("-combine_bands")
     132          79 :         .choices("union", "intersection")
     133          33 :         .action([psOptions](const auto &s)
     134         112 :                 { psOptions->bCombineBandsUnion = s == "union"; })
     135          79 :         .default_value("union")
     136             :         .help(_("Defines how the mask bands of the selected bands are combined "
     137          79 :                 "to generate a single mask band, before being vectorized."));
     138             : 
     139             :     {
     140          79 :         auto &group = argParser->add_mutually_exclusive_group();
     141             : 
     142          79 :         group.add_argument("-ovr")
     143         158 :             .metavar("<index>")
     144          79 :             .scan<'i', int>()
     145          79 :             .store_into(psOptions->nOvrIndex)
     146             :             .help(_("Defines which overview level of source file must be used, "
     147          79 :                     "when overviews are available on the source raster."));
     148             : 
     149          79 :         group.add_argument("-srcnodata")
     150         158 :             .metavar("\"<value>[ <value>]...\"")
     151          79 :             .store_into(psOptions->osSrcNoData)
     152          79 :             .help(_("Set nodata value(s) for input bands."));
     153             :     }
     154             : 
     155          79 :     argParser->add_argument("-t_cs")
     156          79 :         .choices("pixel", "georef")
     157          79 :         .default_value("georef")
     158             :         .action(
     159          40 :             [psOptions](const auto &s)
     160             :             {
     161          20 :                 const bool georefSet{s == "georef"};
     162          20 :                 psOptions->bOutCSGeoref = georefSet;
     163          20 :                 psOptions->bOutCSGeorefRequested = georefSet;
     164          79 :             })
     165          79 :         .help(_("Target coordinate system."));
     166             : 
     167             :     // Note: no store_into (requires post validation)
     168          79 :     argParser->add_argument("-t_srs")
     169         158 :         .metavar("<srs_def>")
     170          79 :         .help(_("Target CRS of the output file.."));
     171             : 
     172          79 :     argParser->add_argument("-split_polys")
     173          79 :         .flag()
     174          79 :         .store_into(psOptions->bSplitPolys)
     175             :         .help(_("Split multipolygons into several features each one with a "
     176          79 :                 "single polygon."));
     177             : 
     178          79 :     argParser->add_argument("-convex_hull")
     179          79 :         .flag()
     180          79 :         .store_into(psOptions->bConvexHull)
     181          79 :         .help(_("Compute the convex hull of the (multi)polygons."));
     182             : 
     183          79 :     argParser->add_argument("-densify")
     184         158 :         .metavar("<value>")
     185          79 :         .scan<'g', double>()
     186          79 :         .store_into(psOptions->dfDensifyDistance)
     187             :         .help(_("The specified value of this option is the maximum distance "
     188          79 :                 "between 2 consecutive points of the output geometry. "));
     189             : 
     190          79 :     argParser->add_argument("-simplify")
     191         158 :         .metavar("<value>")
     192          79 :         .scan<'g', double>()
     193          79 :         .store_into(psOptions->dfSimplifyTolerance)
     194             :         .help(_("The specified value of this option is the tolerance used to "
     195          79 :                 "merge consecutive points of the output geometry."));
     196             : 
     197          79 :     argParser->add_argument("-min_ring_area")
     198         158 :         .metavar("<value>")
     199          79 :         .scan<'g', double>()
     200          79 :         .store_into(psOptions->dfMinRingArea)
     201             :         .help(_("The specified value of this option is the minimum area of a "
     202          79 :                 "ring to be considered."));
     203             : 
     204             :     // Note: no store_into (requires post validation)
     205          79 :     argParser->add_argument("-max_points")
     206         158 :         .metavar("<value>|unlimited")
     207          79 :         .default_value("100")
     208          79 :         .help(_("The maximum number of points in the output geometry."));
     209             : 
     210          79 :     if (psOptionsForBinary)
     211             :     {
     212           7 :         argParser->add_quiet_argument(&psOptionsForBinary->bQuiet);
     213             :         argParser->add_open_options_argument(
     214           7 :             psOptionsForBinary->aosOpenOptions);
     215             :     }
     216             : 
     217          79 :     argParser->add_output_format_argument(psOptions->osFormat);
     218             : 
     219             :     {
     220          79 :         auto &group = argParser->add_mutually_exclusive_group();
     221             : 
     222          79 :         group.add_argument("-location_field_name")
     223         158 :             .metavar("<field_name>")
     224          79 :             .default_value("location")
     225          79 :             .store_into(psOptions->osLocationFieldName)
     226             :             .help(_(
     227             :                 "Specifies the name of the field in the resulting vector "
     228          79 :                 "dataset where the path of the input dataset will be stored."));
     229             : 
     230          79 :         group.add_argument("-no_location")
     231          79 :             .flag()
     232          79 :             .store_into(psOptions->bClearLocation)
     233             :             .help(
     234             :                 _("Turns off the writing of the path of the input dataset as a "
     235          79 :                   "field in the output vector dataset."));
     236             :     }
     237             : 
     238          79 :     argParser->add_argument("-write_absolute_path")
     239          79 :         .flag()
     240          79 :         .store_into(psOptions->bAbsolutePath)
     241          79 :         .help(_("Enables writing the absolute path of the input dataset."));
     242             : 
     243          79 :     argParser->add_layer_creation_options_argument(psOptions->aosLCO);
     244             : 
     245          79 :     argParser->add_dataset_creation_options_argument(psOptions->aosDSCO);
     246             : 
     247          79 :     argParser->add_argument("-lyr_name")
     248         158 :         .metavar("<value>")
     249          79 :         .store_into(psOptions->osDestLayerName)
     250          79 :         .help(_("Name of the target layer."));
     251             : 
     252          79 :     if (psOptionsForBinary)
     253             :     {
     254           7 :         argParser->add_argument("-overwrite")
     255           7 :             .flag()
     256           7 :             .store_into(psOptionsForBinary->bOverwrite)
     257           7 :             .help(_("Overwrite the target layer if it exists."));
     258             : 
     259           7 :         argParser->add_argument("src_filename")
     260          14 :             .metavar("<src_filename>")
     261           7 :             .store_into(psOptionsForBinary->osSource)
     262           7 :             .help(_("Source raster file name."));
     263             : 
     264           7 :         argParser->add_argument("dst_filename")
     265          14 :             .metavar("<dst_filename>")
     266           7 :             .store_into(psOptionsForBinary->osDest)
     267           7 :             .help(_("Destination vector file name."));
     268             :     }
     269             : 
     270          79 :     return argParser;
     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          53 :     explicit GDALFootprintMaskBand(GDALRasterBand *poSrcBand)
     285          53 :         : m_poSrcBand(poSrcBand)
     286             :     {
     287          53 :         nRasterXSize = m_poSrcBand->GetXSize();
     288          53 :         nRasterYSize = m_poSrcBand->GetYSize();
     289          53 :         eDataType = GDT_Byte;
     290          53 :         m_poSrcBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
     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           0 : CPLErr GDALFootprintMaskBand::IReadBlock(int nBlockXOff, int nBlockYOff,
     304             :                                          void *pData)
     305             : {
     306             :     int nWindowXSize;
     307             :     int nWindowYSize;
     308           0 :     m_poSrcBand->GetActualBlockSize(nBlockXOff, nBlockYOff, &nWindowXSize,
     309             :                                     &nWindowYSize);
     310             :     GDALRasterIOExtraArg sExtraArg;
     311           0 :     INIT_RASTERIO_EXTRA_ARG(sExtraArg);
     312           0 :     return IRasterIO(GF_Read, nBlockXOff * nBlockXSize,
     313           0 :                      nBlockYOff * nBlockYSize, nWindowXSize, nWindowYSize,
     314             :                      pData, nWindowXSize, nWindowYSize, GDT_Byte, 1,
     315           0 :                      nBlockXSize, &sExtraArg);
     316             : }
     317             : 
     318        2184 : CPLErr GDALFootprintMaskBand::IRasterIO(
     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        2184 :     if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&
     324        1092 :         eBufType == GDT_Byte && nPixelSpace == 1)
     325             :     {
     326             :         // Request when band seen as the mask band for GDALPolygonize()
     327             : 
     328        1092 :         if (m_poSrcBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize, pData,
     329             :                                   nBufXSize, nBufYSize, eBufType, nPixelSpace,
     330        1092 :                                   nLineSpace, psExtraArg) != CE_None)
     331             :         {
     332           0 :             return CE_Failure;
     333             :         }
     334        1092 :         GByte *pabyData = static_cast<GByte *>(pData);
     335        2184 :         for (int iY = 0; iY < nYSize; ++iY)
     336             :         {
     337       21430 :             for (int iX = 0; iX < nXSize; ++iX)
     338             :             {
     339       20338 :                 if (pabyData[iX])
     340       20108 :                     pabyData[iX] = 1;
     341             :             }
     342        1092 :             pabyData += nLineSpace;
     343             :         }
     344             : 
     345        1092 :         return CE_None;
     346             :     }
     347             : 
     348        1092 :     if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&
     349        1092 :         eBufType == GDT_Int64 &&
     350        1092 :         nPixelSpace == static_cast<int>(sizeof(int64_t)) &&
     351        1092 :         (nLineSpace % nPixelSpace) == 0)
     352             :     {
     353             :         // Request when band seen as the value band for GDALPolygonize()
     354             : 
     355        1092 :         if (m_poSrcBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize, pData,
     356             :                                   nBufXSize, nBufYSize, eBufType, nPixelSpace,
     357        1092 :                                   nLineSpace, psExtraArg) != CE_None)
     358             :         {
     359           0 :             return CE_Failure;
     360             :         }
     361        1092 :         int64_t *panData = static_cast<int64_t *>(pData);
     362        2184 :         for (int iY = 0; iY < nYSize; ++iY)
     363             :         {
     364       21430 :             for (int iX = 0; iX < nXSize; ++iX)
     365             :             {
     366       20338 :                 if (panData[iX])
     367       20108 :                     panData[iX] = 1;
     368             :             }
     369        1092 :             panData += (nLineSpace / nPixelSpace);
     370             :         }
     371             : 
     372        1092 :         return CE_None;
     373             :     }
     374             : 
     375           0 :     return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
     376             :                                      pData, nBufXSize, nBufYSize, eBufType,
     377           0 :                                      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           7 :     explicit GDALFootprintCombinedMaskBand(
     393             :         const std::vector<GDALRasterBand *> &apoSrcBands, bool bUnion)
     394           7 :         : m_apoSrcBands(apoSrcBands), m_bUnion(bUnion)
     395             :     {
     396           7 :         nRasterXSize = m_apoSrcBands[0]->GetXSize();
     397           7 :         nRasterYSize = m_apoSrcBands[0]->GetYSize();
     398           7 :         eDataType = GDT_Byte;
     399           7 :         m_apoSrcBands[0]->GetBlockSize(&nBlockXSize, &nBlockYSize);
     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           0 : CPLErr GDALFootprintCombinedMaskBand::IReadBlock(int nBlockXOff, int nBlockYOff,
     413             :                                                  void *pData)
     414             : {
     415             :     int nWindowXSize;
     416             :     int nWindowYSize;
     417           0 :     m_apoSrcBands[0]->GetActualBlockSize(nBlockXOff, nBlockYOff, &nWindowXSize,
     418             :                                          &nWindowYSize);
     419             :     GDALRasterIOExtraArg sExtraArg;
     420           0 :     INIT_RASTERIO_EXTRA_ARG(sExtraArg);
     421           0 :     return IRasterIO(GF_Read, nBlockXOff * nBlockXSize,
     422           0 :                      nBlockYOff * nBlockYSize, nWindowXSize, nWindowYSize,
     423             :                      pData, nWindowXSize, nWindowYSize, GDT_Byte, 1,
     424           0 :                      nBlockXSize, &sExtraArg);
     425             : }
     426             : 
     427          28 : CPLErr GDALFootprintCombinedMaskBand::IRasterIO(
     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          28 :     if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&
     433          14 :         eBufType == GDT_Byte && nPixelSpace == 1)
     434             :     {
     435             :         // Request when band seen as the mask band for GDALPolygonize()
     436             :         {
     437          14 :             GByte *pabyData = static_cast<GByte *>(pData);
     438          28 :             for (int iY = 0; iY < nYSize; ++iY)
     439             :             {
     440          14 :                 memset(pabyData, m_bUnion ? 0 : 1, nXSize);
     441          14 :                 pabyData += nLineSpace;
     442             :             }
     443             :         }
     444             : 
     445          28 :         std::vector<GByte> abyTmp(static_cast<size_t>(nXSize) * nYSize);
     446          42 :         for (auto poBand : m_apoSrcBands)
     447             :         {
     448          28 :             if (poBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,
     449          28 :                                  abyTmp.data(), nBufXSize, nBufYSize, GDT_Byte,
     450          28 :                                  1, nXSize, psExtraArg) != CE_None)
     451             :             {
     452           0 :                 return CE_Failure;
     453             :             }
     454          28 :             GByte *pabyData = static_cast<GByte *>(pData);
     455          28 :             size_t iTmp = 0;
     456          56 :             for (int iY = 0; iY < nYSize; ++iY)
     457             :             {
     458          28 :                 if (m_bUnion)
     459             :                 {
     460          60 :                     for (int iX = 0; iX < nXSize; ++iX, ++iTmp)
     461             :                     {
     462          44 :                         if (abyTmp[iTmp])
     463          20 :                             pabyData[iX] = 1;
     464             :                     }
     465             :                 }
     466             :                 else
     467             :                 {
     468          48 :                     for (int iX = 0; iX < nXSize; ++iX, ++iTmp)
     469             :                     {
     470          36 :                         if (abyTmp[iTmp] == 0)
     471          18 :                             pabyData[iX] = 0;
     472             :                     }
     473             :                 }
     474          28 :                 pabyData += nLineSpace;
     475             :             }
     476             :         }
     477             : 
     478          14 :         return CE_None;
     479             :     }
     480             : 
     481          14 :     if (eRWFlag == GF_Read && nXSize == nBufXSize && nYSize == nBufYSize &&
     482          14 :         eBufType == GDT_Int64 &&
     483          14 :         nPixelSpace == static_cast<int>(sizeof(int64_t)) &&
     484          14 :         (nLineSpace % nPixelSpace) == 0)
     485             :     {
     486             :         // Request when band seen as the value band for GDALPolygonize()
     487             :         {
     488          14 :             int64_t *panData = static_cast<int64_t *>(pData);
     489          28 :             for (int iY = 0; iY < nYSize; ++iY)
     490             :             {
     491          14 :                 if (m_bUnion)
     492             :                 {
     493           8 :                     memset(panData, 0, nXSize * sizeof(int64_t));
     494             :                 }
     495             :                 else
     496             :                 {
     497           6 :                     int64_t nOne = 1;
     498           6 :                     GDALCopyWords(&nOne, GDT_Int64, 0, panData, GDT_Int64,
     499             :                                   sizeof(int64_t), nXSize);
     500             :                 }
     501          14 :                 panData += (nLineSpace / nPixelSpace);
     502             :             }
     503             :         }
     504             : 
     505          28 :         std::vector<GByte> abyTmp(static_cast<size_t>(nXSize) * nYSize);
     506          42 :         for (auto poBand : m_apoSrcBands)
     507             :         {
     508          28 :             if (poBand->RasterIO(GF_Read, nXOff, nYOff, nXSize, nYSize,
     509          28 :                                  abyTmp.data(), nBufXSize, nBufYSize, GDT_Byte,
     510          28 :                                  1, nXSize, psExtraArg) != CE_None)
     511             :             {
     512           0 :                 return CE_Failure;
     513             :             }
     514          28 :             size_t iTmp = 0;
     515          28 :             int64_t *panData = static_cast<int64_t *>(pData);
     516          56 :             for (int iY = 0; iY < nYSize; ++iY)
     517             :             {
     518          28 :                 if (m_bUnion)
     519             :                 {
     520          60 :                     for (int iX = 0; iX < nXSize; ++iX, ++iTmp)
     521             :                     {
     522          44 :                         if (abyTmp[iTmp])
     523          20 :                             panData[iX] = 1;
     524             :                     }
     525             :                 }
     526             :                 else
     527             :                 {
     528          48 :                     for (int iX = 0; iX < nXSize; ++iX, ++iTmp)
     529             :                     {
     530          36 :                         if (abyTmp[iTmp] == 0)
     531          18 :                             panData[iX] = 0;
     532             :                     }
     533             :                 }
     534          28 :                 panData += (nLineSpace / nPixelSpace);
     535             :             }
     536             :         }
     537          14 :         return CE_None;
     538             :     }
     539             : 
     540           0 :     return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
     541             :                                      pData, nBufXSize, nBufYSize, eBufType,
     542           0 :                                      nPixelSpace, nLineSpace, psExtraArg);
     543             : }
     544             : 
     545             : /************************************************************************/
     546             : /*                    GetOutputLayerAndUpdateDstDS()                    */
     547             : /************************************************************************/
     548             : 
     549             : static OGRLayer *
     550          74 : GetOutputLayerAndUpdateDstDS(const char *pszDest, GDALDatasetH &hDstDS,
     551             :                              GDALDataset *poSrcDS,
     552             :                              const GDALFootprintOptions *psOptions)
     553             : {
     554             : 
     555          74 :     if (pszDest == nullptr)
     556           3 :         pszDest = GDALGetDescription(hDstDS);
     557             : 
     558             :     /* -------------------------------------------------------------------- */
     559             :     /*      Create output dataset if needed                                 */
     560             :     /* -------------------------------------------------------------------- */
     561          74 :     const bool bCreateOutput = psOptions->bCreateOutput || hDstDS == nullptr;
     562             : 
     563          74 :     GDALDriverH hDriver = nullptr;
     564          74 :     if (bCreateOutput)
     565             :     {
     566          68 :         std::string osFormat(psOptions->osFormat);
     567          68 :         if (osFormat.empty())
     568             :         {
     569           8 :             const auto aoDrivers = GetOutputDriversFor(pszDest, GDAL_OF_VECTOR);
     570           8 :             if (aoDrivers.empty())
     571             :             {
     572           1 :                 CPLError(CE_Failure, CPLE_AppDefined,
     573             :                          "Cannot guess driver for %s", pszDest);
     574           1 :                 return nullptr;
     575             :             }
     576             :             else
     577             :             {
     578           7 :                 if (aoDrivers.size() > 1)
     579             :                 {
     580           2 :                     CPLError(CE_Warning, CPLE_AppDefined,
     581             :                              "Several drivers matching %s extension. Using %s",
     582           2 :                              CPLGetExtensionSafe(pszDest).c_str(),
     583           1 :                              aoDrivers[0].c_str());
     584             :                 }
     585           7 :                 osFormat = aoDrivers[0];
     586             :             }
     587             :         }
     588             : 
     589             :         /* ----------------------------------------------------------------- */
     590             :         /*      Find the output driver. */
     591             :         /* ----------------------------------------------------------------- */
     592          67 :         hDriver = GDALGetDriverByName(osFormat.c_str());
     593             :         char **papszDriverMD =
     594          67 :             hDriver ? GDALGetMetadata(hDriver, nullptr) : nullptr;
     595         133 :         if (hDriver == nullptr ||
     596          66 :             !CPLTestBool(CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_VECTOR,
     597         133 :                                               "FALSE")) ||
     598          65 :             !CPLTestBool(
     599             :                 CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_CREATE, "FALSE")))
     600             :         {
     601           2 :             CPLError(CE_Failure, CPLE_NotSupported,
     602             :                      "Output driver `%s' not recognised or does not support "
     603             :                      "direct output file creation.",
     604             :                      osFormat.c_str());
     605           2 :             return nullptr;
     606             :         }
     607             : 
     608          65 :         hDstDS = GDALCreate(hDriver, pszDest, 0, 0, 0, GDT_Unknown,
     609             :                             psOptions->aosDSCO.List());
     610          65 :         if (!hDstDS)
     611             :         {
     612           1 :             return nullptr;
     613             :         }
     614             :     }
     615             : 
     616             :     /* -------------------------------------------------------------------- */
     617             :     /*      Open or create target layer.                                    */
     618             :     /* -------------------------------------------------------------------- */
     619          70 :     auto poDstDS = GDALDataset::FromHandle(hDstDS);
     620          70 :     OGRLayer *poLayer = nullptr;
     621             : 
     622          70 :     if (!bCreateOutput)
     623             :     {
     624          10 :         if (poDstDS->GetLayerCount() == 1 && poDstDS->GetDriver() &&
     625           4 :             EQUAL(poDstDS->GetDriver()->GetDescription(), "ESRI Shapefile"))
     626             :         {
     627           3 :             poLayer = poDstDS->GetLayer(0);
     628             :         }
     629           3 :         else if (!psOptions->osDestLayerName.empty())
     630             :         {
     631             :             poLayer =
     632           3 :                 poDstDS->GetLayerByName(psOptions->osDestLayerName.c_str());
     633           3 :             if (!poLayer)
     634             :             {
     635           1 :                 CPLError(CE_Failure, CPLE_AppDefined, "Cannot find layer %s",
     636             :                          psOptions->osDestLayerName.c_str());
     637           1 :                 return nullptr;
     638             :             }
     639             :         }
     640             :         else
     641             :         {
     642           0 :             poLayer = poDstDS->GetLayerByName(DEFAULT_LAYER_NAME);
     643             :         }
     644             :     }
     645             : 
     646          69 :     if (!poLayer)
     647             :     {
     648          64 :         std::string osDestLayerName = psOptions->osDestLayerName;
     649          64 :         if (osDestLayerName.empty())
     650             :         {
     651         118 :             if (poDstDS->GetDriver() &&
     652          59 :                 EQUAL(poDstDS->GetDriver()->GetDescription(), "ESRI Shapefile"))
     653             :             {
     654           5 :                 osDestLayerName = CPLGetBasenameSafe(pszDest);
     655             :             }
     656             :             else
     657             :             {
     658          54 :                 osDestLayerName = DEFAULT_LAYER_NAME;
     659             :             }
     660             :         }
     661             : 
     662           0 :         std::unique_ptr<OGRSpatialReference, OGRSpatialReferenceReleaser> poSRS;
     663          64 :         if (psOptions->bOutCSGeoref)
     664             :         {
     665          48 :             if (!psOptions->oOutputSRS.IsEmpty())
     666             :             {
     667           3 :                 poSRS.reset(psOptions->oOutputSRS.Clone());
     668             :             }
     669          45 :             else if (auto poSrcSRS = poSrcDS->GetSpatialRef())
     670             :             {
     671          17 :                 poSRS.reset(poSrcSRS->Clone());
     672             :             }
     673             :         }
     674             : 
     675         192 :         poLayer = poDstDS->CreateLayer(
     676          64 :             osDestLayerName.c_str(), poSRS.get(),
     677          64 :             psOptions->bSplitPolys ? wkbPolygon : wkbMultiPolygon,
     678             :             const_cast<char **>(psOptions->aosLCO.List()));
     679             : 
     680          64 :         if (!psOptions->osLocationFieldName.empty())
     681             :         {
     682             :             OGRFieldDefn oFieldDefn(psOptions->osLocationFieldName.c_str(),
     683          62 :                                     OFTString);
     684          62 :             if (poLayer->CreateField(&oFieldDefn) != OGRERR_NONE)
     685           0 :                 return nullptr;
     686             :         }
     687             :     }
     688             : 
     689          69 :     return poLayer;
     690             : }
     691             : 
     692             : /************************************************************************/
     693             : /*                 GeoTransformCoordinateTransformation                 */
     694             : /************************************************************************/
     695             : 
     696             : class GeoTransformCoordinateTransformation final
     697             :     : public OGRCoordinateTransformation
     698             : {
     699             :     const GDALGeoTransform &m_gt;
     700             : 
     701             :   public:
     702          28 :     explicit GeoTransformCoordinateTransformation(const GDALGeoTransform &gt)
     703          28 :         : m_gt(gt)
     704             :     {
     705          28 :     }
     706             : 
     707             :     const OGRSpatialReference *GetSourceCS() const override;
     708             : 
     709          84 :     const OGRSpatialReference *GetTargetCS() const override
     710             :     {
     711          84 :         return nullptr;
     712             :     }
     713             : 
     714           0 :     OGRCoordinateTransformation *Clone() const override
     715             :     {
     716           0 :         return new GeoTransformCoordinateTransformation(m_gt);
     717             :     }
     718             : 
     719           0 :     OGRCoordinateTransformation *GetInverse() const override
     720             :     {
     721           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     722             :                  "GeoTransformCoordinateTransformation::GetInverse() not "
     723             :                  "implemented");
     724           0 :         return nullptr;
     725             :     }
     726             : 
     727          28 :     int Transform(size_t nCount, double *x, double *y, double * /* z */,
     728             :                   double * /* t */, int *pabSuccess) override
     729             :     {
     730         168 :         for (size_t i = 0; i < nCount; ++i)
     731             :         {
     732         140 :             const double X = m_gt[0] + x[i] * m_gt[1] + y[i] * m_gt[2];
     733         140 :             const double Y = m_gt[3] + x[i] * m_gt[4] + y[i] * m_gt[5];
     734         140 :             x[i] = X;
     735         140 :             y[i] = Y;
     736         140 :             if (pabSuccess)
     737         140 :                 pabSuccess[i] = TRUE;
     738             :         }
     739          28 :         return TRUE;
     740             :     }
     741             : };
     742             : 
     743             : const OGRSpatialReference *
     744           0 : GeoTransformCoordinateTransformation::GetSourceCS() const
     745             : {
     746           0 :     return nullptr;
     747             : }
     748             : 
     749             : /************************************************************************/
     750             : /*                             CountPoints()                            */
     751             : /************************************************************************/
     752             : 
     753         142 : static size_t CountPoints(const OGRGeometry *poGeom)
     754             : {
     755         142 :     if (poGeom->getGeometryType() == wkbMultiPolygon)
     756             :     {
     757          49 :         size_t n = 0;
     758          99 :         for (auto *poPoly : poGeom->toMultiPolygon())
     759             :         {
     760          50 :             n += CountPoints(poPoly);
     761             :         }
     762          49 :         return n;
     763             :     }
     764          93 :     else if (poGeom->getGeometryType() == wkbPolygon)
     765             :     {
     766          93 :         size_t n = 0;
     767         187 :         for (auto *poRing : poGeom->toPolygon())
     768             :         {
     769          94 :             n += poRing->getNumPoints() - 1;
     770             :         }
     771          93 :         return n;
     772             :     }
     773           0 :     return 0;
     774             : }
     775             : 
     776             : /************************************************************************/
     777             : /*                   GetMinDistanceBetweenTwoPoints()                   */
     778             : /************************************************************************/
     779             : 
     780           6 : static double GetMinDistanceBetweenTwoPoints(const OGRGeometry *poGeom)
     781             : {
     782           6 :     if (poGeom->getGeometryType() == wkbMultiPolygon)
     783             :     {
     784           2 :         double v = std::numeric_limits<double>::max();
     785           4 :         for (auto *poPoly : poGeom->toMultiPolygon())
     786             :         {
     787           2 :             v = std::min(v, GetMinDistanceBetweenTwoPoints(poPoly));
     788             :         }
     789           2 :         return v;
     790             :     }
     791           4 :     else if (poGeom->getGeometryType() == wkbPolygon)
     792             :     {
     793           2 :         double v = std::numeric_limits<double>::max();
     794           4 :         for (auto *poRing : poGeom->toPolygon())
     795             :         {
     796           2 :             v = std::min(v, GetMinDistanceBetweenTwoPoints(poRing));
     797             :         }
     798           2 :         return v;
     799             :     }
     800           2 :     else if (poGeom->getGeometryType() == wkbLineString)
     801             :     {
     802           2 :         double v = std::numeric_limits<double>::max();
     803           2 :         const auto poLS = poGeom->toLineString();
     804           2 :         const int nNumPoints = poLS->getNumPoints();
     805          44 :         for (int i = 0; i < nNumPoints - 1; ++i)
     806             :         {
     807          42 :             const double x1 = poLS->getX(i);
     808          42 :             const double y1 = poLS->getY(i);
     809          42 :             const double x2 = poLS->getX(i + 1);
     810          42 :             const double y2 = poLS->getY(i + 1);
     811          42 :             const double d = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);
     812          42 :             if (d > 0)
     813          42 :                 v = std::min(v, d);
     814             :         }
     815           2 :         return sqrt(v);
     816             :     }
     817           0 :     return 0;
     818             : }
     819             : 
     820             : /************************************************************************/
     821             : /*                       GDALFootprintProcess()                         */
     822             : /************************************************************************/
     823             : 
     824          69 : static bool GDALFootprintProcess(GDALDataset *poSrcDS, OGRLayer *poDstLayer,
     825             :                                  const GDALFootprintOptions *psOptions)
     826             : {
     827          69 :     std::unique_ptr<OGRCoordinateTransformation> poCT_SRS;
     828          69 :     const OGRSpatialReference *poDstSRS = poDstLayer->GetSpatialRef();
     829          69 :     if (!psOptions->oOutputSRS.IsEmpty())
     830           3 :         poDstSRS = &(psOptions->oOutputSRS);
     831          69 :     if (poDstSRS)
     832             :     {
     833          22 :         auto poSrcSRS = poSrcDS->GetSpatialRef();
     834          22 :         if (!poSrcSRS)
     835             :         {
     836           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     837             :                      "Output layer has CRS, but input is not georeferenced");
     838           1 :             return false;
     839             :         }
     840          21 :         poCT_SRS.reset(OGRCreateCoordinateTransformation(poSrcSRS, poDstSRS));
     841          21 :         if (!poCT_SRS)
     842           0 :             return false;
     843             :     }
     844             : 
     845         136 :     std::vector<int> anBands = psOptions->anBands;
     846          68 :     const int nBandCount = poSrcDS->GetRasterCount();
     847          68 :     if (anBands.empty())
     848             :     {
     849         139 :         for (int i = 1; i <= nBandCount; ++i)
     850          76 :             anBands.push_back(i);
     851             :     }
     852             : 
     853         136 :     std::vector<GDALRasterBand *> apoSrcMaskBands;
     854             :     const CPLStringList aosSrcNoData(
     855         136 :         CSLTokenizeString2(psOptions->osSrcNoData.c_str(), " ", 0));
     856         136 :     std::vector<double> adfSrcNoData;
     857          68 :     if (!psOptions->osSrcNoData.empty())
     858             :     {
     859           6 :         if (aosSrcNoData.size() != 1 &&
     860           2 :             static_cast<size_t>(aosSrcNoData.size()) != anBands.size())
     861             :         {
     862           1 :             CPLError(CE_Failure, CPLE_AppDefined,
     863             :                      "Number of values in -srcnodata should be 1 or the number "
     864             :                      "of bands");
     865           1 :             return false;
     866             :         }
     867           7 :         for (int i = 0; i < aosSrcNoData.size(); ++i)
     868             :         {
     869           4 :             adfSrcNoData.emplace_back(CPLAtof(aosSrcNoData[i]));
     870             :         }
     871             :     }
     872          67 :     bool bGlobalMask = true;
     873         134 :     std::vector<std::unique_ptr<GDALRasterBand>> apoTmpNoDataMaskBands;
     874         142 :     for (size_t i = 0; i < anBands.size(); ++i)
     875             :     {
     876          80 :         const int nBand = anBands[i];
     877          80 :         if (nBand <= 0 || nBand > nBandCount)
     878             :         {
     879           1 :             CPLError(CE_Failure, CPLE_AppDefined, "Invalid band number: %d",
     880             :                      nBand);
     881           5 :             return false;
     882             :         }
     883          79 :         auto poBand = poSrcDS->GetRasterBand(nBand);
     884          79 :         if (!adfSrcNoData.empty())
     885             :         {
     886           4 :             bGlobalMask = false;
     887             :             apoTmpNoDataMaskBands.emplace_back(
     888          12 :                 std::make_unique<GDALNoDataMaskBand>(
     889           6 :                     poBand, adfSrcNoData.size() == 1 ? adfSrcNoData[0]
     890           6 :                                                      : adfSrcNoData[i]));
     891           4 :             apoSrcMaskBands.push_back(apoTmpNoDataMaskBands.back().get());
     892             :         }
     893             :         else
     894             :         {
     895             :             GDALRasterBand *poMaskBand;
     896          75 :             const int nMaskFlags = poBand->GetMaskFlags();
     897          75 :             if (poBand->GetColorInterpretation() == GCI_AlphaBand)
     898             :             {
     899           2 :                 poMaskBand = poBand;
     900             :             }
     901             :             else
     902             :             {
     903          73 :                 if ((nMaskFlags & GMF_PER_DATASET) == 0)
     904             :                 {
     905          65 :                     bGlobalMask = false;
     906             :                 }
     907          73 :                 poMaskBand = poBand->GetMaskBand();
     908             :             }
     909          75 :             if (psOptions->nOvrIndex >= 0)
     910             :             {
     911          11 :                 if (nMaskFlags == GMF_NODATA)
     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           5 :                     auto poOvrBand = poBand->GetOverview(psOptions->nOvrIndex);
     917           5 :                     if (!poOvrBand)
     918             :                     {
     919           2 :                         if (poBand->GetOverviewCount() == 0)
     920             :                         {
     921           1 :                             CPLError(
     922             :                                 CE_Failure, CPLE_AppDefined,
     923             :                                 "Overview index %d invalid for this dataset. "
     924             :                                 "Bands of this dataset have no "
     925             :                                 "precomputed overviews",
     926           1 :                                 psOptions->nOvrIndex);
     927             :                         }
     928             :                         else
     929             :                         {
     930           1 :                             CPLError(
     931             :                                 CE_Failure, CPLE_AppDefined,
     932             :                                 "Overview index %d invalid for this dataset. "
     933             :                                 "Value should be in [0,%d] range",
     934           1 :                                 psOptions->nOvrIndex,
     935           1 :                                 poBand->GetOverviewCount() - 1);
     936             :                         }
     937           4 :                         return false;
     938             :                     }
     939           3 :                     if (poOvrBand->GetMaskFlags() != GMF_NODATA)
     940             :                     {
     941           0 :                         CPLError(CE_Failure, CPLE_AppDefined,
     942             :                                  "poOvrBand->GetMaskFlags() != GMF_NODATA");
     943           0 :                         return false;
     944             :                     }
     945           3 :                     poMaskBand = poOvrBand->GetMaskBand();
     946             :                 }
     947             :                 else
     948             :                 {
     949           6 :                     poMaskBand = poMaskBand->GetOverview(psOptions->nOvrIndex);
     950           6 :                     if (!poMaskBand)
     951             :                     {
     952           2 :                         if (poBand->GetMaskBand()->GetOverviewCount() == 0)
     953             :                         {
     954           1 :                             CPLError(
     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           1 :                                 psOptions->nOvrIndex);
     960             :                         }
     961             :                         else
     962             :                         {
     963           1 :                             CPLError(
     964             :                                 CE_Failure, CPLE_AppDefined,
     965             :                                 "Overview index %d invalid for this dataset. "
     966             :                                 "Value should be in [0,%d] range",
     967           1 :                                 psOptions->nOvrIndex,
     968           1 :                                 poBand->GetMaskBand()->GetOverviewCount() - 1);
     969             :                         }
     970           2 :                         return false;
     971             :                     }
     972             :                 }
     973             :             }
     974          71 :             apoSrcMaskBands.push_back(poMaskBand);
     975             :         }
     976             :     }
     977             : 
     978          62 :     std::unique_ptr<OGRCoordinateTransformation> poCT_GT;
     979          62 :     GDALGeoTransform gt;
     980          62 :     if (psOptions->bOutCSGeoref && poSrcDS->GetGeoTransform(gt) == CE_None)
     981             :     {
     982          25 :         auto poMaskBand = apoSrcMaskBands[0];
     983          25 :         gt.Rescale(double(poSrcDS->GetRasterXSize()) / poMaskBand->GetXSize(),
     984          25 :                    double(poSrcDS->GetRasterYSize()) / poMaskBand->GetYSize());
     985          25 :         poCT_GT = std::make_unique<GeoTransformCoordinateTransformation>(gt);
     986             :     }
     987          37 :     else if (psOptions->bOutCSGeorefRequested)
     988             :     {
     989           2 :         CPLError(CE_Failure, CPLE_AppDefined,
     990             :                  "Georeferenced coordinates requested, but "
     991             :                  "input dataset has no geotransform.");
     992           2 :         return false;
     993             :     }
     994          35 :     else if (psOptions->nOvrIndex >= 0)
     995             :     {
     996             :         // Transform from overview pixel coordinates to full resolution
     997             :         // pixel coordinates
     998           3 :         auto poMaskBand = apoSrcMaskBands[0];
     999           3 :         gt[1] = double(poSrcDS->GetRasterXSize()) / poMaskBand->GetXSize();
    1000           3 :         gt[2] = 0;
    1001           3 :         gt[4] = 0;
    1002           3 :         gt[5] = double(poSrcDS->GetRasterYSize()) / poMaskBand->GetYSize();
    1003           3 :         poCT_GT = std::make_unique<GeoTransformCoordinateTransformation>(gt);
    1004             :     }
    1005             : 
    1006          60 :     std::unique_ptr<GDALRasterBand> poMaskForRasterize;
    1007          60 :     if (bGlobalMask || anBands.size() == 1)
    1008             :     {
    1009             :         poMaskForRasterize =
    1010          53 :             std::make_unique<GDALFootprintMaskBand>(apoSrcMaskBands[0]);
    1011             :     }
    1012             :     else
    1013             :     {
    1014           7 :         poMaskForRasterize = std::make_unique<GDALFootprintCombinedMaskBand>(
    1015           7 :             apoSrcMaskBands, psOptions->bCombineBandsUnion);
    1016             :     }
    1017             : 
    1018          60 :     auto hBand = GDALRasterBand::ToHandle(poMaskForRasterize.get());
    1019         120 :     auto poMemLayer = std::make_unique<OGRMemLayer>("", nullptr, wkbUnknown);
    1020             :     const CPLErr eErr =
    1021          60 :         GDALPolygonize(hBand, hBand, OGRLayer::ToHandle(poMemLayer.get()),
    1022             :                        /* iPixValField = */ -1,
    1023          60 :                        /* papszOptions = */ nullptr, psOptions->pfnProgress,
    1024          60 :                        psOptions->pProgressData);
    1025          60 :     if (eErr != CE_None)
    1026             :     {
    1027           0 :         return false;
    1028             :     }
    1029             : 
    1030          60 :     if (!psOptions->bSplitPolys)
    1031             :     {
    1032         116 :         auto poMP = std::make_unique<OGRMultiPolygon>();
    1033         116 :         for (auto &&poFeature : poMemLayer.get())
    1034             :         {
    1035             :             auto poGeom =
    1036         116 :                 std::unique_ptr<OGRGeometry>(poFeature->StealGeometry());
    1037          58 :             CPLAssert(poGeom);
    1038          58 :             if (poGeom->getGeometryType() == wkbPolygon)
    1039             :             {
    1040          58 :                 poMP->addGeometry(std::move(poGeom));
    1041             :             }
    1042             :         }
    1043          58 :         poMemLayer = std::make_unique<OGRMemLayer>("", nullptr, wkbUnknown);
    1044             :         auto poFeature =
    1045          58 :             std::make_unique<OGRFeature>(poMemLayer->GetLayerDefn());
    1046          58 :         poFeature->SetGeometryDirectly(poMP.release());
    1047          58 :         CPL_IGNORE_RET_VAL(poMemLayer->CreateFeature(poFeature.get()));
    1048             :     }
    1049             : 
    1050         122 :     for (auto &&poFeature : poMemLayer.get())
    1051             :     {
    1052          62 :         auto poGeom = std::unique_ptr<OGRGeometry>(poFeature->StealGeometry());
    1053          62 :         CPLAssert(poGeom);
    1054          62 :         if (poGeom->IsEmpty())
    1055           3 :             continue;
    1056             : 
    1057             :         auto poDstFeature =
    1058          59 :             std::make_unique<OGRFeature>(poDstLayer->GetLayerDefn());
    1059          59 :         poDstFeature->SetFrom(poFeature.get());
    1060             : 
    1061          59 :         if (poCT_GT)
    1062             :         {
    1063          28 :             if (poGeom->transform(poCT_GT.get()) != OGRERR_NONE)
    1064           0 :                 return false;
    1065             :         }
    1066             : 
    1067          59 :         if (psOptions->dfDensifyDistance > 0)
    1068             :         {
    1069           2 :             OGREnvelope sEnvelope;
    1070           2 :             poGeom->getEnvelope(&sEnvelope);
    1071             :             // Some sanity check to avoid insane memory allocations
    1072           2 :             if (sEnvelope.MaxX - sEnvelope.MinX >
    1073           2 :                     1e6 * psOptions->dfDensifyDistance ||
    1074           2 :                 sEnvelope.MaxY - sEnvelope.MinY >
    1075           2 :                     1e6 * psOptions->dfDensifyDistance)
    1076             :             {
    1077           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1078             :                          "Densification distance too small compared to "
    1079             :                          "geometry extent");
    1080           0 :                 return false;
    1081             :             }
    1082           2 :             poGeom->segmentize(psOptions->dfDensifyDistance);
    1083             :         }
    1084             : 
    1085          59 :         if (poCT_SRS)
    1086             :         {
    1087          21 :             if (poGeom->transform(poCT_SRS.get()) != OGRERR_NONE)
    1088           0 :                 return false;
    1089             :         }
    1090             : 
    1091          59 :         if (psOptions->dfMinRingArea != 0)
    1092             :         {
    1093           4 :             if (poGeom->getGeometryType() == wkbMultiPolygon)
    1094             :             {
    1095           8 :                 auto poMP = std::make_unique<OGRMultiPolygon>();
    1096           8 :                 for (auto *poPoly : poGeom->toMultiPolygon())
    1097             :                 {
    1098           8 :                     auto poNewPoly = std::make_unique<OGRPolygon>();
    1099           8 :                     for (auto *poRing : poPoly)
    1100             :                     {
    1101           4 :                         if (poRing->get_Area() >= psOptions->dfMinRingArea)
    1102             :                         {
    1103           2 :                             poNewPoly->addRing(poRing);
    1104             :                         }
    1105             :                     }
    1106           4 :                     if (!poNewPoly->IsEmpty())
    1107           2 :                         poMP->addGeometry(std::move(poNewPoly));
    1108             :                 }
    1109           4 :                 poGeom = std::move(poMP);
    1110             :             }
    1111           0 :             else if (poGeom->getGeometryType() == wkbPolygon)
    1112             :             {
    1113           0 :                 auto poNewPoly = std::make_unique<OGRPolygon>();
    1114           0 :                 for (auto *poRing : poGeom->toPolygon())
    1115             :                 {
    1116           0 :                     if (poRing->get_Area() >= psOptions->dfMinRingArea)
    1117             :                     {
    1118           0 :                         poNewPoly->addRing(poRing);
    1119             :                     }
    1120             :                 }
    1121           0 :                 poGeom = std::move(poNewPoly);
    1122             :             }
    1123           4 :             if (poGeom->IsEmpty())
    1124           2 :                 continue;
    1125             :         }
    1126             : 
    1127          57 :         if (psOptions->bConvexHull)
    1128             :         {
    1129           2 :             poGeom.reset(poGeom->ConvexHull());
    1130           2 :             if (!poGeom || poGeom->IsEmpty())
    1131           0 :                 continue;
    1132             :         }
    1133             : 
    1134          57 :         if (psOptions->dfSimplifyTolerance != 0)
    1135             :         {
    1136           2 :             poGeom.reset(poGeom->Simplify(psOptions->dfSimplifyTolerance));
    1137           2 :             if (!poGeom || poGeom->IsEmpty())
    1138           0 :                 continue;
    1139             :         }
    1140             : 
    1141         114 :         if (psOptions->nMaxPoints > 0 &&
    1142          57 :             CountPoints(poGeom.get()) >
    1143          57 :                 static_cast<size_t>(psOptions->nMaxPoints))
    1144             :         {
    1145           2 :             OGREnvelope sEnvelope;
    1146           2 :             poGeom->getEnvelope(&sEnvelope);
    1147           2 :             double tolMin = GetMinDistanceBetweenTwoPoints(poGeom.get());
    1148           4 :             double tolMax = std::max(sEnvelope.MaxY - sEnvelope.MinY,
    1149           2 :                                      sEnvelope.MaxX - sEnvelope.MinX);
    1150          42 :             for (int i = 0; i < 20; ++i)
    1151             :             {
    1152          40 :                 const double tol = (tolMin + tolMax) / 2;
    1153             :                 std::unique_ptr<OGRGeometry> poSimplifiedGeom(
    1154          40 :                     poGeom->Simplify(tol));
    1155          40 :                 if (!poSimplifiedGeom || poSimplifiedGeom->IsEmpty())
    1156             :                 {
    1157           5 :                     tolMax = tol;
    1158           5 :                     continue;
    1159             :                 }
    1160          35 :                 const auto nPoints = CountPoints(poSimplifiedGeom.get());
    1161          35 :                 if (nPoints == static_cast<size_t>(psOptions->nMaxPoints))
    1162             :                 {
    1163           0 :                     tolMax = tol;
    1164           0 :                     break;
    1165             :                 }
    1166          35 :                 else if (nPoints < static_cast<size_t>(psOptions->nMaxPoints))
    1167             :                 {
    1168          35 :                     tolMax = tol;
    1169             :                 }
    1170             :                 else
    1171             :                 {
    1172           0 :                     tolMin = tol;
    1173             :                 }
    1174             :             }
    1175           2 :             poGeom.reset(poGeom->Simplify(tolMax));
    1176           2 :             if (!poGeom || poGeom->IsEmpty())
    1177           0 :                 continue;
    1178             :         }
    1179             : 
    1180          57 :         if (!psOptions->bSplitPolys && poGeom->getGeometryType() == wkbPolygon)
    1181           6 :             poGeom.reset(
    1182             :                 OGRGeometryFactory::forceToMultiPolygon(poGeom.release()));
    1183             : 
    1184          57 :         poDstFeature->SetGeometryDirectly(poGeom.release());
    1185             : 
    1186          57 :         if (!psOptions->osLocationFieldName.empty())
    1187             :         {
    1188             : 
    1189         110 :             std::string osFilename = poSrcDS->GetDescription();
    1190             :             // Make sure it is a file before building absolute path name.
    1191             :             VSIStatBufL sStatBuf;
    1192         112 :             if (psOptions->bAbsolutePath &&
    1193          57 :                 CPLIsFilenameRelative(osFilename.c_str()) &&
    1194           2 :                 VSIStatL(osFilename.c_str(), &sStatBuf) == 0)
    1195             :             {
    1196           2 :                 char *pszCurDir = CPLGetCurrentDir();
    1197           2 :                 if (pszCurDir)
    1198             :                 {
    1199           4 :                     osFilename = CPLProjectRelativeFilenameSafe(
    1200           2 :                         pszCurDir, osFilename.c_str());
    1201           2 :                     CPLFree(pszCurDir);
    1202             :                 }
    1203             :             }
    1204          55 :             poDstFeature->SetField(psOptions->osLocationFieldName.c_str(),
    1205             :                                    osFilename.c_str());
    1206             :         }
    1207             : 
    1208          57 :         if (poDstLayer->CreateFeature(poDstFeature.get()) != OGRERR_NONE)
    1209             :         {
    1210           0 :             return false;
    1211             :         }
    1212             :     }
    1213             : 
    1214          60 :     return true;
    1215             : }
    1216             : 
    1217             : /************************************************************************/
    1218             : /*                  GDALFootprintAppGetParserUsage()                    */
    1219             : /************************************************************************/
    1220             : 
    1221           0 : std::string GDALFootprintAppGetParserUsage()
    1222             : {
    1223             :     try
    1224             :     {
    1225           0 :         GDALFootprintOptions sOptions;
    1226           0 :         GDALFootprintOptionsForBinary sOptionsForBinary;
    1227             :         auto argParser =
    1228           0 :             GDALFootprintAppOptionsGetParser(&sOptions, &sOptionsForBinary);
    1229           0 :         return argParser->usage();
    1230             :     }
    1231           0 :     catch (const std::exception &err)
    1232             :     {
    1233           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
    1234           0 :                  err.what());
    1235           0 :         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          78 : GDALDatasetH GDALFootprint(const char *pszDest, GDALDatasetH hDstDS,
    1269             :                            GDALDatasetH hSrcDataset,
    1270             :                            const GDALFootprintOptions *psOptionsIn,
    1271             :                            int *pbUsageError)
    1272             : {
    1273          78 :     if (pszDest == nullptr && hDstDS == nullptr)
    1274             :     {
    1275           1 :         CPLError(CE_Failure, CPLE_AppDefined,
    1276             :                  "pszDest == NULL && hDstDS == NULL");
    1277             : 
    1278           1 :         if (pbUsageError)
    1279           0 :             *pbUsageError = TRUE;
    1280           1 :         return nullptr;
    1281             :     }
    1282          77 :     if (hSrcDataset == nullptr)
    1283             :     {
    1284           1 :         CPLError(CE_Failure, CPLE_AppDefined, "hSrcDataset== NULL");
    1285             : 
    1286           1 :         if (pbUsageError)
    1287           0 :             *pbUsageError = TRUE;
    1288           1 :         return nullptr;
    1289             :     }
    1290          76 :     if (hDstDS != nullptr && psOptionsIn && psOptionsIn->bCreateOutput)
    1291             :     {
    1292           1 :         CPLError(CE_Failure, CPLE_AppDefined,
    1293             :                  "hDstDS != NULL but options that imply creating a new dataset "
    1294             :                  "have been set.");
    1295             : 
    1296           1 :         if (pbUsageError)
    1297           0 :             *pbUsageError = TRUE;
    1298           1 :         return nullptr;
    1299             :     }
    1300             : 
    1301          75 :     GDALFootprintOptions *psOptionsToFree = nullptr;
    1302          75 :     const GDALFootprintOptions *psOptions = psOptionsIn;
    1303          75 :     if (psOptions == nullptr)
    1304             :     {
    1305           1 :         psOptionsToFree = GDALFootprintOptionsNew(nullptr, nullptr);
    1306           1 :         psOptions = psOptionsToFree;
    1307             :     }
    1308             : 
    1309          75 :     const bool bCloseOutDSOnError = hDstDS == nullptr;
    1310             : 
    1311          75 :     auto poSrcDS = GDALDataset::FromHandle(hSrcDataset);
    1312          75 :     if (poSrcDS->GetRasterCount() == 0)
    1313             :     {
    1314           1 :         CPLError(CE_Failure, CPLE_AppDefined,
    1315             :                  "Input dataset has no raster band.%s",
    1316           1 :                  poSrcDS->GetMetadata("SUBDATASETS")
    1317             :                      ? " You need to specify one subdataset."
    1318             :                      : "");
    1319           1 :         GDALFootprintOptionsFree(psOptionsToFree);
    1320           1 :         if (bCloseOutDSOnError)
    1321           0 :             GDALClose(hDstDS);
    1322           1 :         return nullptr;
    1323             :     }
    1324             : 
    1325             :     auto poLayer =
    1326          74 :         GetOutputLayerAndUpdateDstDS(pszDest, hDstDS, poSrcDS, psOptions);
    1327          74 :     if (!poLayer)
    1328             :     {
    1329           5 :         GDALFootprintOptionsFree(psOptionsToFree);
    1330           5 :         if (hDstDS && bCloseOutDSOnError)
    1331           0 :             GDALClose(hDstDS);
    1332           5 :         return nullptr;
    1333             :     }
    1334             : 
    1335          69 :     if (!GDALFootprintProcess(poSrcDS, poLayer, psOptions))
    1336             :     {
    1337           9 :         GDALFootprintOptionsFree(psOptionsToFree);
    1338           9 :         if (bCloseOutDSOnError)
    1339           8 :             GDALClose(hDstDS);
    1340           9 :         return nullptr;
    1341             :     }
    1342             : 
    1343          60 :     GDALFootprintOptionsFree(psOptionsToFree);
    1344             : 
    1345          60 :     return hDstDS;
    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          79 : GDALFootprintOptionsNew(char **papszArgv,
    1370             :                         GDALFootprintOptionsForBinary *psOptionsForBinary)
    1371             : {
    1372         158 :     auto psOptions = std::make_unique<GDALFootprintOptions>();
    1373             : 
    1374             :     /* -------------------------------------------------------------------- */
    1375             :     /*      Parse arguments.                                                */
    1376             :     /* -------------------------------------------------------------------- */
    1377             : 
    1378         158 :     CPLStringList aosArgv;
    1379             : 
    1380          79 :     if (papszArgv)
    1381             :     {
    1382          78 :         const int nArgc = CSLCount(papszArgv);
    1383         682 :         for (int i = 0; i < nArgc; i++)
    1384             :         {
    1385         604 :             aosArgv.AddString(papszArgv[i]);
    1386             :         }
    1387             :     }
    1388             : 
    1389             :     try
    1390             :     {
    1391             :         auto argParser = GDALFootprintAppOptionsGetParser(psOptions.get(),
    1392          80 :                                                           psOptionsForBinary);
    1393             : 
    1394          79 :         argParser->parse_args_without_binary_name(aosArgv.List());
    1395             : 
    1396          78 :         if (argParser->is_used("-t_srs"))
    1397             :         {
    1398             : 
    1399           4 :             const std::string osVal(argParser->get<std::string>("-t_srs"));
    1400           4 :             if (psOptions->oOutputSRS.SetFromUserInput(osVal.c_str()) !=
    1401             :                 OGRERR_NONE)
    1402             :             {
    1403           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
    1404             :                          "Failed to process SRS definition: %s", osVal.c_str());
    1405           0 :                 return nullptr;
    1406             :             }
    1407           4 :             psOptions->oOutputSRS.SetAxisMappingStrategy(
    1408             :                 OAMS_TRADITIONAL_GIS_ORDER);
    1409             :         }
    1410             : 
    1411          78 :         if (argParser->is_used("-max_points"))
    1412             :         {
    1413             :             const std::string maxPoints{
    1414          32 :                 argParser->get<std::string>("-max_points")};
    1415          32 :             if (maxPoints == "unlimited")
    1416             :             {
    1417           0 :                 psOptions->nMaxPoints = 0;
    1418             :             }
    1419             :             else
    1420             :             {
    1421          32 :                 psOptions->nMaxPoints = atoi(maxPoints.c_str());
    1422          32 :                 if (psOptions->nMaxPoints > 0 && psOptions->nMaxPoints < 3)
    1423             :                 {
    1424           0 :                     CPLError(CE_Failure, CPLE_NotSupported,
    1425             :                              "Invalid value for -max_points");
    1426           0 :                     return nullptr;
    1427             :                 }
    1428             :             }
    1429             :         }
    1430             : 
    1431          78 :         psOptions->bCreateOutput = !psOptions->osFormat.empty();
    1432             :     }
    1433           1 :     catch (const std::exception &err)
    1434             :     {
    1435           1 :         CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
    1436           1 :                  err.what());
    1437           1 :         return nullptr;
    1438             :     }
    1439             : 
    1440          78 :     if (!psOptions->bOutCSGeoref && !psOptions->oOutputSRS.IsEmpty())
    1441             :     {
    1442           1 :         CPLError(CE_Failure, CPLE_AppDefined,
    1443             :                  "-t_cs pixel and -t_srs are mutually exclusive.");
    1444           1 :         return nullptr;
    1445             :     }
    1446             : 
    1447          77 :     if (psOptions->bClearLocation)
    1448             :     {
    1449           2 :         psOptions->osLocationFieldName.clear();
    1450             :     }
    1451             : 
    1452          77 :     if (psOptionsForBinary)
    1453             :     {
    1454           7 :         psOptionsForBinary->bCreateOutput = psOptions->bCreateOutput;
    1455           7 :         psOptionsForBinary->osFormat = psOptions->osFormat;
    1456           7 :         psOptionsForBinary->osDestLayerName = psOptions->osDestLayerName;
    1457             :     }
    1458             : 
    1459          77 :     return psOptions.release();
    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         150 : void GDALFootprintOptionsFree(GDALFootprintOptions *psOptions)
    1475             : {
    1476         150 :     delete psOptions;
    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          37 : void GDALFootprintOptionsSetProgress(GDALFootprintOptions *psOptions,
    1494             :                                      GDALProgressFunc pfnProgress,
    1495             :                                      void *pProgressData)
    1496             : {
    1497          37 :     psOptions->pfnProgress = pfnProgress ? pfnProgress : GDALDummyProgress;
    1498          37 :     psOptions->pProgressData = pProgressData;
    1499          37 : }

Generated by: LCOV version 1.14