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 : double adfGeoTransform[6];
286 6 : if (poDS->GetGeoTransform(adfGeoTransform) != CE_None)
287 : {
288 1 : CPLError(CE_Failure, CPLE_AppDefined, "Dataset has no geotransform");
289 1 : return false;
290 : }
291 : double adfInvGT[6];
292 5 : if (!GDALInvGeoTransform(adfGeoTransform, adfInvGT))
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 : double adfSrcGT[6];
319 3 : if (poSrcDS->GetGeoTransform(adfSrcGT) != 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 = adfSrcGT[0];
327 2 : const double dfULY = adfSrcGT[3];
328 4 : const double dfLRX = adfSrcGT[0] +
329 2 : poSrcDS->GetRasterXSize() * adfSrcGT[1] +
330 2 : poSrcDS->GetRasterYSize() * adfSrcGT[2];
331 4 : const double dfLRY = adfSrcGT[3] +
332 2 : poSrcDS->GetRasterXSize() * adfSrcGT[4] +
333 2 : poSrcDS->GetRasterYSize() * adfSrcGT[5];
334 2 : const double dfX1 =
335 2 : adfInvGT[0] + adfInvGT[1] * dfULX + adfInvGT[2] * dfULY;
336 2 : const double dfY1 =
337 2 : adfInvGT[3] + adfInvGT[4] * dfULX + adfInvGT[5] * dfULY;
338 2 : const double dfX2 =
339 2 : adfInvGT[0] + adfInvGT[1] * dfLRX + adfInvGT[2] * dfLRY;
340 2 : const double dfY2 =
341 2 : adfInvGT[3] + adfInvGT[4] * dfLRX + adfInvGT[5] * dfLRY;
342 2 : constexpr double EPS = 1e-8;
343 : const int nXOff =
344 2 : static_cast<int>(std::max(0.0, std::min(dfX1, dfX2)) + EPS);
345 : const int nYOff =
346 2 : static_cast<int>(std::max(0.0, std::min(dfY1, dfY2)) + EPS);
347 : const int nXSize =
348 2 : static_cast<int>(
349 2 : std::ceil(std::min(static_cast<double>(poDS->GetRasterXSize()),
350 4 : std::max(dfX1, dfX2)) -
351 : EPS)) -
352 2 : nXOff;
353 : const int nYSize =
354 2 : static_cast<int>(
355 2 : std::ceil(std::min(static_cast<double>(poDS->GetRasterYSize()),
356 4 : std::max(dfY1, dfY2)) -
357 : EPS)) -
358 2 : nYOff;
359 :
360 2 : dfTotalPixels += static_cast<double>(nXSize) * nYSize;
361 4 : Region region;
362 2 : region.osFileName = filename;
363 2 : region.nXOff = nXOff;
364 2 : region.nYOff = nYOff;
365 2 : region.nXSize = nXSize;
366 2 : region.nYSize = nYSize;
367 2 : regions.push_back(std::move(region));
368 : }
369 :
370 2 : bool bRet = true;
371 2 : double dfCurPixels = 0;
372 4 : for (const auto ®ion : regions)
373 : {
374 2 : if (bRet)
375 : {
376 2 : CPLDebug("GDAL", "Refresh from source %s",
377 : region.osFileName.c_str());
378 2 : double dfNextCurPixels =
379 : dfCurPixels +
380 2 : static_cast<double>(region.nXSize) * region.nYSize;
381 : // coverity[divide_by_zero]
382 2 : void *pScaledProgress = GDALCreateScaledProgress(
383 : dfCurPixels / dfTotalPixels, dfNextCurPixels / dfTotalPixels,
384 : pfnProgress, pProgressArg);
385 : bRet =
386 2 : PartialRefresh(poDS, anOvrIndices, pszResampling, region.nXOff,
387 2 : region.nYOff, region.nXSize, region.nYSize,
388 : pScaledProgress ? GDALScaledProgress : nullptr,
389 : pScaledProgress);
390 2 : GDALDestroyScaledProgress(pScaledProgress);
391 2 : dfCurPixels = dfNextCurPixels;
392 : }
393 : }
394 :
395 2 : return bRet;
396 : }
397 :
398 : /************************************************************************/
399 : /* PartialRefreshFromBBOX() */
400 : /************************************************************************/
401 :
402 5 : static bool PartialRefreshFromBBOX(GDALDataset *poDS,
403 : const std::vector<double> &bbox,
404 : const char *pszResampling,
405 : const std::vector<int> &anOvrIndices,
406 : GDALProgressFunc pfnProgress,
407 : void *pProgressArg)
408 : {
409 5 : const double dfULX = bbox[0];
410 5 : const double dfLRY = bbox[1];
411 5 : const double dfLRX = bbox[2];
412 5 : const double dfULY = bbox[3];
413 :
414 : double adfGeoTransform[6];
415 5 : if (poDS->GetGeoTransform(adfGeoTransform) != CE_None)
416 : {
417 1 : CPLError(CE_Failure, CPLE_AppDefined, "Dataset has no geotransform");
418 1 : return false;
419 : }
420 : double adfInvGT[6];
421 4 : if (!GDALInvGeoTransform(adfGeoTransform, adfInvGT))
422 : {
423 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
424 1 : return false;
425 : }
426 3 : const double dfX1 = adfInvGT[0] + adfInvGT[1] * dfULX + adfInvGT[2] * dfULY;
427 3 : const double dfY1 = adfInvGT[3] + adfInvGT[4] * dfULX + adfInvGT[5] * dfULY;
428 3 : const double dfX2 = adfInvGT[0] + adfInvGT[1] * dfLRX + adfInvGT[2] * dfLRY;
429 3 : const double dfY2 = adfInvGT[3] + adfInvGT[4] * dfLRX + adfInvGT[5] * dfLRY;
430 3 : constexpr double EPS = 1e-8;
431 : const int nXOff =
432 3 : static_cast<int>(std::max(0.0, std::min(dfX1, dfX2)) + EPS);
433 : const int nYOff =
434 3 : static_cast<int>(std::max(0.0, std::min(dfY1, dfY2)) + EPS);
435 3 : const int nXSize = static_cast<int>(std::ceil(
436 3 : std::min(static_cast<double>(poDS->GetRasterXSize()),
437 6 : std::max(dfX1, dfX2)) -
438 : EPS)) -
439 3 : nXOff;
440 3 : const int nYSize = static_cast<int>(std::ceil(
441 3 : std::min(static_cast<double>(poDS->GetRasterYSize()),
442 6 : std::max(dfY1, dfY2)) -
443 : EPS)) -
444 3 : nYOff;
445 3 : return PartialRefresh(poDS, anOvrIndices, pszResampling, nXOff, nYOff,
446 3 : nXSize, nYSize, pfnProgress, pProgressArg);
447 : }
448 :
449 : /************************************************************************/
450 : /* GDALRasterOverviewAlgorithmRefresh::RunImpl() */
451 : /************************************************************************/
452 :
453 20 : bool GDALRasterOverviewAlgorithmRefresh::RunImpl(GDALProgressFunc pfnProgress,
454 : void *pProgressData)
455 : {
456 20 : auto poDS = m_dataset.GetDatasetRef();
457 20 : CPLAssert(poDS);
458 20 : if (poDS->GetRasterCount() == 0)
459 : {
460 1 : ReportError(CE_Failure, CPLE_AppDefined, "Dataset has no raster band");
461 1 : return false;
462 : }
463 :
464 19 : auto poBand = poDS->GetRasterBand(1);
465 19 : const int nOvCount = poBand->GetOverviewCount();
466 :
467 38 : std::vector<int> levels = m_levels;
468 :
469 : // If no levels are specified, reuse the potentially existing ones.
470 19 : if (levels.empty())
471 : {
472 38 : for (int iOvr = 0; iOvr < nOvCount; ++iOvr)
473 : {
474 20 : auto poOverview = poBand->GetOverview(iOvr);
475 20 : if (poOverview)
476 : {
477 20 : const int nOvFactor = GDALComputeOvFactor(
478 : poOverview->GetXSize(), poBand->GetXSize(),
479 20 : poOverview->GetYSize(), poBand->GetYSize());
480 20 : levels.push_back(nOvFactor);
481 : }
482 : }
483 : }
484 19 : if (levels.empty())
485 : {
486 1 : ReportError(CE_Failure, CPLE_AppDefined, "No overviews to refresh");
487 1 : return false;
488 : }
489 :
490 36 : std::vector<int> anOvrIndices;
491 38 : for (int nLevel : levels)
492 : {
493 21 : int nIdx = -1;
494 25 : for (int iOvr = 0; iOvr < nOvCount; iOvr++)
495 : {
496 24 : auto poOverview = poBand->GetOverview(iOvr);
497 24 : if (poOverview)
498 : {
499 24 : const int nOvFactor = GDALComputeOvFactor(
500 : poOverview->GetXSize(), poBand->GetXSize(),
501 : poOverview->GetYSize(), poBand->GetYSize());
502 28 : if (nOvFactor == nLevel ||
503 4 : nOvFactor == GDALOvLevelAdjust2(nLevel, poBand->GetXSize(),
504 : poBand->GetYSize()))
505 : {
506 20 : nIdx = iOvr;
507 20 : break;
508 : }
509 : }
510 : }
511 21 : if (nIdx < 0)
512 : {
513 1 : CPLError(CE_Failure, CPLE_AppDefined,
514 : "Cannot find overview level with subsampling factor of %d",
515 : nLevel);
516 1 : return false;
517 : }
518 20 : CPLDebug("GDAL", "Refreshing overview idx %d", nIdx);
519 20 : anOvrIndices.push_back(nIdx);
520 : }
521 :
522 34 : std::string resampling = m_resampling;
523 17 : if (resampling.empty())
524 : {
525 : const char *pszResampling =
526 10 : poBand->GetOverview(0)->GetMetadataItem("RESAMPLING");
527 10 : if (pszResampling)
528 : {
529 2 : resampling = pszResampling;
530 2 : CPLDebug("GDAL",
531 : "Reusing resampling method %s from existing "
532 : "overview",
533 : pszResampling);
534 : }
535 : }
536 17 : if (resampling.empty())
537 8 : resampling = "nearest";
538 :
539 17 : if (m_refreshFromSourceTimestamp)
540 : {
541 5 : return PartialRefreshFromSourceTimestamp(
542 5 : poDS, resampling.c_str(), anOvrIndices, pfnProgress, pProgressData);
543 : }
544 12 : else if (!m_refreshBbox.empty())
545 : {
546 5 : return PartialRefreshFromBBOX(poDS, m_refreshBbox, resampling.c_str(),
547 5 : anOvrIndices, pfnProgress, pProgressData);
548 : }
549 7 : else if (!m_like.empty())
550 : {
551 6 : return PartialRefreshFromSourceExtent(poDS, m_like, resampling.c_str(),
552 : anOvrIndices, pfnProgress,
553 6 : pProgressData);
554 : }
555 : else
556 : {
557 1 : return GDALBuildOverviews(
558 : GDALDataset::ToHandle(poDS), resampling.c_str(),
559 1 : static_cast<int>(levels.size()), levels.data(), 0, nullptr,
560 1 : pfnProgress, pProgressData) == CE_None;
561 : }
562 : }
563 :
564 : //! @endcond
|
