Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "raster update" 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_update.h"
14 :
15 : #include "cpl_conv.h"
16 :
17 : #include "gdal_priv.h"
18 : #include "gdal_utils.h"
19 : #include "gdalalg_raster_reproject.h" // for GDALRasterReprojectUtils
20 : #include "gdalalg_raster_overview_refresh.h"
21 : #include "ogr_spatialref.h"
22 :
23 : #include <algorithm>
24 : #include <array>
25 : #include <cmath>
26 : #include <tuple>
27 :
28 : //! @cond Doxygen_Suppress
29 :
30 : #ifndef _
31 : #define _(x) (x)
32 : #endif
33 :
34 : /************************************************************************/
35 : /* GDALRasterUpdateAlgorithm::GDALRasterUpdateAlgorithm() */
36 : /************************************************************************/
37 :
38 92 : GDALRasterUpdateAlgorithm::GDALRasterUpdateAlgorithm(bool standaloneStep)
39 : : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL,
40 0 : ConstructorOptions()
41 92 : .SetStandaloneStep(standaloneStep)
42 92 : .SetInputDatasetMaxCount(1)
43 184 : .SetAddDefaultArguments(false))
44 : {
45 92 : if (standaloneStep)
46 : {
47 55 : AddProgressArg();
48 55 : AddRasterInputArgs(/* openForMixedRasterVector = */ false,
49 : /* hiddenForCLI = */ false);
50 : }
51 : else
52 : {
53 37 : AddRasterHiddenInputDatasetArg();
54 : }
55 :
56 92 : AddOutputDatasetArg(&m_outputDataset, GDAL_OF_RASTER)
57 92 : .SetDatasetInputFlags(GADV_NAME | GADV_OBJECT);
58 :
59 92 : m_update = true;
60 92 : AddUpdateArg(&m_update).SetDefault(true).SetHidden();
61 :
62 184 : AddArg("geometry", 0, _("Clipping geometry (WKT or GeoJSON)"), &m_geometry)
63 92 : .SetMutualExclusionGroup("bbox-geometry-like");
64 184 : AddArg("geometry-crs", 0, _("CRS of clipping geometry"), &m_geometryCrs)
65 184 : .SetIsCRSArg()
66 92 : .AddHiddenAlias("geometry_srs");
67 :
68 92 : GDALRasterReprojectUtils::AddResamplingArg(this, m_resampling);
69 :
70 92 : GDALRasterReprojectUtils::AddWarpOptTransformOptErrorThresholdArg(
71 92 : this, m_warpOptions, m_transformOptions, m_errorThreshold);
72 :
73 : AddArg("no-update-overviews", 0, _("Do not update existing overviews"),
74 92 : &m_noUpdateOverviews);
75 92 : }
76 :
77 : /************************************************************************/
78 : /* GDALRasterUpdateAlgorithm::RunStep() */
79 : /************************************************************************/
80 :
81 11 : bool GDALRasterUpdateAlgorithm::RunStep(GDALPipelineStepRunContext &stepCtxt)
82 : {
83 11 : auto poSrcDS = m_inputDataset[0].GetDatasetRef();
84 11 : CPLAssert(poSrcDS);
85 :
86 11 : auto poDstDS = m_outputDataset.GetDatasetRef();
87 11 : CPLAssert(poDstDS);
88 11 : CPLAssert(poDstDS->GetAccess() == GA_Update);
89 :
90 11 : std::unique_ptr<OGRGeometry> poClipGeom;
91 22 : std::string errMsg;
92 11 : if (!m_geometry.empty())
93 : {
94 2 : std::tie(poClipGeom, errMsg) = GetClipGeometry();
95 2 : if (!poClipGeom)
96 : {
97 1 : ReportError(CE_Failure, CPLE_AppDefined, "%s", errMsg.c_str());
98 1 : return false;
99 : }
100 : }
101 :
102 10 : auto poSrcDriver = poSrcDS->GetDriver();
103 10 : auto poDstDriver = poDstDS->GetDriver();
104 19 : if (poSrcDS == poDstDS ||
105 9 : (poSrcDriver && poDstDriver &&
106 9 : !EQUAL(poSrcDriver->GetDescription(), "MEM") &&
107 3 : !EQUAL(poDstDriver->GetDescription(), "MEM") &&
108 3 : strcmp(poSrcDS->GetDescription(), poDstDS->GetDescription()) == 0))
109 : {
110 2 : ReportError(CE_Failure, CPLE_NotSupported,
111 : "Source and destination datasets must be different");
112 2 : return false;
113 : }
114 :
115 16 : CPLStringList aosOptions;
116 8 : if (!m_resampling.empty())
117 : {
118 8 : aosOptions.AddString("-r");
119 8 : aosOptions.AddString(m_resampling.c_str());
120 : }
121 9 : for (const std::string &opt : m_warpOptions)
122 : {
123 1 : aosOptions.AddString("-wo");
124 1 : aosOptions.AddString(opt.c_str());
125 : }
126 9 : for (const std::string &opt : m_transformOptions)
127 : {
128 1 : aosOptions.AddString("-to");
129 1 : aosOptions.AddString(opt.c_str());
130 : }
131 8 : if (std::isfinite(m_errorThreshold))
132 : {
133 1 : aosOptions.AddString("-et");
134 1 : aosOptions.AddString(CPLSPrintf("%.17g", m_errorThreshold));
135 : }
136 :
137 8 : if (poClipGeom)
138 : {
139 1 : aosOptions.AddString("-cutline");
140 1 : aosOptions.AddString(poClipGeom->exportToWkt());
141 : }
142 :
143 8 : bool bOvrCanBeUpdated = false;
144 8 : std::vector<double> overviewRefreshBBox;
145 10 : if (poDstDS->GetRasterBand(1)->GetOverviewCount() > 0 &&
146 2 : !m_noUpdateOverviews)
147 : {
148 2 : GDALGeoTransform gt;
149 2 : const auto poSrcCRS = poSrcDS->GetSpatialRef();
150 2 : const auto poDstCRS = poDstDS->GetSpatialRef();
151 2 : const bool bBothCRS = poSrcCRS && poDstCRS;
152 2 : const bool bBothNoCRS = !poSrcCRS && !poDstCRS;
153 2 : if ((bBothCRS || bBothNoCRS) && poSrcDS->GetGeoTransform(gt) == CE_None)
154 : {
155 : auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
156 1 : bBothCRS ? OGRCreateCoordinateTransformation(poSrcCRS, poDstCRS)
157 5 : : nullptr);
158 2 : if (bBothNoCRS || poCT)
159 : {
160 2 : const double dfTLX = gt.xorig;
161 2 : const double dfTLY = gt.yorig;
162 :
163 2 : double dfTRX = 0;
164 2 : double dfTRY = 0;
165 2 : gt.Apply(poSrcDS->GetRasterXSize(), 0, &dfTRX, &dfTRY);
166 :
167 2 : double dfBLX = 0;
168 2 : double dfBLY = 0;
169 2 : gt.Apply(0, poSrcDS->GetRasterYSize(), &dfBLX, &dfBLY);
170 :
171 2 : double dfBRX = 0;
172 2 : double dfBRY = 0;
173 2 : gt.Apply(poSrcDS->GetRasterXSize(), poSrcDS->GetRasterYSize(),
174 : &dfBRX, &dfBRY);
175 :
176 : const double dfXMin =
177 2 : std::min(std::min(dfTLX, dfTRX), std::min(dfBLX, dfBRX));
178 : const double dfYMin =
179 2 : std::min(std::min(dfTLY, dfTRY), std::min(dfBLY, dfBRY));
180 : const double dfXMax =
181 2 : std::max(std::max(dfTLX, dfTRX), std::max(dfBLX, dfBRX));
182 : const double dfYMax =
183 2 : std::max(std::max(dfTLY, dfTRY), std::max(dfBLY, dfBRY));
184 2 : double dfOutXMin = dfXMin;
185 2 : double dfOutYMin = dfYMin;
186 2 : double dfOutXMax = dfXMax;
187 2 : double dfOutYMax = dfYMax;
188 3 : if (!poCT || poCT->TransformBounds(
189 : dfXMin, dfYMin, dfXMax, dfYMax, &dfOutXMin,
190 1 : &dfOutYMin, &dfOutXMax, &dfOutYMax, 21))
191 : {
192 2 : bOvrCanBeUpdated = true;
193 2 : CPLDebug("update",
194 : "Refresh overviews from (%f,%f) to (%f,%f)",
195 : dfOutXMin, dfOutYMin, dfOutXMax, dfOutYMax);
196 4 : overviewRefreshBBox = std::vector<double>{
197 2 : dfOutXMin, dfOutYMin, dfOutXMax, dfOutYMax};
198 : }
199 : }
200 : }
201 2 : if (!bOvrCanBeUpdated)
202 : {
203 0 : ReportError(CE_Warning, CPLE_AppDefined,
204 : "Overviews can not be updated");
205 : }
206 : }
207 :
208 8 : bool bOK = false;
209 : GDALWarpAppOptions *psOptions =
210 8 : GDALWarpAppOptionsNew(aosOptions.List(), nullptr);
211 8 : if (psOptions)
212 : {
213 : std::unique_ptr<void, decltype(&GDALDestroyScaledProgress)> pScaledData(
214 16 : nullptr, GDALDestroyScaledProgress);
215 8 : auto pfnProgress = stepCtxt.m_pfnProgress;
216 8 : void *pProgressData = stepCtxt.m_pProgressData;
217 8 : if (pfnProgress)
218 : {
219 3 : pScaledData.reset(
220 : GDALCreateScaledProgress(0.0, bOvrCanBeUpdated ? 0.75 : 1.0,
221 : pfnProgress, pProgressData));
222 3 : GDALWarpAppOptionsSetProgress(psOptions, GDALScaledProgress,
223 : pScaledData.get());
224 : }
225 :
226 8 : GDALDatasetH hSrcDS = GDALDataset::ToHandle(poSrcDS);
227 8 : GDALDatasetH hDstDS = GDALDataset::ToHandle(poDstDS);
228 8 : auto poRetDS = GDALDataset::FromHandle(
229 : GDALWarp(nullptr, hDstDS, 1, &hSrcDS, psOptions, nullptr));
230 8 : GDALWarpAppOptionsFree(psOptions);
231 :
232 8 : bOK = poRetDS != nullptr;
233 8 : if (bOK && bOvrCanBeUpdated)
234 : {
235 2 : GDALRasterOverviewAlgorithmRefresh refresh;
236 2 : refresh.GetArg("dataset")->Set(poRetDS);
237 2 : if (!m_resampling.empty())
238 2 : refresh.GetArg("resampling")->Set(m_resampling);
239 2 : refresh.GetArg("bbox")->Set(overviewRefreshBBox);
240 2 : pScaledData.reset(GDALCreateScaledProgress(0.75, 1.0, pfnProgress,
241 : pProgressData));
242 2 : bOK = refresh.Run(pScaledData ? GDALScaledProgress : nullptr,
243 : pScaledData.get());
244 : }
245 8 : if (bOK && pfnProgress)
246 3 : pfnProgress(1.0, "", pProgressData);
247 : }
248 :
249 8 : return bOK;
250 : }
251 :
252 : /************************************************************************/
253 : /* GDALRasterUpdateAlgorithm::RunImpl() */
254 : /************************************************************************/
255 :
256 9 : bool GDALRasterUpdateAlgorithm::RunImpl(GDALProgressFunc pfnProgress,
257 : void *pProgressData)
258 : {
259 9 : GDALPipelineStepRunContext stepCtxt;
260 9 : stepCtxt.m_pfnProgress = pfnProgress;
261 9 : stepCtxt.m_pProgressData = pProgressData;
262 18 : return RunStep(stepCtxt);
263 : }
264 :
265 : GDALRasterUpdateAlgorithmStandalone::~GDALRasterUpdateAlgorithmStandalone() =
266 : default;
267 :
268 : //! @endcond
|
