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