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_Byte;
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_Byte, 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_Byte && 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_Byte;
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_Byte, 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_Byte && 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_Byte,
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_Byte,
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 : char **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 0 : std::unique_ptr<OGRSpatialReference, OGRSpatialReferenceReleaser> poSRS;
664 64 : if (psOptions->bOutCSGeoref)
665 : {
666 48 : if (!psOptions->oOutputSRS.IsEmpty())
667 : {
668 3 : poSRS.reset(psOptions->oOutputSRS.Clone());
669 : }
670 45 : else if (auto poSrcSRS = poSrcDS->GetSpatialRef())
671 : {
672 17 : poSRS.reset(poSrcSRS->Clone());
673 : }
674 : }
675 :
676 192 : poLayer = poDstDS->CreateLayer(
677 64 : osDestLayerName.c_str(), poSRS.get(),
678 64 : psOptions->bSplitPolys ? wkbPolygon : wkbMultiPolygon,
679 : const_cast<char **>(psOptions->aosLCO.List()));
680 :
681 64 : if (!psOptions->osLocationFieldName.empty())
682 : {
683 : OGRFieldDefn oFieldDefn(psOptions->osLocationFieldName.c_str(),
684 62 : OFTString);
685 62 : if (poLayer->CreateField(&oFieldDefn) != OGRERR_NONE)
686 0 : return nullptr;
687 : }
688 : }
689 :
690 69 : return poLayer;
691 : }
692 :
693 : /************************************************************************/
694 : /* GeoTransformCoordinateTransformation */
695 : /************************************************************************/
696 :
697 : class GeoTransformCoordinateTransformation final
698 : : public OGRCoordinateTransformation
699 : {
700 : const GDALGeoTransform &m_gt;
701 :
702 : public:
703 28 : explicit GeoTransformCoordinateTransformation(const GDALGeoTransform >)
704 28 : : m_gt(gt)
705 : {
706 28 : }
707 :
708 : const OGRSpatialReference *GetSourceCS() const override;
709 :
710 84 : const OGRSpatialReference *GetTargetCS() const override
711 : {
712 84 : return nullptr;
713 : }
714 :
715 0 : OGRCoordinateTransformation *Clone() const override
716 : {
717 0 : return new GeoTransformCoordinateTransformation(m_gt);
718 : }
719 :
720 0 : OGRCoordinateTransformation *GetInverse() const override
721 : {
722 0 : CPLError(CE_Failure, CPLE_NotSupported,
723 : "GeoTransformCoordinateTransformation::GetInverse() not "
724 : "implemented");
725 0 : return nullptr;
726 : }
727 :
728 28 : int Transform(size_t nCount, double *x, double *y, double * /* z */,
729 : double * /* t */, int *pabSuccess) override
730 : {
731 168 : for (size_t i = 0; i < nCount; ++i)
732 : {
733 140 : const double X = m_gt[0] + x[i] * m_gt[1] + y[i] * m_gt[2];
734 140 : const double Y = m_gt[3] + x[i] * m_gt[4] + y[i] * m_gt[5];
735 140 : x[i] = X;
736 140 : y[i] = Y;
737 140 : if (pabSuccess)
738 140 : pabSuccess[i] = TRUE;
739 : }
740 28 : return TRUE;
741 : }
742 : };
743 :
744 : const OGRSpatialReference *
745 0 : GeoTransformCoordinateTransformation::GetSourceCS() const
746 : {
747 0 : return nullptr;
748 : }
749 :
750 : /************************************************************************/
751 : /* CountPoints() */
752 : /************************************************************************/
753 :
754 142 : static size_t CountPoints(const OGRGeometry *poGeom)
755 : {
756 142 : if (poGeom->getGeometryType() == wkbMultiPolygon)
757 : {
758 49 : size_t n = 0;
759 99 : for (auto *poPoly : poGeom->toMultiPolygon())
760 : {
761 50 : n += CountPoints(poPoly);
762 : }
763 49 : return n;
764 : }
765 93 : else if (poGeom->getGeometryType() == wkbPolygon)
766 : {
767 93 : size_t n = 0;
768 187 : for (auto *poRing : poGeom->toPolygon())
769 : {
770 94 : n += poRing->getNumPoints() - 1;
771 : }
772 93 : return n;
773 : }
774 0 : return 0;
775 : }
776 :
777 : /************************************************************************/
778 : /* GetMinDistanceBetweenTwoPoints() */
779 : /************************************************************************/
780 :
781 6 : static double GetMinDistanceBetweenTwoPoints(const OGRGeometry *poGeom)
782 : {
783 6 : if (poGeom->getGeometryType() == wkbMultiPolygon)
784 : {
785 2 : double v = std::numeric_limits<double>::max();
786 4 : for (auto *poPoly : poGeom->toMultiPolygon())
787 : {
788 2 : v = std::min(v, GetMinDistanceBetweenTwoPoints(poPoly));
789 : }
790 2 : return v;
791 : }
792 4 : else if (poGeom->getGeometryType() == wkbPolygon)
793 : {
794 2 : double v = std::numeric_limits<double>::max();
795 4 : for (auto *poRing : poGeom->toPolygon())
796 : {
797 2 : v = std::min(v, GetMinDistanceBetweenTwoPoints(poRing));
798 : }
799 2 : return v;
800 : }
801 2 : else if (poGeom->getGeometryType() == wkbLineString)
802 : {
803 2 : double v = std::numeric_limits<double>::max();
804 2 : const auto poLS = poGeom->toLineString();
805 2 : const int nNumPoints = poLS->getNumPoints();
806 44 : for (int i = 0; i < nNumPoints - 1; ++i)
807 : {
808 42 : const double x1 = poLS->getX(i);
809 42 : const double y1 = poLS->getY(i);
810 42 : const double x2 = poLS->getX(i + 1);
811 42 : const double y2 = poLS->getY(i + 1);
812 42 : const double d = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);
813 42 : if (d > 0)
814 42 : v = std::min(v, d);
815 : }
816 2 : return sqrt(v);
817 : }
818 0 : return 0;
819 : }
820 :
821 : /************************************************************************/
822 : /* GDALFootprintProcess() */
823 : /************************************************************************/
824 :
825 69 : static bool GDALFootprintProcess(GDALDataset *poSrcDS, OGRLayer *poDstLayer,
826 : const GDALFootprintOptions *psOptions)
827 : {
828 69 : std::unique_ptr<OGRCoordinateTransformation> poCT_SRS;
829 69 : const OGRSpatialReference *poDstSRS = poDstLayer->GetSpatialRef();
830 69 : if (!psOptions->oOutputSRS.IsEmpty())
831 3 : poDstSRS = &(psOptions->oOutputSRS);
832 69 : if (poDstSRS)
833 : {
834 22 : auto poSrcSRS = poSrcDS->GetSpatialRef();
835 22 : if (!poSrcSRS)
836 : {
837 1 : CPLError(CE_Failure, CPLE_AppDefined,
838 : "Output layer has CRS, but input is not georeferenced");
839 1 : return false;
840 : }
841 21 : poCT_SRS.reset(OGRCreateCoordinateTransformation(poSrcSRS, poDstSRS));
842 21 : if (!poCT_SRS)
843 0 : return false;
844 : }
845 :
846 136 : std::vector<int> anBands = psOptions->anBands;
847 68 : const int nBandCount = poSrcDS->GetRasterCount();
848 68 : if (anBands.empty())
849 : {
850 139 : for (int i = 1; i <= nBandCount; ++i)
851 76 : anBands.push_back(i);
852 : }
853 :
854 136 : std::vector<GDALRasterBand *> apoSrcMaskBands;
855 : const CPLStringList aosSrcNoData(
856 136 : CSLTokenizeString2(psOptions->osSrcNoData.c_str(), " ", 0));
857 136 : std::vector<double> adfSrcNoData;
858 68 : if (!psOptions->osSrcNoData.empty())
859 : {
860 6 : if (aosSrcNoData.size() != 1 &&
861 2 : static_cast<size_t>(aosSrcNoData.size()) != anBands.size())
862 : {
863 1 : CPLError(CE_Failure, CPLE_AppDefined,
864 : "Number of values in -srcnodata should be 1 or the number "
865 : "of bands");
866 1 : return false;
867 : }
868 7 : for (int i = 0; i < aosSrcNoData.size(); ++i)
869 : {
870 4 : adfSrcNoData.emplace_back(CPLAtof(aosSrcNoData[i]));
871 : }
872 : }
873 67 : bool bGlobalMask = true;
874 134 : std::vector<std::unique_ptr<GDALRasterBand>> apoTmpNoDataMaskBands;
875 142 : for (size_t i = 0; i < anBands.size(); ++i)
876 : {
877 80 : const int nBand = anBands[i];
878 80 : if (nBand <= 0 || nBand > nBandCount)
879 : {
880 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid band number: %d",
881 : nBand);
882 5 : return false;
883 : }
884 79 : auto poBand = poSrcDS->GetRasterBand(nBand);
885 79 : if (!adfSrcNoData.empty())
886 : {
887 4 : bGlobalMask = false;
888 : apoTmpNoDataMaskBands.emplace_back(
889 12 : std::make_unique<GDALNoDataMaskBand>(
890 6 : poBand, adfSrcNoData.size() == 1 ? adfSrcNoData[0]
891 6 : : adfSrcNoData[i]));
892 4 : apoSrcMaskBands.push_back(apoTmpNoDataMaskBands.back().get());
893 : }
894 : else
895 : {
896 : GDALRasterBand *poMaskBand;
897 75 : const int nMaskFlags = poBand->GetMaskFlags();
898 75 : if (poBand->GetColorInterpretation() == GCI_AlphaBand)
899 : {
900 2 : poMaskBand = poBand;
901 : }
902 : else
903 : {
904 73 : if ((nMaskFlags & GMF_PER_DATASET) == 0)
905 : {
906 65 : bGlobalMask = false;
907 : }
908 73 : poMaskBand = poBand->GetMaskBand();
909 : }
910 75 : if (psOptions->nOvrIndex >= 0)
911 : {
912 11 : if (nMaskFlags == GMF_NODATA)
913 : {
914 : // If the mask band is based on nodata, we don't need
915 : // to check the overviews of the mask band, but we
916 : // can take the mask band of the overviews
917 5 : auto poOvrBand = poBand->GetOverview(psOptions->nOvrIndex);
918 5 : if (!poOvrBand)
919 : {
920 2 : if (poBand->GetOverviewCount() == 0)
921 : {
922 1 : CPLError(
923 : CE_Failure, CPLE_AppDefined,
924 : "Overview index %d invalid for this dataset. "
925 : "Bands of this dataset have no "
926 : "precomputed overviews",
927 1 : psOptions->nOvrIndex);
928 : }
929 : else
930 : {
931 1 : CPLError(
932 : CE_Failure, CPLE_AppDefined,
933 : "Overview index %d invalid for this dataset. "
934 : "Value should be in [0,%d] range",
935 1 : psOptions->nOvrIndex,
936 1 : poBand->GetOverviewCount() - 1);
937 : }
938 4 : return false;
939 : }
940 3 : if (poOvrBand->GetMaskFlags() != GMF_NODATA)
941 : {
942 0 : CPLError(CE_Failure, CPLE_AppDefined,
943 : "poOvrBand->GetMaskFlags() != GMF_NODATA");
944 0 : return false;
945 : }
946 3 : poMaskBand = poOvrBand->GetMaskBand();
947 : }
948 : else
949 : {
950 6 : poMaskBand = poMaskBand->GetOverview(psOptions->nOvrIndex);
951 6 : if (!poMaskBand)
952 : {
953 2 : if (poBand->GetMaskBand()->GetOverviewCount() == 0)
954 : {
955 1 : CPLError(
956 : CE_Failure, CPLE_AppDefined,
957 : "Overview index %d invalid for this dataset. "
958 : "Mask bands of this dataset have no "
959 : "precomputed overviews",
960 1 : psOptions->nOvrIndex);
961 : }
962 : else
963 : {
964 1 : CPLError(
965 : CE_Failure, CPLE_AppDefined,
966 : "Overview index %d invalid for this dataset. "
967 : "Value should be in [0,%d] range",
968 1 : psOptions->nOvrIndex,
969 1 : poBand->GetMaskBand()->GetOverviewCount() - 1);
970 : }
971 2 : return false;
972 : }
973 : }
974 : }
975 71 : apoSrcMaskBands.push_back(poMaskBand);
976 : }
977 : }
978 :
979 62 : std::unique_ptr<OGRCoordinateTransformation> poCT_GT;
980 62 : GDALGeoTransform gt;
981 62 : if (psOptions->bOutCSGeoref && poSrcDS->GetGeoTransform(gt) == CE_None)
982 : {
983 25 : auto poMaskBand = apoSrcMaskBands[0];
984 25 : gt.Rescale(double(poSrcDS->GetRasterXSize()) / poMaskBand->GetXSize(),
985 25 : double(poSrcDS->GetRasterYSize()) / poMaskBand->GetYSize());
986 25 : poCT_GT = std::make_unique<GeoTransformCoordinateTransformation>(gt);
987 : }
988 37 : else if (psOptions->bOutCSGeorefRequested)
989 : {
990 2 : CPLError(CE_Failure, CPLE_AppDefined,
991 : "Georeferenced coordinates requested, but "
992 : "input dataset has no geotransform.");
993 2 : return false;
994 : }
995 35 : else if (psOptions->nOvrIndex >= 0)
996 : {
997 : // Transform from overview pixel coordinates to full resolution
998 : // pixel coordinates
999 3 : auto poMaskBand = apoSrcMaskBands[0];
1000 3 : gt[1] = double(poSrcDS->GetRasterXSize()) / poMaskBand->GetXSize();
1001 3 : gt[2] = 0;
1002 3 : gt[4] = 0;
1003 3 : gt[5] = double(poSrcDS->GetRasterYSize()) / poMaskBand->GetYSize();
1004 3 : poCT_GT = std::make_unique<GeoTransformCoordinateTransformation>(gt);
1005 : }
1006 :
1007 60 : std::unique_ptr<GDALRasterBand> poMaskForRasterize;
1008 60 : if (bGlobalMask || anBands.size() == 1)
1009 : {
1010 : poMaskForRasterize =
1011 53 : std::make_unique<GDALFootprintMaskBand>(apoSrcMaskBands[0]);
1012 : }
1013 : else
1014 : {
1015 7 : poMaskForRasterize = std::make_unique<GDALFootprintCombinedMaskBand>(
1016 7 : apoSrcMaskBands, psOptions->bCombineBandsUnion);
1017 : }
1018 :
1019 60 : auto hBand = GDALRasterBand::ToHandle(poMaskForRasterize.get());
1020 120 : auto poMemLayer = std::make_unique<OGRMemLayer>("", nullptr, wkbUnknown);
1021 : const CPLErr eErr =
1022 60 : GDALPolygonize(hBand, hBand, OGRLayer::ToHandle(poMemLayer.get()),
1023 : /* iPixValField = */ -1,
1024 60 : /* papszOptions = */ nullptr, psOptions->pfnProgress,
1025 60 : psOptions->pProgressData);
1026 60 : if (eErr != CE_None)
1027 : {
1028 0 : return false;
1029 : }
1030 :
1031 60 : if (!psOptions->bSplitPolys)
1032 : {
1033 116 : auto poMP = std::make_unique<OGRMultiPolygon>();
1034 116 : for (auto &&poFeature : poMemLayer.get())
1035 : {
1036 : auto poGeom =
1037 116 : std::unique_ptr<OGRGeometry>(poFeature->StealGeometry());
1038 58 : CPLAssert(poGeom);
1039 58 : if (poGeom->getGeometryType() == wkbPolygon)
1040 : {
1041 58 : poMP->addGeometry(std::move(poGeom));
1042 : }
1043 : }
1044 58 : poMemLayer = std::make_unique<OGRMemLayer>("", nullptr, wkbUnknown);
1045 : auto poFeature =
1046 58 : std::make_unique<OGRFeature>(poMemLayer->GetLayerDefn());
1047 58 : poFeature->SetGeometryDirectly(poMP.release());
1048 58 : CPL_IGNORE_RET_VAL(poMemLayer->CreateFeature(poFeature.get()));
1049 : }
1050 :
1051 122 : for (auto &&poFeature : poMemLayer.get())
1052 : {
1053 62 : auto poGeom = std::unique_ptr<OGRGeometry>(poFeature->StealGeometry());
1054 62 : CPLAssert(poGeom);
1055 62 : if (poGeom->IsEmpty())
1056 3 : continue;
1057 :
1058 : auto poDstFeature =
1059 59 : std::make_unique<OGRFeature>(poDstLayer->GetLayerDefn());
1060 59 : poDstFeature->SetFrom(poFeature.get());
1061 :
1062 59 : if (poCT_GT)
1063 : {
1064 28 : if (poGeom->transform(poCT_GT.get()) != OGRERR_NONE)
1065 0 : return false;
1066 : }
1067 :
1068 59 : if (psOptions->dfDensifyDistance > 0)
1069 : {
1070 2 : OGREnvelope sEnvelope;
1071 2 : poGeom->getEnvelope(&sEnvelope);
1072 : // Some sanity check to avoid insane memory allocations
1073 2 : if (sEnvelope.MaxX - sEnvelope.MinX >
1074 2 : 1e6 * psOptions->dfDensifyDistance ||
1075 2 : sEnvelope.MaxY - sEnvelope.MinY >
1076 2 : 1e6 * psOptions->dfDensifyDistance)
1077 : {
1078 0 : CPLError(CE_Failure, CPLE_AppDefined,
1079 : "Densification distance too small compared to "
1080 : "geometry extent");
1081 0 : return false;
1082 : }
1083 2 : poGeom->segmentize(psOptions->dfDensifyDistance);
1084 : }
1085 :
1086 59 : if (poCT_SRS)
1087 : {
1088 21 : if (poGeom->transform(poCT_SRS.get()) != OGRERR_NONE)
1089 0 : return false;
1090 : }
1091 :
1092 59 : if (psOptions->dfMinRingArea != 0)
1093 : {
1094 4 : if (poGeom->getGeometryType() == wkbMultiPolygon)
1095 : {
1096 8 : auto poMP = std::make_unique<OGRMultiPolygon>();
1097 8 : for (auto *poPoly : poGeom->toMultiPolygon())
1098 : {
1099 8 : auto poNewPoly = std::make_unique<OGRPolygon>();
1100 8 : for (auto *poRing : poPoly)
1101 : {
1102 4 : if (poRing->get_Area() >= psOptions->dfMinRingArea)
1103 : {
1104 2 : poNewPoly->addRing(poRing);
1105 : }
1106 : }
1107 4 : if (!poNewPoly->IsEmpty())
1108 2 : poMP->addGeometry(std::move(poNewPoly));
1109 : }
1110 4 : poGeom = std::move(poMP);
1111 : }
1112 0 : else if (poGeom->getGeometryType() == wkbPolygon)
1113 : {
1114 0 : auto poNewPoly = std::make_unique<OGRPolygon>();
1115 0 : for (auto *poRing : poGeom->toPolygon())
1116 : {
1117 0 : if (poRing->get_Area() >= psOptions->dfMinRingArea)
1118 : {
1119 0 : poNewPoly->addRing(poRing);
1120 : }
1121 : }
1122 0 : poGeom = std::move(poNewPoly);
1123 : }
1124 4 : if (poGeom->IsEmpty())
1125 2 : continue;
1126 : }
1127 :
1128 12 : const auto GeomIsNullOrEmpty = [&poGeom]
1129 6 : { return !poGeom || poGeom->IsEmpty(); };
1130 :
1131 57 : if (psOptions->bConvexHull)
1132 : {
1133 2 : poGeom.reset(poGeom->ConvexHull());
1134 2 : if (GeomIsNullOrEmpty())
1135 0 : continue;
1136 : }
1137 :
1138 16 : const auto DoSimplification = [&poMemLayer, &poGeom](double dfTolerance)
1139 : {
1140 8 : std::string osLastErrorMsg;
1141 : {
1142 8 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
1143 4 : CPLErrorReset();
1144 4 : poGeom.reset(poGeom->Simplify(dfTolerance));
1145 4 : osLastErrorMsg = CPLGetLastErrorMsg();
1146 : }
1147 4 : if (!poGeom || poGeom->IsEmpty())
1148 : {
1149 0 : if (poMemLayer->GetFeatureCount(false) == 1)
1150 : {
1151 0 : if (!osLastErrorMsg.empty())
1152 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s",
1153 : osLastErrorMsg.c_str());
1154 : else
1155 0 : CPLError(CE_Failure, CPLE_AppDefined,
1156 : "Simplification resulted in empty geometry");
1157 0 : return false;
1158 : }
1159 0 : if (!osLastErrorMsg.empty())
1160 0 : CPLError(CE_Warning, CPLE_AppDefined, "%s",
1161 : osLastErrorMsg.c_str());
1162 : }
1163 4 : return true;
1164 57 : };
1165 :
1166 57 : if (psOptions->dfSimplifyTolerance != 0)
1167 : {
1168 2 : if (!DoSimplification(psOptions->dfSimplifyTolerance))
1169 0 : return false;
1170 2 : if (GeomIsNullOrEmpty())
1171 0 : continue;
1172 : }
1173 :
1174 114 : if (psOptions->nMaxPoints > 0 &&
1175 57 : CountPoints(poGeom.get()) >
1176 57 : static_cast<size_t>(psOptions->nMaxPoints))
1177 : {
1178 2 : OGREnvelope sEnvelope;
1179 2 : poGeom->getEnvelope(&sEnvelope);
1180 2 : double tolMin = GetMinDistanceBetweenTwoPoints(poGeom.get());
1181 4 : double tolMax = std::max(sEnvelope.MaxY - sEnvelope.MinY,
1182 2 : sEnvelope.MaxX - sEnvelope.MinX);
1183 42 : for (int i = 0; i < 20; ++i)
1184 : {
1185 40 : const double tol = (tolMin + tolMax) / 2;
1186 : std::unique_ptr<OGRGeometry> poSimplifiedGeom(
1187 40 : poGeom->Simplify(tol));
1188 40 : if (!poSimplifiedGeom || poSimplifiedGeom->IsEmpty())
1189 : {
1190 5 : tolMax = tol;
1191 5 : continue;
1192 : }
1193 35 : const auto nPoints = CountPoints(poSimplifiedGeom.get());
1194 35 : if (nPoints == static_cast<size_t>(psOptions->nMaxPoints))
1195 : {
1196 0 : tolMax = tol;
1197 0 : break;
1198 : }
1199 35 : else if (nPoints < static_cast<size_t>(psOptions->nMaxPoints))
1200 : {
1201 35 : tolMax = tol;
1202 : }
1203 : else
1204 : {
1205 0 : tolMin = tol;
1206 : }
1207 : }
1208 :
1209 2 : if (!DoSimplification(tolMax))
1210 0 : return false;
1211 2 : if (GeomIsNullOrEmpty())
1212 0 : continue;
1213 : }
1214 :
1215 57 : if (!psOptions->bSplitPolys && poGeom->getGeometryType() == wkbPolygon)
1216 6 : poGeom.reset(
1217 : OGRGeometryFactory::forceToMultiPolygon(poGeom.release()));
1218 :
1219 57 : poDstFeature->SetGeometryDirectly(poGeom.release());
1220 :
1221 57 : if (!psOptions->osLocationFieldName.empty())
1222 : {
1223 :
1224 110 : std::string osFilename = poSrcDS->GetDescription();
1225 : // Make sure it is a file before building absolute path name.
1226 : VSIStatBufL sStatBuf;
1227 112 : if (psOptions->bAbsolutePath &&
1228 57 : CPLIsFilenameRelative(osFilename.c_str()) &&
1229 2 : VSIStatL(osFilename.c_str(), &sStatBuf) == 0)
1230 : {
1231 2 : char *pszCurDir = CPLGetCurrentDir();
1232 2 : if (pszCurDir)
1233 : {
1234 4 : osFilename = CPLProjectRelativeFilenameSafe(
1235 2 : pszCurDir, osFilename.c_str());
1236 2 : CPLFree(pszCurDir);
1237 : }
1238 : }
1239 55 : poDstFeature->SetField(psOptions->osLocationFieldName.c_str(),
1240 : osFilename.c_str());
1241 : }
1242 :
1243 57 : if (poDstLayer->CreateFeature(poDstFeature.get()) != OGRERR_NONE)
1244 : {
1245 0 : return false;
1246 : }
1247 : }
1248 :
1249 60 : return true;
1250 : }
1251 :
1252 : /************************************************************************/
1253 : /* GDALFootprintAppGetParserUsage() */
1254 : /************************************************************************/
1255 :
1256 0 : std::string GDALFootprintAppGetParserUsage()
1257 : {
1258 : try
1259 : {
1260 0 : GDALFootprintOptions sOptions;
1261 0 : GDALFootprintOptionsForBinary sOptionsForBinary;
1262 : auto argParser =
1263 0 : GDALFootprintAppOptionsGetParser(&sOptions, &sOptionsForBinary);
1264 0 : return argParser->usage();
1265 : }
1266 0 : catch (const std::exception &err)
1267 : {
1268 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
1269 0 : err.what());
1270 0 : return std::string();
1271 : }
1272 : }
1273 :
1274 : /************************************************************************/
1275 : /* GDALFootprint() */
1276 : /************************************************************************/
1277 :
1278 : /* clang-format off */
1279 : /**
1280 : * Computes the footprint of a raster.
1281 : *
1282 : * This is the equivalent of the
1283 : * <a href="/programs/gdal_footprint.html">gdal_footprint</a> utility.
1284 : *
1285 : * GDALFootprintOptions* must be allocated and freed with
1286 : * GDALFootprintOptionsNew() and GDALFootprintOptionsFree() respectively.
1287 : * pszDest and hDstDS cannot be used at the same time.
1288 : *
1289 : * @param pszDest the vector destination dataset path or NULL.
1290 : * @param hDstDS the vector destination dataset or NULL.
1291 : * @param hSrcDataset the raster source dataset handle.
1292 : * @param psOptionsIn the options struct returned by GDALFootprintOptionsNew()
1293 : * or NULL.
1294 : * @param pbUsageError pointer to a integer output variable to store if any
1295 : * usage error has occurred or NULL.
1296 : * @return the output dataset (new dataset that must be closed using
1297 : * GDALClose(), or hDstDS is not NULL) or NULL in case of error.
1298 : *
1299 : * @since GDAL 3.8
1300 : */
1301 : /* clang-format on */
1302 :
1303 78 : GDALDatasetH GDALFootprint(const char *pszDest, GDALDatasetH hDstDS,
1304 : GDALDatasetH hSrcDataset,
1305 : const GDALFootprintOptions *psOptionsIn,
1306 : int *pbUsageError)
1307 : {
1308 78 : if (pszDest == nullptr && hDstDS == nullptr)
1309 : {
1310 1 : CPLError(CE_Failure, CPLE_AppDefined,
1311 : "pszDest == NULL && hDstDS == NULL");
1312 :
1313 1 : if (pbUsageError)
1314 0 : *pbUsageError = TRUE;
1315 1 : return nullptr;
1316 : }
1317 77 : if (hSrcDataset == nullptr)
1318 : {
1319 1 : CPLError(CE_Failure, CPLE_AppDefined, "hSrcDataset== NULL");
1320 :
1321 1 : if (pbUsageError)
1322 0 : *pbUsageError = TRUE;
1323 1 : return nullptr;
1324 : }
1325 76 : if (hDstDS != nullptr && psOptionsIn && psOptionsIn->bCreateOutput)
1326 : {
1327 1 : CPLError(CE_Failure, CPLE_AppDefined,
1328 : "hDstDS != NULL but options that imply creating a new dataset "
1329 : "have been set.");
1330 :
1331 1 : if (pbUsageError)
1332 0 : *pbUsageError = TRUE;
1333 1 : return nullptr;
1334 : }
1335 :
1336 75 : GDALFootprintOptions *psOptionsToFree = nullptr;
1337 75 : const GDALFootprintOptions *psOptions = psOptionsIn;
1338 75 : if (psOptions == nullptr)
1339 : {
1340 1 : psOptionsToFree = GDALFootprintOptionsNew(nullptr, nullptr);
1341 1 : psOptions = psOptionsToFree;
1342 : }
1343 :
1344 75 : const bool bCloseOutDSOnError = hDstDS == nullptr;
1345 :
1346 75 : auto poSrcDS = GDALDataset::FromHandle(hSrcDataset);
1347 75 : if (poSrcDS->GetRasterCount() == 0)
1348 : {
1349 1 : CPLError(CE_Failure, CPLE_AppDefined,
1350 : "Input dataset has no raster band.%s",
1351 1 : poSrcDS->GetMetadata("SUBDATASETS")
1352 : ? " You need to specify one subdataset."
1353 : : "");
1354 1 : GDALFootprintOptionsFree(psOptionsToFree);
1355 1 : if (bCloseOutDSOnError)
1356 0 : GDALClose(hDstDS);
1357 1 : return nullptr;
1358 : }
1359 :
1360 : auto poLayer =
1361 74 : GetOutputLayerAndUpdateDstDS(pszDest, hDstDS, poSrcDS, psOptions);
1362 74 : if (!poLayer)
1363 : {
1364 5 : GDALFootprintOptionsFree(psOptionsToFree);
1365 5 : if (hDstDS && bCloseOutDSOnError)
1366 0 : GDALClose(hDstDS);
1367 5 : return nullptr;
1368 : }
1369 :
1370 69 : if (!GDALFootprintProcess(poSrcDS, poLayer, psOptions))
1371 : {
1372 9 : GDALFootprintOptionsFree(psOptionsToFree);
1373 9 : if (bCloseOutDSOnError)
1374 8 : GDALClose(hDstDS);
1375 9 : return nullptr;
1376 : }
1377 :
1378 60 : GDALFootprintOptionsFree(psOptionsToFree);
1379 :
1380 60 : return hDstDS;
1381 : }
1382 :
1383 : /************************************************************************/
1384 : /* GDALFootprintOptionsNew() */
1385 : /************************************************************************/
1386 :
1387 : /**
1388 : * Allocates a GDALFootprintOptions struct.
1389 : *
1390 : * @param papszArgv NULL terminated list of options (potentially including
1391 : * filename and open options too), or NULL. The accepted options are the ones of
1392 : * the <a href="/programs/gdal_rasterize.html">gdal_rasterize</a> utility.
1393 : * @param psOptionsForBinary (output) may be NULL (and should generally be
1394 : * NULL), otherwise (gdal_translate_bin.cpp use case) must be allocated with
1395 : * GDALFootprintOptionsForBinaryNew() prior to this function. Will be filled
1396 : * with potentially present filename, open options,...
1397 : * @return pointer to the allocated GDALFootprintOptions struct. Must be freed
1398 : * with GDALFootprintOptionsFree().
1399 : *
1400 : * @since GDAL 3.8
1401 : */
1402 :
1403 : GDALFootprintOptions *
1404 79 : GDALFootprintOptionsNew(char **papszArgv,
1405 : GDALFootprintOptionsForBinary *psOptionsForBinary)
1406 : {
1407 158 : auto psOptions = std::make_unique<GDALFootprintOptions>();
1408 :
1409 : /* -------------------------------------------------------------------- */
1410 : /* Parse arguments. */
1411 : /* -------------------------------------------------------------------- */
1412 :
1413 158 : CPLStringList aosArgv;
1414 :
1415 79 : if (papszArgv)
1416 : {
1417 78 : const int nArgc = CSLCount(papszArgv);
1418 682 : for (int i = 0; i < nArgc; i++)
1419 : {
1420 604 : aosArgv.AddString(papszArgv[i]);
1421 : }
1422 : }
1423 :
1424 : try
1425 : {
1426 : auto argParser = GDALFootprintAppOptionsGetParser(psOptions.get(),
1427 80 : psOptionsForBinary);
1428 :
1429 79 : argParser->parse_args_without_binary_name(aosArgv.List());
1430 :
1431 78 : if (argParser->is_used("-t_srs"))
1432 : {
1433 :
1434 4 : const std::string osVal(argParser->get<std::string>("-t_srs"));
1435 4 : if (psOptions->oOutputSRS.SetFromUserInput(osVal.c_str()) !=
1436 : OGRERR_NONE)
1437 : {
1438 0 : CPLError(CE_Failure, CPLE_AppDefined,
1439 : "Failed to process SRS definition: %s", osVal.c_str());
1440 0 : return nullptr;
1441 : }
1442 4 : psOptions->oOutputSRS.SetAxisMappingStrategy(
1443 : OAMS_TRADITIONAL_GIS_ORDER);
1444 : }
1445 :
1446 78 : if (argParser->is_used("-max_points"))
1447 : {
1448 : const std::string maxPoints{
1449 32 : argParser->get<std::string>("-max_points")};
1450 32 : if (maxPoints == "unlimited")
1451 : {
1452 0 : psOptions->nMaxPoints = 0;
1453 : }
1454 : else
1455 : {
1456 32 : psOptions->nMaxPoints = atoi(maxPoints.c_str());
1457 32 : if (psOptions->nMaxPoints > 0 && psOptions->nMaxPoints < 3)
1458 : {
1459 0 : CPLError(CE_Failure, CPLE_NotSupported,
1460 : "Invalid value for -max_points");
1461 0 : return nullptr;
1462 : }
1463 : }
1464 : }
1465 :
1466 78 : psOptions->bCreateOutput = !psOptions->osFormat.empty();
1467 : }
1468 1 : catch (const std::exception &err)
1469 : {
1470 1 : CPLError(CE_Failure, CPLE_AppDefined, "Unexpected exception: %s",
1471 1 : err.what());
1472 1 : return nullptr;
1473 : }
1474 :
1475 78 : if (!psOptions->bOutCSGeoref && !psOptions->oOutputSRS.IsEmpty())
1476 : {
1477 1 : CPLError(CE_Failure, CPLE_AppDefined,
1478 : "-t_cs pixel and -t_srs are mutually exclusive.");
1479 1 : return nullptr;
1480 : }
1481 :
1482 77 : if (psOptions->bClearLocation)
1483 : {
1484 2 : psOptions->osLocationFieldName.clear();
1485 : }
1486 :
1487 77 : if (psOptionsForBinary)
1488 : {
1489 7 : psOptionsForBinary->bCreateOutput = psOptions->bCreateOutput;
1490 7 : psOptionsForBinary->osFormat = psOptions->osFormat;
1491 7 : psOptionsForBinary->osDestLayerName = psOptions->osDestLayerName;
1492 : }
1493 :
1494 77 : return psOptions.release();
1495 : }
1496 :
1497 : /************************************************************************/
1498 : /* GDALFootprintOptionsFree() */
1499 : /************************************************************************/
1500 :
1501 : /**
1502 : * Frees the GDALFootprintOptions struct.
1503 : *
1504 : * @param psOptions the options struct for GDALFootprint().
1505 : *
1506 : * @since GDAL 3.8
1507 : */
1508 :
1509 150 : void GDALFootprintOptionsFree(GDALFootprintOptions *psOptions)
1510 : {
1511 150 : delete psOptions;
1512 150 : }
1513 :
1514 : /************************************************************************/
1515 : /* GDALFootprintOptionsSetProgress() */
1516 : /************************************************************************/
1517 :
1518 : /**
1519 : * Set a progress function.
1520 : *
1521 : * @param psOptions the options struct for GDALFootprint().
1522 : * @param pfnProgress the progress callback.
1523 : * @param pProgressData the user data for the progress callback.
1524 : *
1525 : * @since GDAL 3.8
1526 : */
1527 :
1528 37 : void GDALFootprintOptionsSetProgress(GDALFootprintOptions *psOptions,
1529 : GDALProgressFunc pfnProgress,
1530 : void *pProgressData)
1531 : {
1532 37 : psOptions->pfnProgress = pfnProgress ? pfnProgress : GDALDummyProgress;
1533 37 : psOptions->pProgressData = pProgressData;
1534 37 : }
|
