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