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