Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "raster overview refresh" subcommand
5 : * Author: Even Rouault <even dot rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2025, Even Rouault <even dot rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_raster_overview_refresh.h"
14 :
15 : #include "cpl_string.h"
16 : #include "gdal_priv.h"
17 : #include "vrtdataset.h"
18 : #include "vrt_priv.h"
19 :
20 : #include <algorithm>
21 : #include <limits>
22 :
23 : //! @cond Doxygen_Suppress
24 :
25 : #ifndef _
26 : #define _(x) (x)
27 : #endif
28 :
29 : /************************************************************************/
30 : /* GDALRasterOverviewAlgorithmRefresh() */
31 : /************************************************************************/
32 :
33 63 : GDALRasterOverviewAlgorithmRefresh::GDALRasterOverviewAlgorithmRefresh()
34 63 : : GDALAlgorithm(NAME, DESCRIPTION, HELP_URL)
35 : {
36 63 : AddProgressArg();
37 63 : AddOpenOptionsArg(&m_openOptions);
38 : AddInputDatasetArg(&m_dataset, GDAL_OF_RASTER | GDAL_OF_UPDATE,
39 : /* positionalAndRequired = */ true,
40 63 : _("Dataset (to be updated in-place, unless --external)"))
41 63 : .AddAlias("dataset");
42 126 : AddArg("external", 0, _("Refresh external overviews"), &m_readOnly)
43 126 : .AddHiddenAlias("ro")
44 63 : .AddHiddenAlias(GDAL_ARG_NAME_READ_ONLY);
45 :
46 126 : AddArg("resampling", 'r', _("Resampling method"), &m_resampling)
47 : .SetChoices("nearest", "average", "cubic", "cubicspline", "lanczos",
48 63 : "bilinear", "gauss", "average_magphase", "rms", "mode")
49 63 : .SetHiddenChoices("near", "none");
50 :
51 126 : AddArg("levels", 0, _("Levels / decimation factors"), &m_levels)
52 63 : .SetMinValueIncluded(2);
53 :
54 63 : AddBBOXArg(&m_refreshBbox, _("Bounding box to refresh"))
55 63 : .SetMutualExclusionGroup("refresh");
56 126 : AddArg("like", 0, _("Use extent of dataset(s)"), &m_like)
57 63 : .SetMutualExclusionGroup("refresh");
58 : AddArg("use-source-timestamp", 0,
59 : _("Use timestamp of VRT or GTI sources as refresh criterion"),
60 126 : &m_refreshFromSourceTimestamp)
61 63 : .SetMutualExclusionGroup("refresh");
62 63 : }
63 :
64 : /************************************************************************/
65 : /* PartialRefresh() */
66 : /************************************************************************/
67 :
68 7 : static bool PartialRefresh(GDALDataset *poDS,
69 : const std::vector<int> &anOvrIndices,
70 : const char *pszResampling, int nXOff, int nYOff,
71 : int nXSize, int nYSize, GDALProgressFunc pfnProgress,
72 : void *pProgressArg)
73 : {
74 7 : int nOvCount = 0;
75 7 : const int nBandCount = poDS->GetRasterCount();
76 14 : for (int i = 0; i < nBandCount; ++i)
77 : {
78 7 : auto poSrcBand = poDS->GetRasterBand(i + 1);
79 7 : if (i == 0)
80 7 : nOvCount = poSrcBand->GetOverviewCount();
81 0 : else if (nOvCount != poSrcBand->GetOverviewCount())
82 : {
83 0 : CPLError(CE_Failure, CPLE_AppDefined,
84 : "Not same number of overviews on all bands");
85 0 : return false;
86 : }
87 : }
88 :
89 14 : std::vector<GDALRasterBand *> apoSrcBands;
90 14 : std::vector<std::vector<GDALRasterBand *>> aapoOverviewBands;
91 14 : for (int i = 0; i < nBandCount; ++i)
92 : {
93 7 : auto poSrcBand = poDS->GetRasterBand(i + 1);
94 7 : apoSrcBands.push_back(poSrcBand);
95 14 : std::vector<GDALRasterBand *> apoOverviewBands;
96 17 : for (int nOvrIdx : anOvrIndices)
97 : {
98 10 : apoOverviewBands.push_back(poSrcBand->GetOverview(nOvrIdx));
99 : }
100 7 : aapoOverviewBands.push_back(std::move(apoOverviewBands));
101 : }
102 :
103 7 : CPLStringList aosOptions;
104 7 : aosOptions.SetNameValue("XOFF", CPLSPrintf("%d", nXOff));
105 7 : aosOptions.SetNameValue("YOFF", CPLSPrintf("%d", nYOff));
106 7 : aosOptions.SetNameValue("XSIZE", CPLSPrintf("%d", nXSize));
107 7 : aosOptions.SetNameValue("YSIZE", CPLSPrintf("%d", nYSize));
108 7 : return GDALRegenerateOverviewsMultiBand(
109 : apoSrcBands, aapoOverviewBands, pszResampling, pfnProgress,
110 14 : pProgressArg, aosOptions.List()) == CE_None;
111 : }
112 :
113 : /************************************************************************/
114 : /* PartialRefreshFromSourceTimestamp() */
115 : /************************************************************************/
116 :
117 : static bool
118 5 : PartialRefreshFromSourceTimestamp(GDALDataset *poDS, const char *pszResampling,
119 : const std::vector<int> &anOvrIndices,
120 : GDALProgressFunc pfnProgress,
121 : void *pProgressArg)
122 : {
123 : VSIStatBufL sStatOvr;
124 15 : std::string osOvr(std::string(poDS->GetDescription()) + ".ovr");
125 5 : if (VSIStatL(osOvr.c_str(), &sStatOvr) != 0)
126 : {
127 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find %s", osOvr.c_str());
128 1 : return false;
129 : }
130 4 : if (sStatOvr.st_mtime == 0)
131 : {
132 1 : CPLError(CE_Failure, CPLE_AppDefined,
133 : "Cannot get modification time of %s", osOvr.c_str());
134 1 : return false;
135 : }
136 :
137 6 : std::vector<GTISourceDesc> regions;
138 :
139 : // init slightly above zero to please Coverity Scan
140 3 : double dfTotalPixels = std::numeric_limits<double>::min();
141 :
142 3 : if (dynamic_cast<VRTDataset *>(poDS))
143 : {
144 : auto poVRTBand =
145 1 : dynamic_cast<VRTSourcedRasterBand *>(poDS->GetRasterBand(1));
146 1 : if (!poVRTBand)
147 : {
148 0 : CPLError(CE_Failure, CPLE_AppDefined,
149 : "Band is not a VRTSourcedRasterBand");
150 0 : return false;
151 : }
152 :
153 3 : for (auto &poSource : poVRTBand->m_papoSources)
154 : {
155 : auto poSimpleSource =
156 2 : dynamic_cast<VRTSimpleSource *>(poSource.get());
157 2 : if (poSimpleSource)
158 : {
159 : VSIStatBufL sStatSource;
160 2 : if (VSIStatL(poSimpleSource->GetSourceDatasetName().c_str(),
161 2 : &sStatSource) == 0)
162 : {
163 2 : if (sStatSource.st_mtime > sStatOvr.st_mtime)
164 : {
165 : double dfXOff, dfYOff, dfXSize, dfYSize;
166 1 : poSimpleSource->GetDstWindow(dfXOff, dfYOff, dfXSize,
167 : dfYSize);
168 1 : constexpr double EPS = 1e-8;
169 1 : int nXOff = static_cast<int>(dfXOff + EPS);
170 1 : int nYOff = static_cast<int>(dfYOff + EPS);
171 1 : int nXSize = static_cast<int>(dfXSize + 0.5);
172 1 : int nYSize = static_cast<int>(dfYSize + 0.5);
173 2 : if (!(nXOff > poDS->GetRasterXSize() ||
174 1 : nYOff > poDS->GetRasterYSize() || nXSize <= 0 ||
175 : nYSize <= 0))
176 : {
177 1 : if (nXOff < 0)
178 : {
179 0 : nXSize += nXOff;
180 0 : nXOff = 0;
181 : }
182 1 : if (nXOff > poDS->GetRasterXSize() - nXSize)
183 : {
184 0 : nXSize = poDS->GetRasterXSize() - nXOff;
185 : }
186 1 : if (nYOff < 0)
187 : {
188 0 : nYSize += nYOff;
189 0 : nYOff = 0;
190 : }
191 1 : if (nYOff > poDS->GetRasterYSize() - nYSize)
192 : {
193 0 : nYSize = poDS->GetRasterYSize() - nYOff;
194 : }
195 :
196 1 : dfTotalPixels +=
197 1 : static_cast<double>(nXSize) * nYSize;
198 2 : GTISourceDesc region;
199 : region.osFilename =
200 1 : poSimpleSource->GetSourceDatasetName();
201 1 : region.nDstXOff = nXOff;
202 1 : region.nDstYOff = nYOff;
203 1 : region.nDstXSize = nXSize;
204 1 : region.nDstYSize = nYSize;
205 1 : regions.push_back(std::move(region));
206 : }
207 : }
208 : }
209 : }
210 : }
211 : }
212 : #ifdef GTI_DRIVER_DISABLED_OR_PLUGIN
213 : else if (poDS->GetDriver() &&
214 : EQUAL(poDS->GetDriver()->GetDescription(), "GTI"))
215 : {
216 : CPLError(CE_Failure, CPLE_NotSupported,
217 : "--use-source-timestamp only works on a GTI "
218 : "dataset if the GTI driver is not built as a plugin, "
219 : "but in core library");
220 : return false;
221 : }
222 : #else
223 2 : else if (auto poGTIDS = GDALDatasetCastToGTIDataset(poDS))
224 : {
225 1 : regions = GTIGetSourcesMoreRecentThan(poGTIDS, sStatOvr.st_mtime);
226 2 : for (const auto ®ion : regions)
227 : {
228 1 : dfTotalPixels +=
229 1 : static_cast<double>(region.nDstXSize) * region.nDstYSize;
230 : }
231 : }
232 : #endif
233 : else
234 : {
235 1 : CPLError(CE_Failure, CPLE_AppDefined,
236 : "--use-source-timestamp only works on a VRT or GTI "
237 : "dataset");
238 1 : return false;
239 : }
240 :
241 2 : bool bRet = true;
242 2 : if (!regions.empty())
243 : {
244 2 : double dfCurPixels = 0;
245 4 : for (const auto ®ion : regions)
246 : {
247 2 : if (bRet)
248 : {
249 2 : CPLDebug("GDAL", "Refresh from source %s",
250 : region.osFilename.c_str());
251 2 : double dfNextCurPixels =
252 : dfCurPixels +
253 2 : static_cast<double>(region.nDstXSize) * region.nDstYSize;
254 2 : void *pScaledProgress = GDALCreateScaledProgress(
255 : dfCurPixels / dfTotalPixels,
256 : dfNextCurPixels / dfTotalPixels, pfnProgress, pProgressArg);
257 2 : bRet = PartialRefresh(
258 2 : poDS, anOvrIndices, pszResampling, region.nDstXOff,
259 2 : region.nDstYOff, region.nDstXSize, region.nDstYSize,
260 : pScaledProgress ? GDALScaledProgress : nullptr,
261 : pScaledProgress);
262 2 : GDALDestroyScaledProgress(pScaledProgress);
263 2 : dfCurPixels = dfNextCurPixels;
264 : }
265 : }
266 : }
267 : else
268 : {
269 0 : CPLDebug("GDAL", "No source is more recent than the overviews");
270 : }
271 :
272 2 : return bRet;
273 : }
274 :
275 : /************************************************************************/
276 : /* PartialRefreshFromSourceExtent() */
277 : /************************************************************************/
278 :
279 6 : static bool PartialRefreshFromSourceExtent(
280 : GDALDataset *poDS, const std::vector<std::string> &sources,
281 : const char *pszResampling, const std::vector<int> &anOvrIndices,
282 : GDALProgressFunc pfnProgress, void *pProgressArg)
283 : {
284 6 : GDALGeoTransform gt;
285 6 : if (poDS->GetGeoTransform(gt) != CE_None)
286 : {
287 1 : CPLError(CE_Failure, CPLE_AppDefined, "Dataset has no geotransform");
288 1 : return false;
289 : }
290 5 : GDALGeoTransform invGT;
291 5 : if (!gt.GetInverse(invGT))
292 : {
293 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
294 1 : return false;
295 : }
296 :
297 : struct Region
298 : {
299 : std::string osFileName{};
300 : int nXOff = 0;
301 : int nYOff = 0;
302 : int nXSize = 0;
303 : int nYSize = 0;
304 : };
305 :
306 8 : std::vector<Region> regions;
307 :
308 : // init slightly above zero to please Coverity Scan
309 4 : double dfTotalPixels = std::numeric_limits<double>::min();
310 6 : for (const std::string &filename : sources)
311 : {
312 : auto poSrcDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
313 4 : filename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR));
314 4 : if (!poSrcDS)
315 1 : return false;
316 :
317 3 : GDALGeoTransform srcGT;
318 3 : if (poSrcDS->GetGeoTransform(srcGT) != CE_None)
319 : {
320 1 : CPLError(CE_Failure, CPLE_AppDefined,
321 : "Source dataset has no geotransform");
322 1 : return false;
323 : }
324 :
325 2 : const double dfULX = srcGT[0];
326 2 : const double dfULY = srcGT[3];
327 2 : const double dfLRX = srcGT[0] + poSrcDS->GetRasterXSize() * srcGT[1] +
328 2 : poSrcDS->GetRasterYSize() * srcGT[2];
329 2 : const double dfLRY = srcGT[3] + poSrcDS->GetRasterXSize() * srcGT[4] +
330 2 : poSrcDS->GetRasterYSize() * srcGT[5];
331 2 : const double dfX1 = invGT[0] + invGT[1] * dfULX + invGT[2] * dfULY;
332 2 : const double dfY1 = invGT[3] + invGT[4] * dfULX + invGT[5] * dfULY;
333 2 : const double dfX2 = invGT[0] + invGT[1] * dfLRX + invGT[2] * dfLRY;
334 2 : const double dfY2 = invGT[3] + invGT[4] * dfLRX + invGT[5] * dfLRY;
335 2 : constexpr double EPS = 1e-8;
336 : const int nXOff =
337 2 : static_cast<int>(std::max(0.0, std::min(dfX1, dfX2)) + EPS);
338 : const int nYOff =
339 2 : static_cast<int>(std::max(0.0, std::min(dfY1, dfY2)) + EPS);
340 : const int nXSize =
341 2 : static_cast<int>(
342 2 : std::ceil(std::min(static_cast<double>(poDS->GetRasterXSize()),
343 4 : std::max(dfX1, dfX2)) -
344 : EPS)) -
345 2 : nXOff;
346 : const int nYSize =
347 2 : static_cast<int>(
348 2 : std::ceil(std::min(static_cast<double>(poDS->GetRasterYSize()),
349 4 : std::max(dfY1, dfY2)) -
350 : EPS)) -
351 2 : nYOff;
352 :
353 2 : dfTotalPixels += static_cast<double>(nXSize) * nYSize;
354 4 : Region region;
355 2 : region.osFileName = filename;
356 2 : region.nXOff = nXOff;
357 2 : region.nYOff = nYOff;
358 2 : region.nXSize = nXSize;
359 2 : region.nYSize = nYSize;
360 2 : regions.push_back(std::move(region));
361 : }
362 :
363 2 : bool bRet = true;
364 2 : double dfCurPixels = 0;
365 4 : for (const auto ®ion : regions)
366 : {
367 2 : if (bRet)
368 : {
369 2 : CPLDebug("GDAL", "Refresh from source %s",
370 : region.osFileName.c_str());
371 2 : double dfNextCurPixels =
372 : dfCurPixels +
373 2 : static_cast<double>(region.nXSize) * region.nYSize;
374 : // coverity[divide_by_zero]
375 2 : void *pScaledProgress = GDALCreateScaledProgress(
376 : dfCurPixels / dfTotalPixels, dfNextCurPixels / dfTotalPixels,
377 : pfnProgress, pProgressArg);
378 : bRet =
379 2 : PartialRefresh(poDS, anOvrIndices, pszResampling, region.nXOff,
380 2 : region.nYOff, region.nXSize, region.nYSize,
381 : pScaledProgress ? GDALScaledProgress : nullptr,
382 : pScaledProgress);
383 2 : GDALDestroyScaledProgress(pScaledProgress);
384 2 : dfCurPixels = dfNextCurPixels;
385 : }
386 : }
387 :
388 2 : return bRet;
389 : }
390 :
391 : /************************************************************************/
392 : /* PartialRefreshFromBBOX() */
393 : /************************************************************************/
394 :
395 5 : static bool PartialRefreshFromBBOX(GDALDataset *poDS,
396 : const std::vector<double> &bbox,
397 : const char *pszResampling,
398 : const std::vector<int> &anOvrIndices,
399 : GDALProgressFunc pfnProgress,
400 : void *pProgressArg)
401 : {
402 5 : const double dfULX = bbox[0];
403 5 : const double dfLRY = bbox[1];
404 5 : const double dfLRX = bbox[2];
405 5 : const double dfULY = bbox[3];
406 :
407 5 : GDALGeoTransform gt;
408 5 : if (poDS->GetGeoTransform(gt) != CE_None)
409 : {
410 1 : CPLError(CE_Failure, CPLE_AppDefined, "Dataset has no geotransform");
411 1 : return false;
412 : }
413 4 : GDALGeoTransform invGT;
414 4 : if (!gt.GetInverse(invGT))
415 : {
416 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
417 1 : return false;
418 : }
419 3 : const double dfX1 = invGT[0] + invGT[1] * dfULX + invGT[2] * dfULY;
420 3 : const double dfY1 = invGT[3] + invGT[4] * dfULX + invGT[5] * dfULY;
421 3 : const double dfX2 = invGT[0] + invGT[1] * dfLRX + invGT[2] * dfLRY;
422 3 : const double dfY2 = invGT[3] + invGT[4] * dfLRX + invGT[5] * dfLRY;
423 3 : constexpr double EPS = 1e-8;
424 : const int nXOff =
425 3 : static_cast<int>(std::max(0.0, std::min(dfX1, dfX2)) + EPS);
426 : const int nYOff =
427 3 : static_cast<int>(std::max(0.0, std::min(dfY1, dfY2)) + EPS);
428 3 : const int nXSize = static_cast<int>(std::ceil(
429 3 : std::min(static_cast<double>(poDS->GetRasterXSize()),
430 6 : std::max(dfX1, dfX2)) -
431 : EPS)) -
432 3 : nXOff;
433 3 : const int nYSize = static_cast<int>(std::ceil(
434 3 : std::min(static_cast<double>(poDS->GetRasterYSize()),
435 6 : std::max(dfY1, dfY2)) -
436 : EPS)) -
437 3 : nYOff;
438 3 : return PartialRefresh(poDS, anOvrIndices, pszResampling, nXOff, nYOff,
439 3 : nXSize, nYSize, pfnProgress, pProgressArg);
440 : }
441 :
442 : /************************************************************************/
443 : /* GDALRasterOverviewAlgorithmRefresh::RunImpl() */
444 : /************************************************************************/
445 :
446 20 : bool GDALRasterOverviewAlgorithmRefresh::RunImpl(GDALProgressFunc pfnProgress,
447 : void *pProgressData)
448 : {
449 20 : auto poDS = m_dataset.GetDatasetRef();
450 20 : CPLAssert(poDS);
451 20 : if (poDS->GetRasterCount() == 0)
452 : {
453 1 : ReportError(CE_Failure, CPLE_AppDefined, "Dataset has no raster band");
454 1 : return false;
455 : }
456 :
457 19 : auto poBand = poDS->GetRasterBand(1);
458 19 : const int nOvCount = poBand->GetOverviewCount();
459 :
460 38 : std::vector<int> levels = m_levels;
461 :
462 : // If no levels are specified, reuse the potentially existing ones.
463 19 : if (levels.empty())
464 : {
465 38 : for (int iOvr = 0; iOvr < nOvCount; ++iOvr)
466 : {
467 20 : auto poOverview = poBand->GetOverview(iOvr);
468 20 : if (poOverview)
469 : {
470 20 : const int nOvFactor = GDALComputeOvFactor(
471 : poOverview->GetXSize(), poBand->GetXSize(),
472 20 : poOverview->GetYSize(), poBand->GetYSize());
473 20 : levels.push_back(nOvFactor);
474 : }
475 : }
476 : }
477 19 : if (levels.empty())
478 : {
479 1 : ReportError(CE_Failure, CPLE_AppDefined, "No overviews to refresh");
480 1 : return false;
481 : }
482 :
483 36 : std::vector<int> anOvrIndices;
484 38 : for (int nLevel : levels)
485 : {
486 21 : int nIdx = -1;
487 25 : for (int iOvr = 0; iOvr < nOvCount; iOvr++)
488 : {
489 24 : auto poOverview = poBand->GetOverview(iOvr);
490 24 : if (poOverview)
491 : {
492 24 : const int nOvFactor = GDALComputeOvFactor(
493 : poOverview->GetXSize(), poBand->GetXSize(),
494 : poOverview->GetYSize(), poBand->GetYSize());
495 28 : if (nOvFactor == nLevel ||
496 4 : nOvFactor == GDALOvLevelAdjust2(nLevel, poBand->GetXSize(),
497 : poBand->GetYSize()))
498 : {
499 20 : nIdx = iOvr;
500 20 : break;
501 : }
502 : }
503 : }
504 21 : if (nIdx < 0)
505 : {
506 1 : CPLError(CE_Failure, CPLE_AppDefined,
507 : "Cannot find overview level with subsampling factor of %d",
508 : nLevel);
509 1 : return false;
510 : }
511 20 : CPLDebug("GDAL", "Refreshing overview idx %d", nIdx);
512 20 : anOvrIndices.push_back(nIdx);
513 : }
514 :
515 34 : std::string resampling = m_resampling;
516 17 : if (resampling.empty())
517 : {
518 : const char *pszResampling =
519 10 : poBand->GetOverview(0)->GetMetadataItem("RESAMPLING");
520 10 : if (pszResampling)
521 : {
522 2 : resampling = pszResampling;
523 2 : CPLDebug("GDAL",
524 : "Reusing resampling method %s from existing "
525 : "overview",
526 : pszResampling);
527 : }
528 : }
529 17 : if (resampling.empty())
530 8 : resampling = "nearest";
531 :
532 17 : if (m_refreshFromSourceTimestamp)
533 : {
534 5 : return PartialRefreshFromSourceTimestamp(
535 5 : poDS, resampling.c_str(), anOvrIndices, pfnProgress, pProgressData);
536 : }
537 12 : else if (!m_refreshBbox.empty())
538 : {
539 5 : return PartialRefreshFromBBOX(poDS, m_refreshBbox, resampling.c_str(),
540 5 : anOvrIndices, pfnProgress, pProgressData);
541 : }
542 7 : else if (!m_like.empty())
543 : {
544 6 : return PartialRefreshFromSourceExtent(poDS, m_like, resampling.c_str(),
545 : anOvrIndices, pfnProgress,
546 6 : pProgressData);
547 : }
548 : else
549 : {
550 1 : return GDALBuildOverviews(
551 : GDALDataset::ToHandle(poDS), resampling.c_str(),
552 1 : static_cast<int>(levels.size()), levels.data(), 0, nullptr,
553 1 : pfnProgress, pProgressData) == CE_None;
554 : }
555 : }
556 :
557 : //! @endcond
|
