Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: gdal "raster shift-longitude" subcommand
5 : * Author: Dan Baston
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2026, ISciences LLC
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdalalg_raster_shift_longitude.h"
14 :
15 : #include "cpl_conv.h"
16 : #include "gdal_priv.h"
17 : #include "gdal_priv_templates.hpp"
18 : #include "../frmts/vrt/gdal_vrt.h"
19 : #include "../frmts/vrt/vrtdataset.h"
20 :
21 : #include <algorithm>
22 :
23 : //! @cond Doxygen_Suppress
24 :
25 : #ifndef _
26 : #define _(x) (x)
27 : #endif
28 :
29 : /************************************************************************/
30 : /*GDALRasterShiftLongitudeAlgorithm::GDALRasterShiftLongitudeAlgorithm()*/
31 : /************************************************************************/
32 :
33 99 : GDALRasterShiftLongitudeAlgorithm::GDALRasterShiftLongitudeAlgorithm(
34 99 : bool bStandalone)
35 99 : : GDALRasterPipelineStepAlgorithm(NAME, DESCRIPTION, HELP_URL, bStandalone)
36 : {
37 198 : AddArg("min-x", 0, _("Sets the minimum longitude value"), &m_minX)
38 99 : .SetRequired();
39 198 : AddArg("max-x", 0, _("Sets the maximum longitude value"), &m_maxX)
40 99 : .SetRequired();
41 99 : AddNodataArg(&m_nodata, /* noneAllowed = */ true, "output-nodata");
42 99 : }
43 :
44 : /************************************************************************/
45 : /* GDALRasterShiftLongitudeAlgorithm::RunStep() */
46 : /************************************************************************/
47 :
48 11 : static double NormalizeLongitude(double lon)
49 : {
50 11 : while (lon < -180.0)
51 3 : lon += 360.0;
52 8 : while (lon > 180.0)
53 0 : lon -= 360.0;
54 8 : return lon;
55 : }
56 :
57 19 : bool GDALRasterShiftLongitudeAlgorithm::RunStep(GDALPipelineStepRunContext &)
58 : {
59 19 : CPLAssert(!m_outputDataset.GetDatasetRef());
60 :
61 19 : const auto poSrcDS = m_inputDataset[0].GetDatasetRef();
62 :
63 19 : GDALGeoTransform srcGT;
64 19 : if (poSrcDS->GetGeoTransform(srcGT) != CE_None)
65 : {
66 1 : ReportError(CE_Failure, CPLE_AppDefined,
67 : "Input dataset does not have a geotransform");
68 1 : return false;
69 : }
70 :
71 18 : if (!srcGT.IsAxisAligned())
72 : {
73 1 : ReportError(CE_Failure, CPLE_AppDefined,
74 : "Input dataset geotransform cannot have a rotation term");
75 1 : return false;
76 : }
77 :
78 17 : const auto nSrcXSize = poSrcDS->GetRasterXSize();
79 :
80 17 : const double dfSrcMinX = srcGT.xorig;
81 17 : const double dfSrcMaxX = srcGT.xorig + srcGT.xscale * nSrcXSize;
82 :
83 : // Snap m_minX and m_maxX to pixel boundaries of the input raster
84 : {
85 : const double dfXOff =
86 34 : std::ceil(std::abs(dfSrcMinX - m_minX) / srcGT.xscale) *
87 17 : (m_minX < dfSrcMinX ? -1 : 1);
88 17 : if (!GDALIsValueInRange<int>(dfXOff))
89 : {
90 1 : CPLError(CE_Failure, CPLE_AppDefined,
91 : "Failed to compute pixel offset between minimum x "
92 : "coordinate of input and output.");
93 1 : return false;
94 : }
95 16 : const int nXOff = static_cast<int>(dfXOff);
96 16 : m_minX = dfSrcMinX + nXOff * srcGT.xscale;
97 : }
98 16 : const double dfDstXSize = std::ceil((m_maxX - m_minX) / srcGT.xscale);
99 16 : if (!GDALIsValueInRange<int>(dfDstXSize))
100 : {
101 1 : CPLError(CE_Failure, CPLE_AppDefined,
102 : "Size of output raster exceeds maximum allowable size.");
103 1 : return false;
104 : }
105 15 : const int nDstXSize = static_cast<int>(dfDstXSize);
106 :
107 15 : if (nDstXSize <= 0)
108 : {
109 2 : ReportError(CE_Failure, CPLE_AppDefined,
110 : "--max-x must be greater than --min-x");
111 2 : return false;
112 : }
113 :
114 13 : m_maxX = m_minX + nDstXSize * srcGT.xscale;
115 13 : const int nYSize = poSrcDS->GetRasterYSize();
116 :
117 26 : auto poDstDS = std::make_unique<VRTDataset>(nDstXSize, nYSize);
118 : {
119 13 : GDALGeoTransform dstGT = srcGT;
120 13 : dstGT.xorig = m_minX;
121 13 : poDstDS->SetGeoTransform(dstGT);
122 : }
123 13 : poDstDS->SetSpatialRef(poSrcDS->GetSpatialRef());
124 :
125 26 : std::vector<GDALRasterWindow> aosSrcWindows;
126 38 : for (int nDstXOff = 0; nDstXOff < nDstXSize;)
127 : {
128 25 : double dfDstChunkMinX = m_minX + nDstXOff * srcGT.xscale;
129 :
130 31 : while (dfDstChunkMinX < dfSrcMinX)
131 : {
132 6 : dfDstChunkMinX += 360;
133 : }
134 38 : while (dfDstChunkMinX >= dfSrcMaxX)
135 : {
136 13 : dfDstChunkMinX -= 360;
137 : }
138 :
139 : const double dfSrcXOff =
140 25 : std::round((dfDstChunkMinX - dfSrcMinX) / srcGT.xscale);
141 25 : if (!GDALIsValueInRange<int>(dfSrcXOff))
142 : {
143 0 : CPLError(CE_Failure, CPLE_AppDefined,
144 : "Failed to calculate source pixels for longitude range "
145 : "beginning at %g",
146 : dfDstChunkMinX);
147 0 : return false;
148 : }
149 25 : const int nSrcXOff = static_cast<int>(dfSrcXOff);
150 :
151 25 : if (nSrcXOff < 0)
152 : {
153 : const double dfMissingRangeMinX =
154 4 : NormalizeLongitude(dfDstChunkMinX);
155 : const int nMissingPx =
156 4 : std::min(std::abs(nSrcXOff), nDstXSize - nDstXOff);
157 : double dfMissingRangeMaxX =
158 4 : NormalizeLongitude(dfDstChunkMinX + srcGT.xscale * nMissingPx);
159 5 : while (dfMissingRangeMaxX < dfMissingRangeMinX)
160 : {
161 1 : dfMissingRangeMaxX += 360;
162 : }
163 :
164 4 : CPLError(
165 : CE_Warning, CPLE_AppDefined,
166 : "No source data available for output longitude range %g to %g",
167 : dfMissingRangeMinX, dfMissingRangeMaxX);
168 4 : nDstXOff += nMissingPx;
169 :
170 4 : const GDALRasterWindow chunk{-1, -1, nMissingPx, nYSize};
171 4 : aosSrcWindows.push_back(chunk);
172 4 : continue;
173 : }
174 :
175 21 : const int nColumnsWanted = nSrcXSize - nSrcXOff;
176 21 : const int nColumnsAvailable = nDstXSize - nDstXOff;
177 :
178 21 : const GDALRasterWindow chunk{
179 21 : nSrcXOff, 0, std::min(nColumnsWanted, nColumnsAvailable), nYSize};
180 :
181 21 : CPLAssert(chunk.nXSize > 0);
182 :
183 21 : const double dstChunkMaxX = dfDstChunkMinX + chunk.nXSize;
184 21 : CPLDebug("ShiftLongitude", "Src %g - %g, Dst %g - %g", dfDstChunkMinX,
185 21 : dstChunkMaxX, m_minX + chunk.nXOff * srcGT.xscale,
186 21 : m_minX + (chunk.nXOff + chunk.nXSize) * srcGT.xscale);
187 :
188 21 : aosSrcWindows.push_back(chunk);
189 21 : nDstXOff += chunk.nXSize;
190 : }
191 :
192 24 : for (int iBand = 1; iBand <= poSrcDS->GetRasterCount(); ++iBand)
193 : {
194 13 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(iBand);
195 13 : poDstDS->AddBand(poSrcBand->GetRasterDataType());
196 : VRTSourcedRasterBand *poDstBand =
197 13 : cpl::down_cast<VRTSourcedRasterBand *>(
198 13 : poDstDS->GetRasterBand(iBand));
199 13 : poDstBand->CopyCommonInfoFrom(poSrcBand);
200 :
201 : bool bHasNoData;
202 13 : double dfSrcNoData{0};
203 :
204 : {
205 : int bSuccess;
206 13 : dfSrcNoData = poSrcBand->GetNoDataValue(&bSuccess);
207 13 : if (!bSuccess)
208 : {
209 11 : dfSrcNoData = VRT_NODATA_UNSET;
210 : }
211 : }
212 :
213 13 : if (!m_nodata.empty())
214 : {
215 7 : bool inexact = false;
216 7 : if (poDstBand->SetNoDataValueAsString(m_nodata.c_str(), &inexact) !=
217 : CE_None)
218 : {
219 2 : CPLError(CE_Failure, CPLE_AppDefined,
220 : "Invalid NoData value: %s", m_nodata.c_str());
221 2 : return false;
222 : }
223 5 : if (inexact)
224 : {
225 0 : CPLError(CE_Failure, CPLE_AppDefined,
226 : "Specified NoData value cannot be represented in the "
227 : "output data type.");
228 0 : return false;
229 : }
230 5 : bHasNoData = true;
231 : }
232 : else
233 : {
234 6 : bHasNoData = GDALCopyNoDataValue(poDstBand, poSrcBand, nullptr);
235 : }
236 :
237 11 : int nDstXOff = 0;
238 34 : for (const auto &chunk : aosSrcWindows)
239 : {
240 23 : const bool bDataAvailable = chunk.nXOff >= 0;
241 :
242 23 : if (bDataAvailable)
243 : {
244 19 : CPLDebug("ShiftLongitude", "Adding source chunk %d,%d %d %d",
245 19 : chunk.nXOff, chunk.nYOff, chunk.nXSize, chunk.nYSize);
246 :
247 : // Scale and offset are propagated to the output band, not handled by the ComplexSource
248 19 : constexpr double dfSourceScale = 1.0;
249 19 : constexpr double dfSourceOffset = 0.0;
250 19 : constexpr int nSrcYOff = 0;
251 19 : constexpr int nDstYOff = 0;
252 38 : if (poDstBand->AddComplexSource(
253 19 : poSrcBand, chunk.nXOff, nSrcYOff, chunk.nXSize,
254 19 : chunk.nYSize, nDstXOff, nDstYOff, chunk.nXSize,
255 19 : chunk.nYSize, dfSourceOffset, dfSourceScale,
256 19 : dfSrcNoData) != CE_None)
257 : {
258 0 : return false;
259 : }
260 : }
261 4 : else if (!bHasNoData)
262 : {
263 0 : CPLErrorOnce(
264 : CE_Warning, CPLE_AppDefined,
265 : "Some regions of the output dataset have no source data "
266 : "available, but a NoData value has not been defined. You "
267 : "can set a value using --output-nodata");
268 : }
269 23 : nDstXOff += chunk.nXSize;
270 : }
271 : }
272 :
273 11 : m_outputDataset.Set(std::move(poDstDS));
274 :
275 11 : return true;
276 : }
277 :
278 : GDALRasterShiftLongitudeAlgorithmStandalone::
279 : ~GDALRasterShiftLongitudeAlgorithmStandalone() = default;
280 :
281 : //! @endcond
|
