Line data Source code
1 : /******************************************************************************
2 : * Name: gdalmultiDim_meshgrid.cpp
3 : * Project: GDAL Core
4 : * Purpose: Return a vector of coordinate matrices from coordinate vectors.
5 : * Author: Even Rouault <even.rouault at spatialys.com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2024, Even Rouault <even.rouault at spatialys.com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "gdal_priv.h"
14 :
15 : #include <algorithm>
16 : #include <limits>
17 :
18 : /************************************************************************/
19 : /* GetConcatenatedNames() */
20 : /************************************************************************/
21 :
22 : static std::string
23 20 : GetConcatenatedNames(const std::vector<std::shared_ptr<GDALMDArray>> &apoArrays)
24 : {
25 20 : std::string ret;
26 72 : for (const auto &poArray : apoArrays)
27 : {
28 52 : if (!ret.empty())
29 32 : ret += ", ";
30 52 : ret += poArray->GetFullName();
31 : }
32 20 : return ret;
33 : }
34 :
35 : /************************************************************************/
36 : /* GDALMDArrayMeshGrid */
37 : /************************************************************************/
38 :
39 : class GDALMDArrayMeshGrid final : public GDALMDArray
40 : {
41 : const std::vector<std::shared_ptr<GDALMDArray>> m_apoArrays;
42 : std::vector<std::shared_ptr<GDALDimension>> m_apoDims{};
43 : const size_t m_iDim;
44 : const bool m_bIJIndexing;
45 :
46 : protected:
47 10 : explicit GDALMDArrayMeshGrid(
48 : const std::vector<std::shared_ptr<GDALMDArray>> &apoArrays,
49 : const std::vector<std::shared_ptr<GDALDimension>> &apoDims, size_t iDim,
50 : bool bIJIndexing)
51 20 : : GDALAbstractMDArray(std::string(),
52 20 : "Mesh grid view of " +
53 10 : GetConcatenatedNames(apoArrays)),
54 20 : GDALMDArray(std::string(),
55 20 : "Mesh grid view of " + GetConcatenatedNames(apoArrays)),
56 : m_apoArrays(apoArrays), m_apoDims(apoDims), m_iDim(iDim),
57 50 : m_bIJIndexing(bIJIndexing)
58 : {
59 10 : }
60 :
61 : bool IRead(const GUInt64 *arrayStartIdx, const size_t *count,
62 : const GInt64 *arrayStep, const GPtrDiff_t *bufferStride,
63 : const GDALExtendedDataType &bufferDataType,
64 : void *pDstBuffer) const override;
65 :
66 : public:
67 : static std::shared_ptr<GDALMDArrayMeshGrid>
68 10 : Create(const std::vector<std::shared_ptr<GDALMDArray>> &apoArrays,
69 : size_t iDim, bool bIJIndexing)
70 : {
71 20 : std::vector<std::shared_ptr<GDALDimension>> apoDims;
72 36 : for (size_t i = 0; i < apoArrays.size(); ++i)
73 : {
74 26 : const size_t iTranslatedDim = (!bIJIndexing && i <= 1) ? 1 - i : i;
75 26 : apoDims.push_back(apoArrays[iTranslatedDim]->GetDimensions()[0]);
76 : }
77 : auto newAr(std::shared_ptr<GDALMDArrayMeshGrid>(
78 10 : new GDALMDArrayMeshGrid(apoArrays, apoDims, iDim, bIJIndexing)));
79 10 : newAr->SetSelf(newAr);
80 20 : return newAr;
81 : }
82 :
83 1 : bool IsWritable() const override
84 : {
85 1 : return false;
86 : }
87 :
88 13 : const std::string &GetFilename() const override
89 : {
90 13 : return m_apoArrays[m_iDim]->GetFilename();
91 : }
92 :
93 : const std::vector<std::shared_ptr<GDALDimension>> &
94 165 : GetDimensions() const override
95 : {
96 165 : return m_apoDims;
97 : }
98 :
99 60 : const GDALExtendedDataType &GetDataType() const override
100 : {
101 60 : return m_apoArrays[m_iDim]->GetDataType();
102 : }
103 :
104 : std::shared_ptr<GDALAttribute>
105 1 : GetAttribute(const std::string &osName) const override
106 : {
107 1 : return m_apoArrays[m_iDim]->GetAttribute(osName);
108 : }
109 :
110 : std::vector<std::shared_ptr<GDALAttribute>>
111 2 : GetAttributes(CSLConstList papszOptions = nullptr) const override
112 : {
113 2 : return m_apoArrays[m_iDim]->GetAttributes(papszOptions);
114 : }
115 :
116 1 : const std::string &GetUnit() const override
117 : {
118 1 : return m_apoArrays[m_iDim]->GetUnit();
119 : }
120 :
121 1 : const void *GetRawNoDataValue() const override
122 : {
123 1 : return m_apoArrays[m_iDim]->GetRawNoDataValue();
124 : }
125 :
126 1 : double GetOffset(bool *pbHasOffset,
127 : GDALDataType *peStorageType) const override
128 : {
129 1 : return m_apoArrays[m_iDim]->GetOffset(pbHasOffset, peStorageType);
130 : }
131 :
132 1 : double GetScale(bool *pbHasScale,
133 : GDALDataType *peStorageType) const override
134 : {
135 1 : return m_apoArrays[m_iDim]->GetScale(pbHasScale, peStorageType);
136 : }
137 : };
138 :
139 : /************************************************************************/
140 : /* IRead() */
141 : /************************************************************************/
142 :
143 20 : bool GDALMDArrayMeshGrid::IRead(const GUInt64 *arrayStartIdx,
144 : const size_t *count, const GInt64 *arrayStep,
145 : const GPtrDiff_t *bufferStride,
146 : const GDALExtendedDataType &bufferDataType,
147 : void *pDstBuffer) const
148 : {
149 20 : const size_t nBufferDTSize = bufferDataType.GetSize();
150 20 : const size_t iTranslatedDim =
151 20 : (!m_bIJIndexing && m_iDim <= 1) ? 1 - m_iDim : m_iDim;
152 40 : std::vector<GByte> abyTmpData(nBufferDTSize * count[iTranslatedDim]);
153 20 : const GPtrDiff_t strideOne[] = {1};
154 40 : if (!m_apoArrays[m_iDim]->Read(&arrayStartIdx[iTranslatedDim],
155 20 : &count[iTranslatedDim],
156 20 : &arrayStep[iTranslatedDim], strideOne,
157 20 : bufferDataType, abyTmpData.data()))
158 0 : return false;
159 :
160 20 : const auto nDims = GetDimensionCount();
161 :
162 : struct Stack
163 : {
164 : size_t nIters = 0;
165 : GByte *dst_ptr = nullptr;
166 : GPtrDiff_t dst_inc_offset = 0;
167 : };
168 :
169 : // +1 to avoid -Werror=null-dereference
170 20 : std::vector<Stack> stack(nDims + 1);
171 72 : for (size_t i = 0; i < nDims; i++)
172 : {
173 52 : stack[i].dst_inc_offset =
174 52 : static_cast<GPtrDiff_t>(bufferStride[i] * nBufferDTSize);
175 : }
176 20 : stack[0].dst_ptr = static_cast<GByte *>(pDstBuffer);
177 20 : size_t dimIdx = 0;
178 20 : size_t valIdx = 0;
179 :
180 142 : lbl_next_depth:
181 142 : if (dimIdx == nDims - 1)
182 : {
183 92 : auto nIters = count[dimIdx];
184 92 : GByte *dst_ptr = stack[dimIdx].dst_ptr;
185 92 : if (dimIdx == iTranslatedDim)
186 : {
187 34 : valIdx = 0;
188 : while (true)
189 : {
190 120 : GDALExtendedDataType::CopyValue(
191 120 : &abyTmpData[nBufferDTSize * valIdx], bufferDataType,
192 : dst_ptr, bufferDataType);
193 120 : if ((--nIters) == 0)
194 34 : break;
195 86 : ++valIdx;
196 86 : dst_ptr += stack[dimIdx].dst_inc_offset;
197 : }
198 : }
199 : else
200 : {
201 : while (true)
202 : {
203 216 : GDALExtendedDataType::CopyValue(
204 216 : &abyTmpData[nBufferDTSize * valIdx], bufferDataType,
205 : dst_ptr, bufferDataType);
206 216 : if ((--nIters) == 0)
207 58 : break;
208 158 : dst_ptr += stack[dimIdx].dst_inc_offset;
209 : }
210 : }
211 : }
212 : else
213 : {
214 50 : if (dimIdx == iTranslatedDim)
215 18 : valIdx = 0;
216 50 : stack[dimIdx].nIters = count[dimIdx];
217 : while (true)
218 : {
219 122 : dimIdx++;
220 122 : stack[dimIdx].dst_ptr = stack[dimIdx - 1].dst_ptr;
221 122 : goto lbl_next_depth;
222 122 : lbl_return_to_caller:
223 122 : dimIdx--;
224 122 : if ((--stack[dimIdx].nIters) == 0)
225 50 : break;
226 72 : if (dimIdx == iTranslatedDim)
227 26 : ++valIdx;
228 72 : stack[dimIdx].dst_ptr += stack[dimIdx].dst_inc_offset;
229 : }
230 : }
231 142 : if (dimIdx > 0)
232 : {
233 122 : goto lbl_return_to_caller;
234 : }
235 :
236 20 : if (bufferDataType.NeedsFreeDynamicMemory())
237 : {
238 20 : for (size_t i = 0; i < count[iTranslatedDim]; ++i)
239 : {
240 16 : bufferDataType.FreeDynamicMemory(&abyTmpData[i * nBufferDTSize]);
241 : }
242 : }
243 :
244 20 : return true;
245 : }
246 :
247 : /************************************************************************/
248 : /* GDALMDArrayGetMeshGrid() */
249 : /************************************************************************/
250 :
251 : /** Return a list of multidimensional arrays from a list of one-dimensional
252 : * arrays.
253 : *
254 : * This is typically used to transform one-dimensional longitude, latitude
255 : * arrays into 2D ones.
256 : *
257 : * More formally, for one-dimensional arrays x1, x2,..., xn with lengths
258 : * Ni=len(xi), returns (N1, N2, ..., Nn) shaped arrays if indexing="ij" or
259 : * (N2, N1, ..., Nn) shaped arrays if indexing="xy" with the elements of xi
260 : * repeated to fill the matrix along the first dimension for x1, the second
261 : * for x2 and so on.
262 : *
263 : * For example, if x = [1, 2], and y = [3, 4, 5],
264 : * GetMeshGrid([x, y], ["INDEXING=xy"]) will return [xm, ym] such that
265 : * xm=[[1, 2],[1, 2],[1, 2]] and ym=[[3, 3],[4, 4],[5, 5]],
266 : * or more generally xm[any index][i] = x[i] and ym[i][any index]=y[i]
267 : *
268 : * and
269 : * GetMeshGrid([x, y], ["INDEXING=ij"]) will return [xm, ym] such that
270 : * xm=[[1, 1, 1],[2, 2, 2]] and ym=[[3, 4, 5],[3, 4, 5]],
271 : * or more generally xm[i][any index] = x[i] and ym[any index][i]=y[i]
272 : *
273 : * The currently supported options are:
274 : * <ul>
275 : * <li>INDEXING=xy/ij: Cartesian ("xy", default) or matrix ("ij") indexing of
276 : * output.
277 : * </li>
278 : * </ul>
279 : *
280 : * This is the same as
281 : * <a href="https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html">numpy.meshgrid()</a>
282 : * function.
283 : *
284 : * This is the same as the C function GDALMDArrayGetMeshGrid()
285 : *
286 : * @param apoArrays Input arrays
287 : * @param papszOptions NULL, or NULL terminated list of options.
288 : *
289 : * @return an array of coordinate matrices
290 : * @since 3.10
291 : */
292 :
293 7 : /* static */ std::vector<std::shared_ptr<GDALMDArray>> GDALMDArray::GetMeshGrid(
294 : const std::vector<std::shared_ptr<GDALMDArray>> &apoArrays,
295 : CSLConstList papszOptions)
296 : {
297 7 : std::vector<std::shared_ptr<GDALMDArray>> ret;
298 19 : for (const auto &poArray : apoArrays)
299 : {
300 13 : if (poArray->GetDimensionCount() != 1)
301 : {
302 1 : CPLError(CE_Failure, CPLE_NotSupported,
303 : "Only 1-D input arrays are accepted");
304 1 : return ret;
305 : }
306 : }
307 :
308 : const char *pszIndexing =
309 6 : CSLFetchNameValueDef(papszOptions, "INDEXING", "xy");
310 6 : if (!EQUAL(pszIndexing, "xy") && !EQUAL(pszIndexing, "ij"))
311 : {
312 1 : CPLError(CE_Failure, CPLE_NotSupported,
313 : "Only INDEXING=xy or ij is accepted");
314 1 : return ret;
315 : }
316 5 : const bool bIJIndexing = EQUAL(pszIndexing, "ij");
317 :
318 15 : for (size_t i = 0; i < apoArrays.size(); ++i)
319 : {
320 10 : ret.push_back(GDALMDArrayMeshGrid::Create(apoArrays, i, bIJIndexing));
321 : }
322 :
323 5 : return ret;
324 : }
|
