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