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