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