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