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