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