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