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