Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Mapinfo Image Warper
4 : * Purpose: Implementation of one or more GDALTrasformerFunc types, including
5 : * the GenImgProj (general image reprojector) transformer.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2002, i3 - information integration and imaging
10 : * Fort Collin, CO
11 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
12 : * Copyright (c) 2021, CLS
13 : *
14 : * SPDX-License-Identifier: MIT
15 : ****************************************************************************/
16 :
17 : #include "cpl_port.h"
18 : #include "gdal_alg.h"
19 : #include "gdal_alg_priv.h"
20 :
21 : #include <climits>
22 : #include <cmath>
23 : #include <cstddef>
24 : #include <cstdlib>
25 : #include <cstring>
26 :
27 : #include <algorithm>
28 : #include <limits>
29 : #include <utility>
30 :
31 : #include "cpl_conv.h"
32 : #include "cpl_error.h"
33 : #include "cpl_list.h"
34 : #include "cpl_minixml.h"
35 : #include "cpl_multiproc.h"
36 : #include "cpl_string.h"
37 : #include "cpl_vsi.h"
38 : #include "gdal.h"
39 : #include "gdal_priv.h"
40 : #include "ogr_core.h"
41 : #include "ogr_spatialref.h"
42 : #include "ogr_srs_api.h"
43 :
44 : CPL_C_START
45 : void *GDALDeserializeGCPTransformer(CPLXMLNode *psTree);
46 : void *GDALDeserializeTPSTransformer(CPLXMLNode *psTree);
47 : void *GDALDeserializeGeoLocTransformer(CPLXMLNode *psTree);
48 : void *GDALDeserializeRPCTransformer(CPLXMLNode *psTree);
49 : CPL_C_END
50 :
51 : static CPLXMLNode *GDALSerializeReprojectionTransformer(void *pTransformArg);
52 : static void *GDALDeserializeReprojectionTransformer(CPLXMLNode *psTree);
53 :
54 : static CPLXMLNode *GDALSerializeGenImgProjTransformer(void *pTransformArg);
55 : static void *GDALDeserializeGenImgProjTransformer(CPLXMLNode *psTree);
56 :
57 : static void *GDALCreateApproxTransformer2(GDALTransformerFunc pfnRawTransformer,
58 : void *pRawTransformerArg,
59 : double dfMaxErrorForward,
60 : double dfMaxErrorReverse);
61 :
62 : /************************************************************************/
63 : /* GDALIsTransformer() */
64 : /************************************************************************/
65 :
66 10598 : bool GDALIsTransformer(void *hTransformerArg, const char *pszClassName)
67 : {
68 10598 : if (!hTransformerArg)
69 0 : return false;
70 : // All transformers should have a GDALTransformerInfo member as their first members
71 10598 : GDALTransformerInfo *psInfo =
72 : static_cast<GDALTransformerInfo *>(hTransformerArg);
73 10598 : return memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
74 20302 : strlen(GDAL_GTI2_SIGNATURE)) == 0 &&
75 20302 : strcmp(psInfo->pszClassName, pszClassName) == 0;
76 : }
77 :
78 : /************************************************************************/
79 : /* GDALTransformFunc */
80 : /* */
81 : /* Documentation for GDALTransformFunc typedef. */
82 : /************************************************************************/
83 :
84 : /*!
85 :
86 : \typedef typedef int (*GDALTransformerFunc)( void *pTransformerArg, int
87 : bDstToSrc, int nPointCount, double *x, double *y, double *z, int *panSuccess );
88 :
89 : Generic signature for spatial point transformers.
90 :
91 : This function signature is used for a variety of functions that accept
92 : passed in functions used to transform point locations between two coordinate
93 : spaces.
94 :
95 : The GDALCreateGenImgProjTransformer(), GDALCreateReprojectionTransformerEx(),
96 : GDALCreateGCPTransformer() and GDALCreateApproxTransformer() functions can
97 : be used to prepare argument data for some built-in transformers. As well,
98 : applications can implement their own transformers to the following signature.
99 :
100 : \code
101 : typedef int
102 : (*GDALTransformerFunc)( void *pTransformerArg,
103 : int bDstToSrc, int nPointCount,
104 : double *x, double *y, double *z, int *panSuccess );
105 : \endcode
106 :
107 : @param pTransformerArg application supplied callback data used by the
108 : transformer.
109 :
110 : @param bDstToSrc if TRUE the transformation will be from the destination
111 : coordinate space to the source coordinate system, otherwise the transformation
112 : will be from the source coordinate system to the destination coordinate system.
113 :
114 : @param nPointCount number of points in the x, y and z arrays.
115 :
116 : @param x input X coordinates. Results returned in same array.
117 :
118 : @param y input Y coordinates. Results returned in same array.
119 :
120 : @param z input Z coordinates. Results returned in same array.
121 :
122 : @param panSuccess array of ints in which success (TRUE) or failure (FALSE)
123 : flags are returned for the translation of each point.
124 :
125 : @return TRUE if the overall transformation succeeds (though some individual
126 : points may have failed) or FALSE if the overall transformation fails.
127 :
128 : */
129 :
130 : /************************************************************************/
131 : /* GDALSuggestedWarpOutput() */
132 : /************************************************************************/
133 :
134 : /**
135 : * Suggest output file size.
136 : *
137 : * This function is used to suggest the size, and georeferenced extents
138 : * appropriate given the indicated transformation and input file. It walks
139 : * the edges of the input file (approximately 20 sample points along each
140 : * edge) transforming into output coordinates in order to get an extents box.
141 : *
142 : * Then a resolution is computed with the intent that the length of the
143 : * distance from the top left corner of the output imagery to the bottom right
144 : * corner would represent the same number of pixels as in the source image.
145 : * Note that if the image is somewhat rotated the diagonal taken isn't of the
146 : * whole output bounding rectangle, but instead of the locations where the
147 : * top/left and bottom/right corners transform. The output pixel size is
148 : * always square. This is intended to approximately preserve the resolution
149 : * of the input data in the output file.
150 : *
151 : * The values returned in padfGeoTransformOut, pnPixels and pnLines are
152 : * the suggested number of pixels and lines for the output file, and the
153 : * geotransform relating those pixels to the output georeferenced coordinates.
154 : *
155 : * The trickiest part of using the function is ensuring that the
156 : * transformer created is from source file pixel/line coordinates to
157 : * output file georeferenced coordinates. This can be accomplished with
158 : * GDALCreateGenImgProjTransformer() by passing a NULL for the hDstDS.
159 : *
160 : * @param hSrcDS the input image (it is assumed the whole input image is
161 : * being transformed).
162 : * @param pfnTransformer the transformer function.
163 : * @param pTransformArg the callback data for the transformer function.
164 : * @param padfGeoTransformOut the array of six doubles in which the suggested
165 : * geotransform is returned.
166 : * @param pnPixels int in which the suggest pixel width of output is returned.
167 : * @param pnLines int in which the suggest pixel height of output is returned.
168 : *
169 : * @return CE_None if successful or CE_Failure otherwise.
170 : */
171 :
172 40 : CPLErr CPL_STDCALL GDALSuggestedWarpOutput(GDALDatasetH hSrcDS,
173 : GDALTransformerFunc pfnTransformer,
174 : void *pTransformArg,
175 : double *padfGeoTransformOut,
176 : int *pnPixels, int *pnLines)
177 :
178 : {
179 40 : VALIDATE_POINTER1(hSrcDS, "GDALSuggestedWarpOutput", CE_Failure);
180 :
181 40 : double adfExtent[4] = {};
182 :
183 40 : return GDALSuggestedWarpOutput2(hSrcDS, pfnTransformer, pTransformArg,
184 : padfGeoTransformOut, pnPixels, pnLines,
185 40 : adfExtent, 0);
186 : }
187 :
188 488 : static bool GDALSuggestedWarpOutput2_MustAdjustForRightBorder(
189 : GDALTransformerFunc pfnTransformer, void *pTransformArg, double *padfExtent,
190 : int /* nPixels*/, int nLines, double dfPixelSizeX, double dfPixelSizeY)
191 : {
192 488 : double adfX[21] = {};
193 488 : double adfY[21] = {};
194 :
195 488 : const double dfMaxXOut = padfExtent[2];
196 488 : const double dfMaxYOut = padfExtent[3];
197 :
198 : // Take 20 steps.
199 488 : int nSamplePoints = 0;
200 10736 : for (double dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05)
201 : {
202 : // Ensure we end exactly at the end.
203 10248 : if (dfRatio > 0.99)
204 488 : dfRatio = 1.0;
205 :
206 : // Along right.
207 10248 : adfX[nSamplePoints] = dfMaxXOut;
208 10248 : adfY[nSamplePoints] = dfMaxYOut - dfPixelSizeY * dfRatio * nLines;
209 10248 : nSamplePoints++;
210 : }
211 488 : double adfZ[21] = {};
212 :
213 488 : int abSuccess[21] = {};
214 :
215 488 : pfnTransformer(pTransformArg, TRUE, nSamplePoints, adfX, adfY, adfZ,
216 : abSuccess);
217 :
218 488 : int abSuccess2[21] = {};
219 :
220 488 : pfnTransformer(pTransformArg, FALSE, nSamplePoints, adfX, adfY, adfZ,
221 : abSuccess2);
222 :
223 488 : nSamplePoints = 0;
224 488 : int nBadCount = 0;
225 10736 : for (double dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05)
226 : {
227 10248 : const double expected_x = dfMaxXOut;
228 10248 : const double expected_y = dfMaxYOut - dfPixelSizeY * dfRatio * nLines;
229 10248 : if (!abSuccess[nSamplePoints] || !abSuccess2[nSamplePoints] ||
230 7644 : fabs(adfX[nSamplePoints] - expected_x) > dfPixelSizeX ||
231 6018 : fabs(adfY[nSamplePoints] - expected_y) > dfPixelSizeY)
232 : {
233 4230 : nBadCount++;
234 : }
235 10248 : nSamplePoints++;
236 : }
237 :
238 488 : return nBadCount == nSamplePoints;
239 : }
240 :
241 387 : static bool GDALSuggestedWarpOutput2_MustAdjustForBottomBorder(
242 : GDALTransformerFunc pfnTransformer, void *pTransformArg, double *padfExtent,
243 : int nPixels, int /* nLines */, double dfPixelSizeX, double dfPixelSizeY)
244 : {
245 387 : double adfX[21] = {};
246 387 : double adfY[21] = {};
247 :
248 387 : const double dfMinXOut = padfExtent[0];
249 387 : const double dfMinYOut = padfExtent[1];
250 :
251 : // Take 20 steps.
252 387 : int nSamplePoints = 0;
253 8514 : for (double dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05)
254 : {
255 : // Ensure we end exactly at the end.
256 8127 : if (dfRatio > 0.99)
257 387 : dfRatio = 1.0;
258 :
259 : // Along right.
260 8127 : adfX[nSamplePoints] = dfMinXOut + dfPixelSizeX * dfRatio * nPixels;
261 8127 : adfY[nSamplePoints] = dfMinYOut;
262 8127 : nSamplePoints++;
263 : }
264 387 : double adfZ[21] = {};
265 :
266 387 : int abSuccess[21] = {};
267 :
268 387 : pfnTransformer(pTransformArg, TRUE, nSamplePoints, adfX, adfY, adfZ,
269 : abSuccess);
270 :
271 387 : int abSuccess2[21] = {};
272 :
273 387 : pfnTransformer(pTransformArg, FALSE, nSamplePoints, adfX, adfY, adfZ,
274 : abSuccess2);
275 :
276 387 : nSamplePoints = 0;
277 387 : int nBadCount = 0;
278 8514 : for (double dfRatio = 0.0; dfRatio <= 1.01; dfRatio += 0.05)
279 : {
280 8127 : const double expected_x = dfMinXOut + dfPixelSizeX * dfRatio * nPixels;
281 8127 : const double expected_y = dfMinYOut;
282 8127 : if (!abSuccess[nSamplePoints] || !abSuccess2[nSamplePoints] ||
283 6241 : fabs(adfX[nSamplePoints] - expected_x) > dfPixelSizeX ||
284 5952 : fabs(adfY[nSamplePoints] - expected_y) > dfPixelSizeY)
285 : {
286 2175 : nBadCount++;
287 : }
288 8127 : nSamplePoints++;
289 : }
290 :
291 387 : return nBadCount == nSamplePoints;
292 : }
293 :
294 : /************************************************************************/
295 : /* GDALSuggestedWarpOutput2() */
296 : /************************************************************************/
297 :
298 : /**
299 : * Suggest output file size.
300 : *
301 : * This function is used to suggest the size, and georeferenced extents
302 : * appropriate given the indicated transformation and input file. It walks
303 : * the edges of the input file (approximately 20 sample points along each
304 : * edge) transforming into output coordinates in order to get an extents box.
305 : *
306 : * Then a resolution is computed with the intent that the length of the
307 : * distance from the top left corner of the output imagery to the bottom right
308 : * corner would represent the same number of pixels as in the source image.
309 : * Note that if the image is somewhat rotated the diagonal taken isn't of the
310 : * whole output bounding rectangle, but instead of the locations where the
311 : * top/left and bottom/right corners transform. The output pixel size is
312 : * always square. This is intended to approximately preserve the resolution
313 : * of the input data in the output file.
314 : *
315 : * The values returned in padfGeoTransformOut, pnPixels and pnLines are
316 : * the suggested number of pixels and lines for the output file, and the
317 : * geotransform relating those pixels to the output georeferenced coordinates.
318 : *
319 : * The trickiest part of using the function is ensuring that the
320 : * transformer created is from source file pixel/line coordinates to
321 : * output file georeferenced coordinates. This can be accomplished with
322 : * GDALCreateGenImgProjTransformer() by passing a NULL for the hDstDS.
323 : *
324 : * @param hSrcDS the input image (it is assumed the whole input image is
325 : * being transformed).
326 : * @param pfnTransformer the transformer function.
327 : * @param pTransformArg the callback data for the transformer function.
328 : * @param padfGeoTransformOut the array of six doubles in which the suggested
329 : * geotransform is returned.
330 : * @param pnPixels int in which the suggest pixel width of output is returned.
331 : * @param pnLines int in which the suggest pixel height of output is returned.
332 : * @param padfExtent Four entry array to return extents as (xmin, ymin, xmax,
333 : * ymax).
334 : * @param nOptions Options flags. Zero or GDAL_SWO_ROUND_UP_SIZE to ask *pnPixels
335 : * and *pnLines to be rounded up instead of being rounded to the closes integer, or
336 : * GDAL_SWO_FORCE_SQUARE_PIXEL to indicate that the generated pixel size is a square.
337 : *
338 : * @return CE_None if successful or CE_Failure otherwise.
339 : */
340 :
341 916 : CPLErr CPL_STDCALL GDALSuggestedWarpOutput2(GDALDatasetH hSrcDS,
342 : GDALTransformerFunc pfnTransformer,
343 : void *pTransformArg,
344 : double *padfGeoTransformOut,
345 : int *pnPixels, int *pnLines,
346 : double *padfExtent, int nOptions)
347 : {
348 916 : VALIDATE_POINTER1(hSrcDS, "GDALSuggestedWarpOutput2", CE_Failure);
349 :
350 : const bool bIsGDALGenImgProjTransform{
351 1832 : pTransformArg &&
352 916 : GDALIsTransformer(pTransformArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME)};
353 :
354 : /* -------------------------------------------------------------------- */
355 : /* Setup sample points all around the edge of the input raster. */
356 : /* -------------------------------------------------------------------- */
357 916 : if (bIsGDALGenImgProjTransform)
358 : {
359 : // In case CHECK_WITH_INVERT_PROJ has been modified.
360 916 : GDALRefreshGenImgProjTransformer(pTransformArg);
361 : }
362 0 : else if (GDALIsTransformer(pTransformArg,
363 : GDAL_APPROX_TRANSFORMER_CLASS_NAME))
364 : {
365 : // In case CHECK_WITH_INVERT_PROJ has been modified.
366 0 : GDALRefreshApproxTransformer(pTransformArg);
367 : }
368 :
369 916 : const int nInXSize = GDALGetRasterXSize(hSrcDS);
370 916 : const int nInYSize = GDALGetRasterYSize(hSrcDS);
371 :
372 : /* ------------------------------------------------------------- */
373 : /* Special case for warping on the same (or null) CRS. */
374 : /* ------------------------------------------------------------- */
375 916 : if ((!nOptions || (nOptions & GDAL_SWO_FORCE_SQUARE_PIXEL) == 0) &&
376 915 : pTransformArg && bIsGDALGenImgProjTransform)
377 : {
378 915 : const GDALGenImgProjTransformInfo *psInfo =
379 : static_cast<const GDALGenImgProjTransformInfo *>(pTransformArg);
380 :
381 915 : if (!psInfo->pSrcTransformer &&
382 842 : !psInfo->bHasCustomTransformationPipeline &&
383 838 : !psInfo->pDstTransformer && psInfo->adfSrcGeoTransform[2] == 0 &&
384 838 : psInfo->adfSrcGeoTransform[4] == 0 &&
385 838 : psInfo->adfDstGeoTransform[0] == 0 &&
386 826 : psInfo->adfDstGeoTransform[1] == 1 &&
387 826 : psInfo->adfDstGeoTransform[2] == 0 &&
388 826 : psInfo->adfDstGeoTransform[3] == 0 &&
389 826 : psInfo->adfDstGeoTransform[4] == 0 &&
390 826 : psInfo->adfDstGeoTransform[5] == 1)
391 : {
392 826 : const OGRSpatialReference *poSourceCRS = nullptr;
393 826 : const OGRSpatialReference *poTargetCRS = nullptr;
394 :
395 826 : if (psInfo->pReprojectArg)
396 : {
397 576 : const GDALReprojectionTransformInfo *psRTI =
398 : static_cast<const GDALReprojectionTransformInfo *>(
399 : psInfo->pReprojectArg);
400 576 : poSourceCRS = psRTI->poForwardTransform->GetSourceCS();
401 576 : poTargetCRS = psRTI->poForwardTransform->GetTargetCS();
402 : }
403 :
404 1402 : if ((!poSourceCRS && !poTargetCRS) ||
405 576 : (poSourceCRS && poTargetCRS &&
406 576 : poSourceCRS->IsSame(poTargetCRS)))
407 : {
408 :
409 595 : const bool bNorthUp{psInfo->adfSrcGeoTransform[5] < 0.0};
410 :
411 595 : memcpy(padfGeoTransformOut, psInfo->adfSrcGeoTransform,
412 : sizeof(double) * 6);
413 :
414 595 : if (!bNorthUp)
415 : {
416 58 : padfGeoTransformOut[3] = padfGeoTransformOut[3] +
417 58 : nInYSize * padfGeoTransformOut[5];
418 58 : padfGeoTransformOut[5] = -padfGeoTransformOut[5];
419 : }
420 :
421 595 : *pnPixels = nInXSize;
422 595 : *pnLines = nInYSize;
423 :
424 : // Calculate extent from hSrcDS
425 595 : if (padfExtent)
426 : {
427 595 : padfExtent[0] = psInfo->adfSrcGeoTransform[0];
428 595 : padfExtent[1] = psInfo->adfSrcGeoTransform[3] +
429 595 : nInYSize * psInfo->adfSrcGeoTransform[5];
430 595 : padfExtent[2] = psInfo->adfSrcGeoTransform[0] +
431 595 : nInXSize * psInfo->adfSrcGeoTransform[1];
432 595 : padfExtent[3] = psInfo->adfSrcGeoTransform[3];
433 595 : if (!bNorthUp)
434 : {
435 58 : std::swap(padfExtent[1], padfExtent[3]);
436 : }
437 : }
438 595 : return CE_None;
439 : }
440 : }
441 : }
442 :
443 321 : const int N_PIXELSTEP = 50;
444 : int nSteps = static_cast<int>(
445 321 : static_cast<double>(std::min(nInYSize, nInXSize)) / N_PIXELSTEP + 0.5);
446 321 : if (nSteps < 20)
447 296 : nSteps = 20;
448 25 : else if (nSteps > 100)
449 10 : nSteps = 100;
450 :
451 : // TODO(rouault): How is this goto retry supposed to work? Added in r20537.
452 : // Does redoing the same malloc multiple times work? If it is needed, can
453 : // it be converted to a tigher while loop around the MALLOC3s and free? Is
454 : // the point to try with the full requested steps. Then, if there is not
455 : // enough memory, back off and try with just 20 steps?
456 321 : retry:
457 321 : int nStepsPlusOne = nSteps + 1;
458 321 : int nSampleMax = nStepsPlusOne * nStepsPlusOne;
459 :
460 321 : double dfStep = 1.0 / nSteps;
461 321 : double *padfY = nullptr;
462 321 : double *padfZ = nullptr;
463 321 : double *padfYRevert = nullptr;
464 321 : double *padfZRevert = nullptr;
465 :
466 : int *pabSuccess = static_cast<int *>(
467 321 : VSI_MALLOC3_VERBOSE(sizeof(int), nStepsPlusOne, nStepsPlusOne));
468 : double *padfX = static_cast<double *>(
469 321 : VSI_MALLOC3_VERBOSE(sizeof(double) * 3, nStepsPlusOne, nStepsPlusOne));
470 : double *padfXRevert = static_cast<double *>(
471 321 : VSI_MALLOC3_VERBOSE(sizeof(double) * 3, nStepsPlusOne, nStepsPlusOne));
472 321 : if (pabSuccess == nullptr || padfX == nullptr || padfXRevert == nullptr)
473 : {
474 0 : CPLFree(padfX);
475 0 : CPLFree(padfXRevert);
476 0 : CPLFree(pabSuccess);
477 0 : if (nSteps > 20)
478 : {
479 0 : nSteps = 20;
480 0 : goto retry;
481 : }
482 0 : return CE_Failure;
483 : }
484 :
485 321 : padfY = padfX + nSampleMax;
486 321 : padfZ = padfX + nSampleMax * 2;
487 321 : padfYRevert = padfXRevert + nSampleMax;
488 321 : padfZRevert = padfXRevert + nSampleMax * 2;
489 :
490 : // Take N_STEPS steps.
491 7965 : for (int iStep = 0; iStep <= nSteps; iStep++)
492 : {
493 7644 : double dfRatio = (iStep == nSteps) ? 1.0 : iStep * dfStep;
494 7644 : int iStep2 = iStep;
495 :
496 : // Along top.
497 7644 : padfX[iStep2] = dfRatio * nInXSize;
498 7644 : padfY[iStep2] = 0.0;
499 7644 : padfZ[iStep2] = 0.0;
500 :
501 : // Along bottom.
502 7644 : iStep2 += nStepsPlusOne;
503 7644 : padfX[iStep2] = dfRatio * nInXSize;
504 7644 : padfY[iStep2] = nInYSize;
505 7644 : padfZ[iStep2] = 0.0;
506 :
507 : // Along left.
508 7644 : iStep2 += nStepsPlusOne;
509 7644 : padfX[iStep2] = 0.0;
510 7644 : padfY[iStep2] = dfRatio * nInYSize;
511 7644 : padfZ[iStep2] = 0.0;
512 :
513 : // Along right.
514 7644 : iStep2 += nStepsPlusOne;
515 7644 : padfX[iStep2] = nInXSize;
516 7644 : padfY[iStep2] = dfRatio * nInYSize;
517 7644 : padfZ[iStep2] = 0.0;
518 : }
519 :
520 321 : int nSamplePoints = 4 * nStepsPlusOne;
521 :
522 321 : memset(pabSuccess, 1, sizeof(int) * nSampleMax);
523 :
524 : /* -------------------------------------------------------------------- */
525 : /* Transform them to the output coordinate system. */
526 : /* -------------------------------------------------------------------- */
527 321 : CPLTurnFailureIntoWarning(true);
528 321 : pfnTransformer(pTransformArg, FALSE, nSamplePoints, padfX, padfY, padfZ,
529 : pabSuccess);
530 321 : CPLTurnFailureIntoWarning(false);
531 :
532 321 : constexpr int SIGN_FINAL_UNINIT = -2;
533 321 : constexpr int SIGN_FINAL_INVALID = 0;
534 321 : int iSignDiscontinuity = SIGN_FINAL_UNINIT;
535 321 : int nFailedCount = 0;
536 321 : const int iSignArray[2] = {-1, 1};
537 30897 : for (int i = 0; i < nSamplePoints; i++)
538 : {
539 30576 : if (pabSuccess[i])
540 : {
541 : // Fix for https://trac.osgeo.org/gdal/ticket/7243
542 : // where echo "-2050000.000 2050000.000" |
543 : // gdaltransform -s_srs EPSG:3411 -t_srs EPSG:4326
544 : // gives "-180 63.691332898492"
545 : // but we would rather like 180
546 27612 : if (iSignDiscontinuity == 1 || iSignDiscontinuity == -1)
547 : {
548 9777 : if (!((iSignDiscontinuity * padfX[i] > 0 &&
549 9728 : iSignDiscontinuity * padfX[i] <= 180.0) ||
550 50 : (fabs(padfX[i] - iSignDiscontinuity * -180.0) < 1e-8)))
551 : {
552 27 : iSignDiscontinuity = SIGN_FINAL_INVALID;
553 : }
554 : }
555 17835 : else if (iSignDiscontinuity == SIGN_FINAL_UNINIT)
556 : {
557 741 : for (const auto &iSign : iSignArray)
558 : {
559 554 : if ((iSign * padfX[i] > 0 && iSign * padfX[i] <= 180.0) ||
560 421 : (fabs(padfX[i] - iSign * -180.0) < 1e-8))
561 : {
562 133 : iSignDiscontinuity = iSign;
563 133 : break;
564 : }
565 : }
566 320 : if (iSignDiscontinuity == SIGN_FINAL_UNINIT)
567 : {
568 187 : iSignDiscontinuity = SIGN_FINAL_INVALID;
569 : }
570 : }
571 : }
572 : else
573 : {
574 2964 : nFailedCount++;
575 : }
576 : }
577 :
578 321 : if (iSignDiscontinuity == 1 || iSignDiscontinuity == -1)
579 : {
580 10338 : for (int i = 0; i < nSamplePoints; i++)
581 : {
582 10232 : if (pabSuccess[i])
583 : {
584 9615 : if (fabs(padfX[i] - iSignDiscontinuity * -180.0) < 1e-8)
585 : {
586 2 : double axTemp[2] = {iSignDiscontinuity * -180.0,
587 2 : iSignDiscontinuity * 180.0};
588 2 : double ayTemp[2] = {padfY[i], padfY[i]};
589 2 : double azTemp[2] = {padfZ[i], padfZ[i]};
590 2 : int abSuccess[2] = {FALSE, FALSE};
591 2 : CPLTurnFailureIntoWarning(true);
592 2 : if (pfnTransformer(pTransformArg, TRUE, 2, axTemp, ayTemp,
593 2 : azTemp, abSuccess) &&
594 4 : fabs(axTemp[0] - axTemp[1]) < 1e-8 &&
595 2 : fabs(ayTemp[0] - ayTemp[1]) < 1e-8)
596 : {
597 2 : padfX[i] = iSignDiscontinuity * 180.0;
598 : }
599 2 : CPLTurnFailureIntoWarning(false);
600 : }
601 : }
602 : }
603 : }
604 :
605 : /* -------------------------------------------------------------------- */
606 : /* Check if the computed target coordinates are revertable. */
607 : /* If not, try the detailed grid sampling. */
608 : /* -------------------------------------------------------------------- */
609 321 : if (nFailedCount)
610 : {
611 49 : CPLDebug("WARP", "At least one point failed after direct transform");
612 : }
613 : else
614 : {
615 272 : memcpy(padfXRevert, padfX, nSamplePoints * sizeof(double));
616 272 : memcpy(padfYRevert, padfY, nSamplePoints * sizeof(double));
617 272 : memcpy(padfZRevert, padfZ, nSamplePoints * sizeof(double));
618 272 : CPLTurnFailureIntoWarning(true);
619 272 : pfnTransformer(pTransformArg, TRUE, nSamplePoints, padfXRevert,
620 : padfYRevert, padfZRevert, pabSuccess);
621 272 : CPLTurnFailureIntoWarning(false);
622 :
623 24863 : for (int i = 0; nFailedCount == 0 && i < nSamplePoints; i++)
624 : {
625 24607 : if (!pabSuccess[i])
626 : {
627 16 : nFailedCount++;
628 16 : break;
629 : }
630 :
631 24591 : double dfRatio = (i % nStepsPlusOne) * dfStep;
632 24591 : if (dfRatio > 0.99)
633 1022 : dfRatio = 1.0;
634 :
635 24591 : double dfExpectedX = 0.0;
636 24591 : double dfExpectedY = 0.0;
637 24591 : if (i < nStepsPlusOne)
638 : {
639 6279 : dfExpectedX = dfRatio * nInXSize;
640 : }
641 18312 : else if (i < 2 * nStepsPlusOne)
642 : {
643 6232 : dfExpectedX = dfRatio * nInXSize;
644 6232 : dfExpectedY = nInYSize;
645 : }
646 12080 : else if (i < 3 * nStepsPlusOne)
647 : {
648 6102 : dfExpectedY = dfRatio * nInYSize;
649 : }
650 : else
651 : {
652 5978 : dfExpectedX = nInXSize;
653 5978 : dfExpectedY = dfRatio * nInYSize;
654 : }
655 :
656 24591 : if (fabs(padfXRevert[i] - dfExpectedX) >
657 24591 : nInXSize / static_cast<double>(nSteps) ||
658 24582 : fabs(padfYRevert[i] - dfExpectedY) >
659 24582 : nInYSize / static_cast<double>(nSteps))
660 9 : nFailedCount++;
661 : }
662 272 : if (nFailedCount != 0)
663 25 : CPLDebug("WARP",
664 : "At least one point failed after revert transform");
665 : }
666 :
667 : /* -------------------------------------------------------------------- */
668 : /* If any of the edge points failed to transform, we need to */
669 : /* build a fairly detailed internal grid of points instead to */
670 : /* help identify the area that is transformable. */
671 : /* -------------------------------------------------------------------- */
672 321 : if (nFailedCount)
673 : {
674 74 : nSamplePoints = 0;
675 :
676 : // Take N_STEPS steps.
677 1782 : for (int iStep = 0; iStep <= nSteps; iStep++)
678 : {
679 1708 : double dfRatio = (iStep == nSteps) ? 1.0 : iStep * dfStep;
680 :
681 51102 : for (int iStep2 = 0; iStep2 <= nSteps; iStep2++)
682 : {
683 49394 : const double dfRatio2 =
684 49394 : iStep2 == nSteps ? 1.0 : iStep2 * dfStep;
685 :
686 : // From top to bottom, from left to right.
687 49394 : padfX[nSamplePoints] = dfRatio2 * nInXSize;
688 49394 : padfY[nSamplePoints] = dfRatio * nInYSize;
689 49394 : padfZ[nSamplePoints] = 0.0;
690 49394 : nSamplePoints++;
691 : }
692 : }
693 :
694 74 : CPLAssert(nSamplePoints == nSampleMax);
695 :
696 74 : CPLTurnFailureIntoWarning(true);
697 74 : pfnTransformer(pTransformArg, FALSE, nSamplePoints, padfX, padfY, padfZ,
698 : pabSuccess);
699 74 : CPLTurnFailureIntoWarning(false);
700 : }
701 :
702 : /* -------------------------------------------------------------------- */
703 : /* Collect the bounds, ignoring any failed points. */
704 : /* -------------------------------------------------------------------- */
705 321 : double dfMinXOut = 0.0;
706 321 : double dfMinYOut = 0.0;
707 321 : double dfMaxXOut = 0.0;
708 321 : double dfMaxYOut = 0.0;
709 321 : bool bGotInitialPoint = false;
710 :
711 321 : nFailedCount = 0;
712 73459 : for (int i = 0; i < nSamplePoints; i++)
713 : {
714 73138 : int x_i = 0;
715 73138 : int y_i = 0;
716 :
717 73138 : if (nSamplePoints == nSampleMax)
718 : {
719 49394 : x_i = i % nStepsPlusOne;
720 49394 : y_i = i / nStepsPlusOne;
721 : }
722 : else
723 : {
724 23744 : if (i < 2 * nStepsPlusOne)
725 : {
726 11872 : x_i = i % nStepsPlusOne;
727 11872 : y_i = (i < nStepsPlusOne) ? 0 : nSteps;
728 : }
729 : }
730 :
731 73138 : if (x_i > 0 && (pabSuccess[i - 1] || pabSuccess[i]))
732 : {
733 47905 : double x_out_before = padfX[i - 1];
734 47905 : double x_out_after = padfX[i];
735 47905 : int nIter = 0;
736 47905 : double x_in_before =
737 47905 : static_cast<double>(x_i - 1) * nInXSize / nSteps;
738 47905 : double x_in_after = static_cast<double>(x_i) * nInXSize / nSteps;
739 47905 : int invalid_before = !(pabSuccess[i - 1]);
740 47905 : int invalid_after = !(pabSuccess[i]);
741 :
742 : // Detect discontinuity in target coordinates when the target x
743 : // coordinates change sign. This may be a false positive when the
744 : // target tx is around 0 Dichotomic search to reduce the interval
745 : // to near the discontinuity and get a better out extent.
746 63581 : while ((invalid_before || invalid_after ||
747 133758 : x_out_before * x_out_after < 0.0) &&
748 : nIter < 16)
749 : {
750 22272 : double x = (x_in_before + x_in_after) / 2.0;
751 22272 : double y = static_cast<double>(y_i) * nInYSize / nSteps;
752 22272 : double z = 0.0;
753 22272 : int bSuccess = TRUE;
754 22272 : if (pfnTransformer(pTransformArg, FALSE, 1, &x, &y, &z,
755 39062 : &bSuccess) &&
756 16790 : bSuccess)
757 : {
758 16790 : if (bGotInitialPoint)
759 : {
760 16772 : dfMinXOut = std::min(dfMinXOut, x);
761 16772 : dfMinYOut = std::min(dfMinYOut, y);
762 16772 : dfMaxXOut = std::max(dfMaxXOut, x);
763 16772 : dfMaxYOut = std::max(dfMaxYOut, y);
764 : }
765 : else
766 : {
767 18 : bGotInitialPoint = true;
768 18 : dfMinXOut = x;
769 18 : dfMaxXOut = x;
770 18 : dfMinYOut = y;
771 18 : dfMaxYOut = y;
772 : }
773 :
774 16790 : if (invalid_before || x_out_before * x < 0)
775 : {
776 9033 : invalid_after = FALSE;
777 9033 : x_in_after = (x_in_before + x_in_after) / 2.0;
778 9033 : x_out_after = x;
779 : }
780 : else
781 : {
782 7757 : invalid_before = FALSE;
783 7757 : x_out_before = x;
784 7757 : x_in_before = (x_in_before + x_in_after) / 2.0;
785 : }
786 : }
787 : else
788 : {
789 5482 : if (invalid_before)
790 : {
791 2736 : x_in_before = (x_in_before + x_in_after) / 2.0;
792 : }
793 2746 : else if (invalid_after)
794 : {
795 2746 : x_in_after = (x_in_before + x_in_after) / 2.0;
796 : }
797 : else
798 : {
799 0 : break;
800 : }
801 : }
802 22272 : nIter++;
803 : }
804 : }
805 :
806 73138 : if (!pabSuccess[i])
807 : {
808 12036 : nFailedCount++;
809 12036 : continue;
810 : }
811 :
812 61102 : if (bGotInitialPoint)
813 : {
814 60800 : dfMinXOut = std::min(dfMinXOut, padfX[i]);
815 60800 : dfMinYOut = std::min(dfMinYOut, padfY[i]);
816 60800 : dfMaxXOut = std::max(dfMaxXOut, padfX[i]);
817 60800 : dfMaxYOut = std::max(dfMaxYOut, padfY[i]);
818 : }
819 : else
820 : {
821 302 : bGotInitialPoint = true;
822 302 : dfMinXOut = padfX[i];
823 302 : dfMaxXOut = padfX[i];
824 302 : dfMinYOut = padfY[i];
825 302 : dfMaxYOut = padfY[i];
826 : }
827 : }
828 :
829 321 : if (nFailedCount > nSamplePoints - 10)
830 : {
831 4 : CPLError(CE_Failure, CPLE_AppDefined,
832 : "Too many points (%d out of %d) failed to transform, "
833 : "unable to compute output bounds.",
834 : nFailedCount, nSamplePoints);
835 :
836 4 : CPLFree(padfX);
837 4 : CPLFree(padfXRevert);
838 4 : CPLFree(pabSuccess);
839 :
840 4 : return CE_Failure;
841 : }
842 :
843 317 : if (nFailedCount)
844 45 : CPLDebug("GDAL",
845 : "GDALSuggestedWarpOutput(): %d out of %d points failed to "
846 : "transform.",
847 : nFailedCount, nSamplePoints);
848 :
849 317 : bool bIsGeographicCoordsDeg = false;
850 317 : if (bIsGDALGenImgProjTransform)
851 : {
852 317 : const GDALGenImgProjTransformInfo *pGIPTI =
853 : static_cast<const GDALGenImgProjTransformInfo *>(pTransformArg);
854 317 : if (pGIPTI->pSrcTransformer == GDALGeoLocTransform &&
855 28 : pGIPTI->pDstTransformer == nullptr &&
856 28 : pGIPTI->adfDstGeoTransform[0] == 0 &&
857 26 : pGIPTI->adfDstGeoTransform[1] == 1 &&
858 26 : pGIPTI->adfDstGeoTransform[2] == 0 &&
859 26 : pGIPTI->adfDstGeoTransform[3] == 0 &&
860 26 : pGIPTI->adfDstGeoTransform[4] == 0 &&
861 26 : pGIPTI->adfDstGeoTransform[5] == 1)
862 : {
863 : /* --------------------------------------------------------------------
864 : */
865 : /* Special case for geolocation array, to quickly find the
866 : * bounds. */
867 : /* --------------------------------------------------------------------
868 : */
869 26 : const GDALGeoLocTransformInfo *pGLTI =
870 : static_cast<const GDALGeoLocTransformInfo *>(
871 : pGIPTI->pSrcTransformArg);
872 :
873 26 : if (pGIPTI->pReproject == nullptr)
874 : {
875 : const char *pszGLSRS =
876 26 : CSLFetchNameValue(pGLTI->papszGeolocationInfo, "SRS");
877 26 : if (pszGLSRS == nullptr)
878 : {
879 4 : bIsGeographicCoordsDeg = true;
880 : }
881 : else
882 : {
883 44 : OGRSpatialReference oSRS;
884 22 : if (oSRS.SetFromUserInput(pszGLSRS) == OGRERR_NONE &&
885 42 : oSRS.IsGeographic() &&
886 20 : std::fabs(oSRS.GetAngularUnits() -
887 20 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-9)
888 : {
889 20 : bIsGeographicCoordsDeg = true;
890 : }
891 : }
892 : }
893 :
894 208 : for (const auto &xy :
895 26 : {std::pair<double, double>(pGLTI->dfMinX, pGLTI->dfYAtMinX),
896 26 : std::pair<double, double>(pGLTI->dfXAtMinY, pGLTI->dfMinY),
897 26 : std::pair<double, double>(pGLTI->dfMaxX, pGLTI->dfYAtMaxX),
898 130 : std::pair<double, double>(pGLTI->dfXAtMaxY, pGLTI->dfMaxY)})
899 : {
900 104 : double x = xy.first;
901 104 : double y = xy.second;
902 104 : if (pGLTI->bSwapXY)
903 : {
904 4 : std::swap(x, y);
905 : }
906 104 : double xOut = std::numeric_limits<double>::quiet_NaN();
907 104 : double yOut = std::numeric_limits<double>::quiet_NaN();
908 104 : if (pGIPTI->pReproject == nullptr ||
909 0 : pGIPTI->pReproject(pGIPTI->pReprojectArg, false, 1, &x, &y,
910 : nullptr, nullptr))
911 : {
912 104 : xOut = x;
913 104 : yOut = y;
914 : }
915 104 : dfMinXOut = std::min(dfMinXOut, xOut);
916 104 : dfMinYOut = std::min(dfMinYOut, yOut);
917 104 : dfMaxXOut = std::max(dfMaxXOut, xOut);
918 104 : dfMaxYOut = std::max(dfMaxYOut, yOut);
919 26 : }
920 : }
921 291 : else if (pGIPTI->pSrcTransformer == nullptr &&
922 248 : pGIPTI->pDstTransformer == nullptr &&
923 248 : pGIPTI->pReproject == GDALReprojectionTransform &&
924 237 : pGIPTI->adfDstGeoTransform[0] == 0 &&
925 235 : pGIPTI->adfDstGeoTransform[1] == 1 &&
926 235 : pGIPTI->adfDstGeoTransform[2] == 0 &&
927 235 : pGIPTI->adfDstGeoTransform[3] == 0 &&
928 235 : pGIPTI->adfDstGeoTransform[4] == 0 &&
929 235 : pGIPTI->adfDstGeoTransform[5] == 1)
930 : {
931 : /* ------------------------------------------------------------- */
932 : /* Special case for warping using source geotransform and */
933 : /* reprojection to deal with the poles. */
934 : /* ------------------------------------------------------------- */
935 235 : const GDALReprojectionTransformInfo *psRTI =
936 : static_cast<const GDALReprojectionTransformInfo *>(
937 : pGIPTI->pReprojectArg);
938 : const OGRSpatialReference *poSourceCRS =
939 235 : psRTI->poForwardTransform->GetSourceCS();
940 : const OGRSpatialReference *poTargetCRS =
941 235 : psRTI->poForwardTransform->GetTargetCS();
942 469 : if (poTargetCRS != nullptr &&
943 234 : psRTI->poReverseTransform != nullptr &&
944 234 : poTargetCRS->IsGeographic() &&
945 90 : fabs(poTargetCRS->GetAngularUnits() -
946 559 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-9 &&
947 90 : (!poSourceCRS || !poSourceCRS->IsGeographic()))
948 : {
949 83 : bIsGeographicCoordsDeg = true;
950 :
951 83 : std::unique_ptr<CPLConfigOptionSetter> poSetter;
952 83 : if (pGIPTI->bCheckWithInvertPROJ)
953 : {
954 : // CHECK_WITH_INVERT_PROJ=YES prevent reliable
955 : // transformation of poles.
956 4 : poSetter = std::make_unique<CPLConfigOptionSetter>(
957 4 : "CHECK_WITH_INVERT_PROJ", "NO", false);
958 4 : GDALRefreshGenImgProjTransformer(pTransformArg);
959 : // GDALRefreshGenImgProjTransformer() has invalidated psRTI
960 4 : psRTI = static_cast<const GDALReprojectionTransformInfo *>(
961 : pGIPTI->pReprojectArg);
962 : }
963 :
964 249 : for (const auto &sign : iSignArray)
965 : {
966 166 : double X = 0.0;
967 166 : const double Yinit = 90.0 * sign;
968 166 : double Y = Yinit;
969 166 : if (psRTI->poReverseTransform->Transform(1, &X, &Y))
970 : {
971 112 : const auto invGT = pGIPTI->adfSrcInvGeoTransform;
972 112 : const double x = invGT[0] + X * invGT[1] + Y * invGT[2];
973 112 : const double y = invGT[3] + X * invGT[4] + Y * invGT[5];
974 112 : constexpr double EPSILON = 1e-5;
975 112 : if (x >= -EPSILON && x <= nInXSize + EPSILON &&
976 26 : y >= -EPSILON && y <= nInYSize + EPSILON)
977 : {
978 6 : if (psRTI->poForwardTransform->Transform(1, &X,
979 12 : &Y) &&
980 6 : fabs(Y - Yinit) <= 1e-6)
981 : {
982 6 : bool bMinXMaxXSet = false;
983 6 : if (poSourceCRS)
984 : {
985 : const char *pszProjection =
986 6 : poSourceCRS->GetAttrValue("PROJECTION");
987 6 : if (pszProjection &&
988 6 : EQUAL(pszProjection,
989 : SRS_PT_ORTHOGRAPHIC))
990 : {
991 : const double dfLon0 =
992 4 : poSourceCRS->GetNormProjParm(
993 : SRS_PP_CENTRAL_MERIDIAN, 0.0);
994 4 : dfMinXOut = dfLon0 - 90;
995 4 : dfMaxXOut = dfLon0 + 90;
996 4 : bMinXMaxXSet = true;
997 : }
998 : }
999 6 : if (!bMinXMaxXSet)
1000 : {
1001 2 : dfMinXOut = -180;
1002 2 : dfMaxXOut = 180;
1003 : }
1004 6 : if (sign < 0)
1005 2 : dfMinYOut = Yinit;
1006 : else
1007 4 : dfMaxYOut = Yinit;
1008 : }
1009 : }
1010 : }
1011 : }
1012 :
1013 83 : if (poSetter)
1014 : {
1015 4 : poSetter.reset();
1016 4 : GDALRefreshGenImgProjTransformer(pTransformArg);
1017 4 : pGIPTI = static_cast<const GDALGenImgProjTransformInfo *>(
1018 : pTransformArg);
1019 4 : psRTI = static_cast<const GDALReprojectionTransformInfo *>(
1020 : pGIPTI->pReprojectArg);
1021 4 : poSourceCRS = psRTI->poForwardTransform->GetSourceCS();
1022 4 : poTargetCRS = psRTI->poForwardTransform->GetTargetCS();
1023 : }
1024 : }
1025 :
1026 : // Use TransformBounds() to handle more particular cases
1027 235 : if (poSourceCRS != nullptr && poTargetCRS != nullptr &&
1028 234 : pGIPTI->adfSrcGeoTransform[1] != 0 &&
1029 234 : pGIPTI->adfSrcGeoTransform[2] == 0 &&
1030 234 : pGIPTI->adfSrcGeoTransform[4] == 0 &&
1031 234 : pGIPTI->adfSrcGeoTransform[5] != 0)
1032 : {
1033 234 : const double dfULX = pGIPTI->adfSrcGeoTransform[0];
1034 234 : const double dfULY = pGIPTI->adfSrcGeoTransform[3];
1035 234 : const double dfLRX =
1036 234 : dfULX + pGIPTI->adfSrcGeoTransform[1] * nInXSize;
1037 234 : const double dfLRY =
1038 234 : dfULY + pGIPTI->adfSrcGeoTransform[5] * nInYSize;
1039 234 : const double dfMinSrcX = std::min(dfULX, dfLRX);
1040 234 : const double dfMinSrcY = std::min(dfULY, dfLRY);
1041 234 : const double dfMaxSrcX = std::max(dfULX, dfLRX);
1042 234 : const double dfMaxSrcY = std::max(dfULY, dfLRY);
1043 234 : double dfTmpMinXOut = std::numeric_limits<double>::max();
1044 234 : double dfTmpMinYOut = std::numeric_limits<double>::max();
1045 234 : double dfTmpMaxXOut = std::numeric_limits<double>::min();
1046 234 : double dfTmpMaxYOut = std::numeric_limits<double>::min();
1047 468 : if (psRTI->poForwardTransform->TransformBounds(
1048 : dfMinSrcX, dfMinSrcY, dfMaxSrcX, dfMaxSrcY,
1049 : &dfTmpMinXOut, &dfTmpMinYOut, &dfTmpMaxXOut,
1050 : &dfTmpMaxYOut,
1051 234 : 2)) // minimum number of points as we already have a
1052 : // logic above to sample
1053 : {
1054 228 : dfMinXOut = std::min(dfMinXOut, dfTmpMinXOut);
1055 228 : dfMinYOut = std::min(dfMinYOut, dfTmpMinYOut);
1056 228 : dfMaxXOut = std::max(dfMaxXOut, dfTmpMaxXOut);
1057 228 : dfMaxYOut = std::max(dfMaxYOut, dfTmpMaxYOut);
1058 : }
1059 : }
1060 : }
1061 : }
1062 :
1063 : /* -------------------------------------------------------------------- */
1064 : /* Compute the distance in "georeferenced" units from the top */
1065 : /* corner of the transformed input image to the bottom left */
1066 : /* corner of the transformed input. Use this distance to */
1067 : /* compute an approximate pixel size in the output */
1068 : /* georeferenced coordinates. */
1069 : /* -------------------------------------------------------------------- */
1070 317 : double dfDiagonalDist = 0.0;
1071 317 : double dfDeltaX = 0.0;
1072 317 : double dfDeltaY = 0.0;
1073 :
1074 317 : if (pabSuccess[0] && pabSuccess[nSamplePoints - 1])
1075 : {
1076 275 : dfDeltaX = padfX[nSamplePoints - 1] - padfX[0];
1077 275 : dfDeltaY = padfY[nSamplePoints - 1] - padfY[0];
1078 : // In some cases this can result in 0 values. See #5980
1079 : // Fallback to safer method in that case.
1080 : }
1081 317 : if (dfDeltaX == 0.0 || dfDeltaY == 0.0)
1082 : {
1083 47 : dfDeltaX = dfMaxXOut - dfMinXOut;
1084 47 : dfDeltaY = dfMaxYOut - dfMinYOut;
1085 : }
1086 :
1087 317 : dfDiagonalDist = sqrt(dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY);
1088 :
1089 : /* -------------------------------------------------------------------- */
1090 : /* Compute a pixel size from this. */
1091 : /* -------------------------------------------------------------------- */
1092 : const double dfPixelSize =
1093 317 : dfDiagonalDist / sqrt(static_cast<double>(nInXSize) * nInXSize +
1094 317 : static_cast<double>(nInYSize) * nInYSize);
1095 :
1096 317 : const double dfPixels = (dfMaxXOut - dfMinXOut) / dfPixelSize;
1097 317 : const double dfLines = (dfMaxYOut - dfMinYOut) / dfPixelSize;
1098 :
1099 317 : const int knIntMaxMinusOne = std::numeric_limits<int>::max() - 1;
1100 317 : if (dfPixels > knIntMaxMinusOne || dfLines > knIntMaxMinusOne)
1101 : {
1102 0 : CPLError(CE_Failure, CPLE_AppDefined,
1103 : "Computed dimensions are too big : %.0f x %.0f",
1104 : dfPixels + 0.5, dfLines + 0.5);
1105 :
1106 0 : CPLFree(padfX);
1107 0 : CPLFree(padfXRevert);
1108 0 : CPLFree(pabSuccess);
1109 :
1110 0 : return CE_Failure;
1111 : }
1112 :
1113 317 : if ((nOptions & GDAL_SWO_ROUND_UP_SIZE) != 0)
1114 : {
1115 8 : constexpr double EPS = 1e-5;
1116 8 : *pnPixels = static_cast<int>(std::ceil(dfPixels - EPS));
1117 8 : *pnLines = static_cast<int>(std::ceil(dfLines - EPS));
1118 : }
1119 : else
1120 : {
1121 309 : *pnPixels = static_cast<int>(dfPixels + 0.5);
1122 309 : *pnLines = static_cast<int>(dfLines + 0.5);
1123 : }
1124 :
1125 317 : double dfPixelSizeX = dfPixelSize;
1126 317 : double dfPixelSizeY = dfPixelSize;
1127 :
1128 317 : const double adfRatioArray[] = {0.000, 0.001, 0.010, 0.100, 1.000};
1129 :
1130 : /* -------------------------------------------------------------------- */
1131 : /* Check that the right border is not completely out of source */
1132 : /* image. If so, adjust the x pixel size a bit in the hope it will */
1133 : /* fit. */
1134 : /* -------------------------------------------------------------------- */
1135 496 : for (const auto &dfRatio : adfRatioArray)
1136 : {
1137 488 : const double dfTryPixelSizeX =
1138 488 : dfPixelSizeX - dfPixelSizeX * dfRatio / *pnPixels;
1139 488 : double adfExtent[4] = {dfMinXOut, dfMaxYOut - (*pnLines) * dfPixelSizeY,
1140 488 : dfMinXOut + (*pnPixels) * dfTryPixelSizeX,
1141 488 : dfMaxYOut};
1142 488 : if (!GDALSuggestedWarpOutput2_MustAdjustForRightBorder(
1143 : pfnTransformer, pTransformArg, adfExtent, *pnPixels, *pnLines,
1144 : dfTryPixelSizeX, dfPixelSizeY))
1145 : {
1146 309 : dfPixelSizeX = dfTryPixelSizeX;
1147 309 : break;
1148 : }
1149 : }
1150 :
1151 : /* -------------------------------------------------------------------- */
1152 : /* Check that the bottom border is not completely out of source */
1153 : /* image. If so, adjust the y pixel size a bit in the hope it will */
1154 : /* fit. */
1155 : /* -------------------------------------------------------------------- */
1156 397 : for (const auto &dfRatio : adfRatioArray)
1157 : {
1158 387 : const double dfTryPixelSizeY =
1159 387 : dfPixelSizeY - dfPixelSizeY * dfRatio / *pnLines;
1160 : double adfExtent[4] = {
1161 387 : dfMinXOut, dfMaxYOut - (*pnLines) * dfTryPixelSizeY,
1162 387 : dfMinXOut + (*pnPixels) * dfPixelSizeX, dfMaxYOut};
1163 387 : if (!GDALSuggestedWarpOutput2_MustAdjustForBottomBorder(
1164 : pfnTransformer, pTransformArg, adfExtent, *pnPixels, *pnLines,
1165 : dfPixelSizeX, dfTryPixelSizeY))
1166 : {
1167 307 : dfPixelSizeY = dfTryPixelSizeY;
1168 307 : break;
1169 : }
1170 : }
1171 :
1172 : /* -------------------------------------------------------------------- */
1173 : /* Recompute some bounds so that all return values are consistent */
1174 : /* -------------------------------------------------------------------- */
1175 317 : double dfMaxXOutNew = dfMinXOut + (*pnPixels) * dfPixelSizeX;
1176 317 : if (bIsGeographicCoordsDeg &&
1177 107 : ((dfMaxXOut <= 180 && dfMaxXOutNew > 180) || dfMaxXOut == 180))
1178 : {
1179 3 : dfMaxXOut = 180;
1180 3 : dfPixelSizeX = (dfMaxXOut - dfMinXOut) / *pnPixels;
1181 : }
1182 : else
1183 : {
1184 314 : dfMaxXOut = dfMaxXOutNew;
1185 : }
1186 :
1187 317 : double dfMinYOutNew = dfMaxYOut - (*pnLines) * dfPixelSizeY;
1188 317 : if (bIsGeographicCoordsDeg && dfMinYOut >= -90 && dfMinYOutNew < -90)
1189 : {
1190 0 : dfMinYOut = -90;
1191 0 : dfPixelSizeY = (dfMaxYOut - dfMinYOut) / *pnLines;
1192 : }
1193 : else
1194 : {
1195 317 : dfMinYOut = dfMinYOutNew;
1196 : }
1197 :
1198 : /* -------------------------------------------------------------------- */
1199 : /* Return raw extents. */
1200 : /* -------------------------------------------------------------------- */
1201 317 : padfExtent[0] = dfMinXOut;
1202 317 : padfExtent[1] = dfMinYOut;
1203 317 : padfExtent[2] = dfMaxXOut;
1204 317 : padfExtent[3] = dfMaxYOut;
1205 :
1206 : /* -------------------------------------------------------------------- */
1207 : /* Set the output geotransform. */
1208 : /* -------------------------------------------------------------------- */
1209 317 : padfGeoTransformOut[0] = dfMinXOut;
1210 317 : padfGeoTransformOut[1] = dfPixelSizeX;
1211 317 : padfGeoTransformOut[2] = 0.0;
1212 317 : padfGeoTransformOut[3] = dfMaxYOut;
1213 317 : padfGeoTransformOut[4] = 0.0;
1214 317 : padfGeoTransformOut[5] = -dfPixelSizeY;
1215 :
1216 317 : CPLFree(padfX);
1217 317 : CPLFree(padfXRevert);
1218 317 : CPLFree(pabSuccess);
1219 :
1220 317 : return CE_None;
1221 : }
1222 :
1223 : /************************************************************************/
1224 : /* GetCurrentCheckWithInvertPROJ() */
1225 : /************************************************************************/
1226 :
1227 2654 : static bool GetCurrentCheckWithInvertPROJ()
1228 : {
1229 2654 : return CPLTestBool(CPLGetConfigOption("CHECK_WITH_INVERT_PROJ", "NO"));
1230 : }
1231 :
1232 : /************************************************************************/
1233 : /* GDALCreateGenImgProjTransformerInternal() */
1234 : /************************************************************************/
1235 :
1236 : static void *GDALCreateSimilarGenImgProjTransformer(void *hTransformArg,
1237 : double dfRatioX,
1238 : double dfRatioY);
1239 :
1240 1659 : static GDALGenImgProjTransformInfo *GDALCreateGenImgProjTransformerInternal()
1241 : {
1242 : /* -------------------------------------------------------------------- */
1243 : /* Initialize the transform info. */
1244 : /* -------------------------------------------------------------------- */
1245 : GDALGenImgProjTransformInfo *psInfo =
1246 : static_cast<GDALGenImgProjTransformInfo *>(
1247 1659 : CPLCalloc(sizeof(GDALGenImgProjTransformInfo), 1));
1248 :
1249 1659 : memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
1250 : strlen(GDAL_GTI2_SIGNATURE));
1251 1659 : psInfo->sTI.pszClassName = GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME;
1252 1659 : psInfo->sTI.pfnTransform = GDALGenImgProjTransform;
1253 1659 : psInfo->sTI.pfnCleanup = GDALDestroyGenImgProjTransformer;
1254 1659 : psInfo->sTI.pfnSerialize = GDALSerializeGenImgProjTransformer;
1255 1659 : psInfo->sTI.pfnCreateSimilar = GDALCreateSimilarGenImgProjTransformer;
1256 :
1257 1659 : psInfo->bCheckWithInvertPROJ = GetCurrentCheckWithInvertPROJ();
1258 1659 : psInfo->bHasCustomTransformationPipeline = false;
1259 :
1260 1659 : return psInfo;
1261 : }
1262 :
1263 : /************************************************************************/
1264 : /* GDALCreateSimilarGenImgProjTransformer() */
1265 : /************************************************************************/
1266 :
1267 32 : static void *GDALCreateSimilarGenImgProjTransformer(void *hTransformArg,
1268 : double dfRatioX,
1269 : double dfRatioY)
1270 : {
1271 32 : VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarGenImgProjTransformer",
1272 : nullptr);
1273 :
1274 32 : GDALGenImgProjTransformInfo *psInfo =
1275 : static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
1276 :
1277 : GDALGenImgProjTransformInfo *psClonedInfo =
1278 32 : GDALCreateGenImgProjTransformerInternal();
1279 :
1280 32 : memcpy(psClonedInfo, psInfo, sizeof(GDALGenImgProjTransformInfo));
1281 :
1282 32 : psClonedInfo->bCheckWithInvertPROJ = GetCurrentCheckWithInvertPROJ();
1283 :
1284 32 : if (psClonedInfo->pSrcTransformArg)
1285 7 : psClonedInfo->pSrcTransformArg = GDALCreateSimilarTransformer(
1286 : psInfo->pSrcTransformArg, dfRatioX, dfRatioY);
1287 25 : else if (dfRatioX != 1.0 || dfRatioY != 1.0)
1288 : {
1289 10 : if (psClonedInfo->adfSrcGeoTransform[2] == 0.0 &&
1290 10 : psClonedInfo->adfSrcGeoTransform[4] == 0.0)
1291 : {
1292 10 : psClonedInfo->adfSrcGeoTransform[1] *= dfRatioX;
1293 10 : psClonedInfo->adfSrcGeoTransform[5] *= dfRatioY;
1294 : }
1295 : else
1296 : {
1297 : // If the x and y ratios are not equal, then we cannot really
1298 : // compute a geotransform.
1299 0 : psClonedInfo->adfSrcGeoTransform[1] *= dfRatioX;
1300 0 : psClonedInfo->adfSrcGeoTransform[2] *= dfRatioX;
1301 0 : psClonedInfo->adfSrcGeoTransform[4] *= dfRatioX;
1302 0 : psClonedInfo->adfSrcGeoTransform[5] *= dfRatioX;
1303 : }
1304 10 : if (!GDALInvGeoTransform(psClonedInfo->adfSrcGeoTransform,
1305 10 : psClonedInfo->adfSrcInvGeoTransform))
1306 : {
1307 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
1308 0 : GDALDestroyGenImgProjTransformer(psClonedInfo);
1309 0 : return nullptr;
1310 : }
1311 : }
1312 :
1313 32 : if (psClonedInfo->pReprojectArg)
1314 11 : psClonedInfo->pReprojectArg =
1315 11 : GDALCloneTransformer(psInfo->pReprojectArg);
1316 :
1317 32 : if (psClonedInfo->pDstTransformArg)
1318 0 : psClonedInfo->pDstTransformArg =
1319 0 : GDALCloneTransformer(psInfo->pDstTransformArg);
1320 :
1321 32 : return psClonedInfo;
1322 : }
1323 :
1324 : /************************************************************************/
1325 : /* GDALCreateGenImgProjTransformer() */
1326 : /************************************************************************/
1327 :
1328 : /**
1329 : * Create image to image transformer.
1330 : *
1331 : * This function creates a transformation object that maps from pixel/line
1332 : * coordinates on one image to pixel/line coordinates on another image. The
1333 : * images may potentially be georeferenced in different coordinate systems,
1334 : * and may used GCPs to map between their pixel/line coordinates and
1335 : * georeferenced coordinates (as opposed to the default assumption that their
1336 : * geotransform should be used).
1337 : *
1338 : * This transformer potentially performs three concatenated transformations.
1339 : *
1340 : * The first stage is from source image pixel/line coordinates to source
1341 : * image georeferenced coordinates, and may be done using the geotransform,
1342 : * or if not defined using a polynomial model derived from GCPs. If GCPs
1343 : * are used this stage is accomplished using GDALGCPTransform().
1344 : *
1345 : * The second stage is to change projections from the source coordinate system
1346 : * to the destination coordinate system, assuming they differ. This is
1347 : * accomplished internally using GDALReprojectionTransform().
1348 : *
1349 : * The third stage is converting from destination image georeferenced
1350 : * coordinates to destination image coordinates. This is done using the
1351 : * destination image geotransform, or if not available, using a polynomial
1352 : * model derived from GCPs. If GCPs are used this stage is accomplished using
1353 : * GDALGCPTransform(). This stage is skipped if hDstDS is NULL when the
1354 : * transformation is created.
1355 : *
1356 : * @param hSrcDS source dataset, or NULL.
1357 : * @param pszSrcWKT the coordinate system for the source dataset. If NULL,
1358 : * it will be read from the dataset itself.
1359 : * @param hDstDS destination dataset (or NULL).
1360 : * @param pszDstWKT the coordinate system for the destination dataset. If
1361 : * NULL, and hDstDS not NULL, it will be read from the destination dataset.
1362 : * @param bGCPUseOK TRUE if GCPs should be used if the geotransform is not
1363 : * available on the source dataset (not destination).
1364 : * @param dfGCPErrorThreshold ignored/deprecated.
1365 : * @param nOrder the maximum order to use for GCP derived polynomials if
1366 : * possible. Use 0 to autoselect, or -1 for thin plate splines.
1367 : *
1368 : * @return handle suitable for use GDALGenImgProjTransform(), and to be
1369 : * deallocated with GDALDestroyGenImgProjTransformer().
1370 : */
1371 :
1372 90 : void *GDALCreateGenImgProjTransformer(GDALDatasetH hSrcDS,
1373 : const char *pszSrcWKT,
1374 : GDALDatasetH hDstDS,
1375 : const char *pszDstWKT, int bGCPUseOK,
1376 : CPL_UNUSED double dfGCPErrorThreshold,
1377 : int nOrder)
1378 : {
1379 90 : char **papszOptions = nullptr;
1380 :
1381 90 : if (pszSrcWKT != nullptr)
1382 31 : papszOptions = CSLSetNameValue(papszOptions, "SRC_SRS", pszSrcWKT);
1383 90 : if (pszDstWKT != nullptr)
1384 31 : papszOptions = CSLSetNameValue(papszOptions, "DST_SRS", pszDstWKT);
1385 90 : if (!bGCPUseOK)
1386 0 : papszOptions = CSLSetNameValue(papszOptions, "GCPS_OK", "FALSE");
1387 90 : if (nOrder != 0)
1388 0 : papszOptions = CSLSetNameValue(papszOptions, "MAX_GCP_ORDER",
1389 0 : CPLString().Printf("%d", nOrder));
1390 :
1391 90 : void *pRet = GDALCreateGenImgProjTransformer2(hSrcDS, hDstDS, papszOptions);
1392 90 : CSLDestroy(papszOptions);
1393 :
1394 90 : return pRet;
1395 : }
1396 :
1397 : /************************************************************************/
1398 : /* InsertCenterLong() */
1399 : /* */
1400 : /* Insert a CENTER_LONG Extension entry on a GEOGCS to indicate */
1401 : /* the center longitude of the dataset for wrapping purposes. */
1402 : /************************************************************************/
1403 :
1404 715 : static void InsertCenterLong(GDALDatasetH hDS, OGRSpatialReference *poSRS,
1405 : CPLStringList &aosOptions)
1406 :
1407 : {
1408 1252 : if (!poSRS->IsGeographic() || std::fabs(poSRS->GetAngularUnits() -
1409 537 : CPLAtof(SRS_UA_DEGREE_CONV)) > 1e-9)
1410 : {
1411 179 : return;
1412 : }
1413 :
1414 536 : if (poSRS->GetExtension(nullptr, "CENTER_LONG"))
1415 0 : return;
1416 :
1417 : /* -------------------------------------------------------------------- */
1418 : /* For now we only do this if we have a geotransform since */
1419 : /* other forms require a bunch of extra work. */
1420 : /* -------------------------------------------------------------------- */
1421 536 : double adfGeoTransform[6] = {};
1422 :
1423 536 : if (GDALGetGeoTransform(hDS, adfGeoTransform) != CE_None)
1424 0 : return;
1425 :
1426 : /* -------------------------------------------------------------------- */
1427 : /* Compute min/max longitude based on testing the four corners. */
1428 : /* -------------------------------------------------------------------- */
1429 536 : const int nXSize = GDALGetRasterXSize(hDS);
1430 536 : const int nYSize = GDALGetRasterYSize(hDS);
1431 :
1432 : const double dfMinLong =
1433 1072 : std::min(std::min(adfGeoTransform[0] + 0 * adfGeoTransform[1] +
1434 536 : 0 * adfGeoTransform[2],
1435 1072 : adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
1436 536 : 0 * adfGeoTransform[2]),
1437 1072 : std::min(adfGeoTransform[0] + 0 * adfGeoTransform[1] +
1438 536 : nYSize * adfGeoTransform[2],
1439 1072 : adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
1440 536 : nYSize * adfGeoTransform[2]));
1441 : const double dfMaxLong =
1442 1072 : std::max(std::max(adfGeoTransform[0] + 0 * adfGeoTransform[1] +
1443 536 : 0 * adfGeoTransform[2],
1444 1072 : adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
1445 536 : 0 * adfGeoTransform[2]),
1446 1072 : std::max(adfGeoTransform[0] + 0 * adfGeoTransform[1] +
1447 536 : nYSize * adfGeoTransform[2],
1448 1072 : adfGeoTransform[0] + nXSize * adfGeoTransform[1] +
1449 536 : nYSize * adfGeoTransform[2]));
1450 :
1451 : const double dfEpsilon =
1452 536 : std::max(std::fabs(adfGeoTransform[1]), std::fabs(adfGeoTransform[2]));
1453 : // If the raster covers more than 360 degree (allow an extra pixel),
1454 : // give up
1455 536 : constexpr double RELATIVE_EPSILON = 0.05; // for numeric precision issues
1456 536 : if (dfMaxLong - dfMinLong > 360.0 + dfEpsilon * (1 + RELATIVE_EPSILON))
1457 0 : return;
1458 :
1459 : /* -------------------------------------------------------------------- */
1460 : /* Insert center long. */
1461 : /* -------------------------------------------------------------------- */
1462 536 : const double dfCenterLong = (dfMaxLong + dfMinLong) / 2.0;
1463 536 : aosOptions.SetNameValue("CENTER_LONG", CPLSPrintf("%g", dfCenterLong));
1464 : }
1465 :
1466 : /************************************************************************/
1467 : /* GDALComputeAreaOfInterest() */
1468 : /************************************************************************/
1469 :
1470 1006 : bool GDALComputeAreaOfInterest(OGRSpatialReference *poSRS, double adfGT[6],
1471 : int nXSize, int nYSize,
1472 : double &dfWestLongitudeDeg,
1473 : double &dfSouthLatitudeDeg,
1474 : double &dfEastLongitudeDeg,
1475 : double &dfNorthLatitudeDeg)
1476 : {
1477 1006 : bool ret = false;
1478 :
1479 1006 : if (!poSRS)
1480 0 : return false;
1481 :
1482 1006 : OGRSpatialReference oSrcSRSHoriz(*poSRS);
1483 1006 : if (oSrcSRSHoriz.IsCompound())
1484 : {
1485 17 : oSrcSRSHoriz.StripVertical();
1486 : }
1487 :
1488 1006 : OGRSpatialReference *poGeog = oSrcSRSHoriz.CloneGeogCS();
1489 1006 : if (poGeog)
1490 : {
1491 1006 : poGeog->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1492 1006 : poGeog->SetAngularUnits(SRS_UA_DEGREE, CPLAtof(SRS_UA_DEGREE_CONV));
1493 :
1494 1006 : auto poCT = OGRCreateCoordinateTransformation(&oSrcSRSHoriz, poGeog);
1495 1006 : if (poCT)
1496 : {
1497 1006 : poCT->SetEmitErrors(false);
1498 :
1499 : double x[4], y[4];
1500 1006 : x[0] = adfGT[0];
1501 1006 : y[0] = adfGT[3];
1502 1006 : x[1] = adfGT[0] + nXSize * adfGT[1];
1503 1006 : y[1] = adfGT[3];
1504 1006 : x[2] = adfGT[0];
1505 1006 : y[2] = adfGT[3] + nYSize * adfGT[5];
1506 1006 : x[3] = x[1];
1507 1006 : y[3] = y[2];
1508 1006 : int validity[4] = {false, false, false, false};
1509 1006 : poCT->Transform(4, x, y, nullptr, validity);
1510 1006 : dfWestLongitudeDeg = std::numeric_limits<double>::max();
1511 1006 : dfSouthLatitudeDeg = std::numeric_limits<double>::max();
1512 1006 : dfEastLongitudeDeg = -std::numeric_limits<double>::max();
1513 1006 : dfNorthLatitudeDeg = -std::numeric_limits<double>::max();
1514 5030 : for (int i = 0; i < 4; i++)
1515 : {
1516 4024 : if (validity[i])
1517 : {
1518 4012 : ret = true;
1519 4012 : dfWestLongitudeDeg = std::min(dfWestLongitudeDeg, x[i]);
1520 4012 : dfSouthLatitudeDeg = std::min(dfSouthLatitudeDeg, y[i]);
1521 4012 : dfEastLongitudeDeg = std::max(dfEastLongitudeDeg, x[i]);
1522 4012 : dfNorthLatitudeDeg = std::max(dfNorthLatitudeDeg, y[i]);
1523 : }
1524 : }
1525 1006 : if (validity[0] && validity[1] && x[0] > x[1])
1526 : {
1527 3 : dfWestLongitudeDeg = x[0];
1528 3 : dfEastLongitudeDeg = x[1];
1529 : }
1530 1006 : if (ret && std::fabs(dfWestLongitudeDeg) <= 180 &&
1531 1002 : std::fabs(dfEastLongitudeDeg) <= 180 &&
1532 998 : std::fabs(dfSouthLatitudeDeg) <= 90 &&
1533 994 : std::fabs(dfNorthLatitudeDeg) <= 90)
1534 : {
1535 994 : CPLDebug("GDAL", "Computing area of interest: %g, %g, %g, %g",
1536 : dfWestLongitudeDeg, dfSouthLatitudeDeg,
1537 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
1538 : }
1539 : else
1540 : {
1541 12 : CPLDebug("GDAL", "Could not compute area of interest");
1542 12 : dfWestLongitudeDeg = 0;
1543 12 : dfSouthLatitudeDeg = 0;
1544 12 : dfEastLongitudeDeg = 0;
1545 12 : dfNorthLatitudeDeg = 0;
1546 : }
1547 1006 : OGRCoordinateTransformation::DestroyCT(poCT);
1548 : }
1549 :
1550 1006 : delete poGeog;
1551 : }
1552 :
1553 1006 : return ret;
1554 : }
1555 :
1556 3 : bool GDALComputeAreaOfInterest(OGRSpatialReference *poSRS, double dfX1,
1557 : double dfY1, double dfX2, double dfY2,
1558 : double &dfWestLongitudeDeg,
1559 : double &dfSouthLatitudeDeg,
1560 : double &dfEastLongitudeDeg,
1561 : double &dfNorthLatitudeDeg)
1562 : {
1563 3 : bool ret = false;
1564 :
1565 3 : if (!poSRS)
1566 0 : return false;
1567 :
1568 3 : OGRSpatialReference oSrcSRSHoriz(*poSRS);
1569 3 : if (oSrcSRSHoriz.IsCompound())
1570 : {
1571 0 : oSrcSRSHoriz.StripVertical();
1572 : }
1573 :
1574 3 : OGRSpatialReference *poGeog = oSrcSRSHoriz.CloneGeogCS();
1575 3 : if (poGeog)
1576 : {
1577 3 : poGeog->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1578 :
1579 3 : auto poCT = OGRCreateCoordinateTransformation(&oSrcSRSHoriz, poGeog);
1580 3 : if (poCT)
1581 : {
1582 : double x[4], y[4];
1583 3 : x[0] = dfX1;
1584 3 : y[0] = dfY1;
1585 3 : x[1] = dfX2;
1586 3 : y[1] = dfY1;
1587 3 : x[2] = dfX1;
1588 3 : y[2] = dfY2;
1589 3 : x[3] = dfX2;
1590 3 : y[3] = dfY2;
1591 3 : int validity[4] = {false, false, false, false};
1592 3 : poCT->Transform(4, x, y, nullptr, validity);
1593 3 : dfWestLongitudeDeg = std::numeric_limits<double>::max();
1594 3 : dfSouthLatitudeDeg = std::numeric_limits<double>::max();
1595 3 : dfEastLongitudeDeg = -std::numeric_limits<double>::max();
1596 3 : dfNorthLatitudeDeg = -std::numeric_limits<double>::max();
1597 15 : for (int i = 0; i < 4; i++)
1598 : {
1599 12 : if (validity[i])
1600 : {
1601 12 : ret = true;
1602 12 : dfWestLongitudeDeg = std::min(dfWestLongitudeDeg, x[i]);
1603 12 : dfSouthLatitudeDeg = std::min(dfSouthLatitudeDeg, y[i]);
1604 12 : dfEastLongitudeDeg = std::max(dfEastLongitudeDeg, x[i]);
1605 12 : dfNorthLatitudeDeg = std::max(dfNorthLatitudeDeg, y[i]);
1606 : }
1607 : }
1608 3 : if (validity[0] && validity[1] && (dfX1 - dfX2) * (x[0] - x[1]) < 0)
1609 : {
1610 0 : dfWestLongitudeDeg = x[0];
1611 0 : dfEastLongitudeDeg = x[1];
1612 : }
1613 3 : if (ret)
1614 : {
1615 3 : CPLDebug("GDAL", "Computing area of interest: %g, %g, %g, %g",
1616 : dfWestLongitudeDeg, dfSouthLatitudeDeg,
1617 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
1618 : }
1619 : else
1620 : {
1621 0 : CPLDebug("GDAL", "Could not compute area of interest");
1622 0 : dfWestLongitudeDeg = 0;
1623 0 : dfSouthLatitudeDeg = 0;
1624 0 : dfEastLongitudeDeg = 0;
1625 0 : dfNorthLatitudeDeg = 0;
1626 : }
1627 3 : delete poCT;
1628 : }
1629 :
1630 3 : delete poGeog;
1631 : }
1632 :
1633 3 : return ret;
1634 : }
1635 :
1636 : /************************************************************************/
1637 : /* GDALGCPAntimeridianUnwrap() */
1638 : /************************************************************************/
1639 :
1640 : /* Deal with discontinuties of dfGCPX longitudes around the anti-meridian.
1641 : * Cf https://github.com/OSGeo/gdal/issues/8371
1642 : */
1643 40 : static void GDALGCPAntimeridianUnwrap(int nGCPCount, GDAL_GCP *pasGCPList,
1644 : const OGRSpatialReference &oSRS,
1645 : CSLConstList papszOptions)
1646 : {
1647 : const char *pszGCPAntimeridianUnwrap =
1648 40 : CSLFetchNameValueDef(papszOptions, "GCP_ANTIMERIDIAN_UNWRAP", "AUTO");
1649 119 : const bool bForced = EQUAL(pszGCPAntimeridianUnwrap, "YES") ||
1650 39 : EQUAL(pszGCPAntimeridianUnwrap, "ON") ||
1651 118 : EQUAL(pszGCPAntimeridianUnwrap, "TRUE") ||
1652 39 : EQUAL(pszGCPAntimeridianUnwrap, "1");
1653 46 : if (bForced || (!oSRS.IsEmpty() && oSRS.IsGeographic() &&
1654 6 : fabs(oSRS.GetAngularUnits(nullptr) -
1655 6 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-8 &&
1656 6 : EQUAL(pszGCPAntimeridianUnwrap, "AUTO")))
1657 : {
1658 6 : if (!bForced)
1659 : {
1660 : // Proceed to unwrapping only if the longitudes are within
1661 : // [-180, -170] or [170, 180]
1662 425 : for (int i = 0; i < nGCPCount; ++i)
1663 : {
1664 423 : const double dfLongAbs = fabs(pasGCPList[i].dfGCPX);
1665 423 : if (dfLongAbs > 180 || dfLongAbs < 170)
1666 : {
1667 3 : return;
1668 : }
1669 : }
1670 : }
1671 :
1672 3 : bool bDone = false;
1673 633 : for (int i = 0; i < nGCPCount; ++i)
1674 : {
1675 630 : if (pasGCPList[i].dfGCPX < 0)
1676 : {
1677 48 : if (!bDone)
1678 : {
1679 3 : bDone = true;
1680 3 : CPLDebug("WARP", "GCP longitude unwrapping");
1681 : }
1682 48 : pasGCPList[i].dfGCPX += 360;
1683 : }
1684 : }
1685 : }
1686 : }
1687 :
1688 : /************************************************************************/
1689 : /* GDALCreateGenImgProjTransformer2() */
1690 : /************************************************************************/
1691 :
1692 : /* clang-format off */
1693 : /**
1694 : * Create image to image transformer.
1695 : *
1696 : * This function creates a transformation object that maps from pixel/line
1697 : * coordinates on one image to pixel/line coordinates on another image. The
1698 : * images may potentially be georeferenced in different coordinate systems,
1699 : * and may used GCPs to map between their pixel/line coordinates and
1700 : * georeferenced coordinates (as opposed to the default assumption that their
1701 : * geotransform should be used).
1702 : *
1703 : * This transformer potentially performs three concatenated transformations.
1704 : *
1705 : * The first stage is from source image pixel/line coordinates to source
1706 : * image georeferenced coordinates, and may be done using the geotransform,
1707 : * or if not defined using a polynomial model derived from GCPs. If GCPs
1708 : * are used this stage is accomplished using GDALGCPTransform().
1709 : *
1710 : * The second stage is to change projections from the source coordinate system
1711 : * to the destination coordinate system, assuming they differ. This is
1712 : * accomplished internally using GDALReprojectionTransform().
1713 : *
1714 : * The third stage is converting from destination image georeferenced
1715 : * coordinates to destination image coordinates. This is done using the
1716 : * destination image geotransform, or if not available, using a polynomial
1717 : * model derived from GCPs. If GCPs are used this stage is accomplished using
1718 : * GDALGCPTransform(). This stage is skipped if hDstDS is NULL when the
1719 : * transformation is created.
1720 : *
1721 : * Supported Options (specified with the -to switch of gdalwarp for example):
1722 : * <ul>
1723 : * <li> SRC_SRS: WKT SRS, or any string recognized by
1724 : * OGRSpatialReference::SetFromUserInput(), to be used as an override for
1725 : * hSrcDS.</li>
1726 : * <li> DST_SRS: WKT SRS, or any string recognized by
1727 : * OGRSpatialReference::SetFromUserInput(), to be used as an override for
1728 : * hDstDS.
1729 : * </li>
1730 : * <li> COORDINATE_OPERATION: (GDAL >= 3.0) Coordinate operation, as
1731 : * a PROJ or WKT string, used as an override over the normally computed
1732 : * pipeline. The pipeline must take into account the axis order of the source
1733 : * and target SRS. <li> COORDINATE_EPOCH: (GDAL >= 3.0) Coordinate epoch,
1734 : * expressed as a decimal year. Useful for time-dependent coordinate operations.
1735 : * </li>
1736 : * <li> SRC_COORDINATE_EPOCH: (GDAL >= 3.4) Coordinate epoch of source CRS,
1737 : * expressed as a decimal year. Useful for time-dependent coordinate operations.
1738 : * </li>
1739 : * <li> DST_COORDINATE_EPOCH: (GDAL >= 3.4) Coordinate epoch of target CRS,
1740 : * expressed as a decimal year. Useful for time-dependent coordinate operations.
1741 : * </li>
1742 : * <li> GCPS_OK: If false, GCPs will not be used, default is TRUE.
1743 : * </li>
1744 : * <li> REFINE_MINIMUM_GCPS: The minimum amount of GCPs that should be available
1745 : * after the refinement.
1746 : * </li>
1747 : * <li> REFINE_TOLERANCE: The tolerance that specifies when a GCP will be
1748 : * eliminated.
1749 : * </li>
1750 : * <li> MAX_GCP_ORDER: the maximum order to use for GCP derived polynomials if
1751 : * possible. The default is to autoselect based on the number of GCPs.
1752 : * A value of -1 triggers use of Thin Plate Spline instead of polynomials.
1753 : * </li>
1754 : * <li>GCP_ANTIMERIDIAN_UNWRAP=AUTO/YES/NO. (GDAL >= 3.8) Whether to
1755 : * "unwrap" longitudes of ground control points that span the antimeridian.
1756 : * For datasets with GCPs in longitude/latitude coordinate space spanning the
1757 : * antimeridian, longitudes will have a discontinuity on +/- 180 deg, and
1758 : * will result in a subset of the GCPs with longitude in the [-180,-170] range
1759 : * and another subset in [170, 180]. By default (AUTO), that situation will be
1760 : * detected and longitudes in [-180,-170] will be shifted to [180, 190] to get
1761 : * a continuous set. This option can be set to YES to force that behavior
1762 : * (useful if no SRS information is available), or to NO to disable it.
1763 : * </li>
1764 : * <li> SRC_METHOD: may have a value which is one of GEOTRANSFORM,
1765 : * GCP_POLYNOMIAL, GCP_TPS, GEOLOC_ARRAY, RPC to force only one geolocation
1766 : * method to be considered on the source dataset. Will be used for pixel/line
1767 : * to georef transformation on the source dataset. NO_GEOTRANSFORM can be
1768 : * used to specify the identity geotransform (ungeoreference image)
1769 : * </li>
1770 : * <li> DST_METHOD: may have a value which is one of GEOTRANSFORM,
1771 : * GCP_POLYNOMIAL, GCP_TPS, GEOLOC_ARRAY (added in 3.5), RPC to force only one
1772 : * geolocation method to be considered on the target dataset. Will be used for
1773 : * pixel/line to georef transformation on the destination dataset.
1774 : * NO_GEOTRANSFORM can be used to specify the identity geotransform
1775 : * (ungeoreference image)
1776 : * </li>
1777 : * <li> RPC_HEIGHT: A fixed height to be used with RPC
1778 : * calculations.
1779 : * </li>
1780 : * <li> RPC_DEM: The name of a DEM file to be used with RPC
1781 : * calculations. See GDALCreateRPCTransformerV2() for more details.
1782 : * </li>
1783 : * <li> Other RPC related options. See GDALCreateRPCTransformerV2()
1784 : * </li>
1785 : * <li>
1786 : * INSERT_CENTER_LONG: May be set to FALSE to disable setting up a CENTER_LONG
1787 : * value on the coordinate system to rewrap things around the center of the
1788 : * image.
1789 : * </li>
1790 : * <li> SRC_APPROX_ERROR_IN_SRS_UNIT=err_threshold_in_SRS_units. (GDAL
1791 : * >= 2.2) Use an approximate transformer for the source transformer. Must be
1792 : * defined together with SRC_APPROX_ERROR_IN_PIXEL to be taken into account.
1793 : * </li>
1794 : * <li> SRC_APPROX_ERROR_IN_PIXEL=err_threshold_in_pixel. (GDAL >= 2.2) Use
1795 : * an approximate transformer for the source transformer.. Must be defined
1796 : * together with SRC_APPROX_ERROR_IN_SRS_UNIT to be taken into account.
1797 : * </li>
1798 : * <li>
1799 : * DST_APPROX_ERROR_IN_SRS_UNIT=err_threshold_in_SRS_units. (GDAL >= 2.2) Use
1800 : * an approximate transformer for the destination transformer. Must be defined
1801 : * together with DST_APPROX_ERROR_IN_PIXEL to be taken into account.
1802 : * </li>
1803 : * <li>
1804 : * DST_APPROX_ERROR_IN_PIXEL=err_threshold_in_pixel. (GDAL >= 2.2) Use an
1805 : * approximate transformer for the destination transformer. Must be defined
1806 : * together with DST_APPROX_ERROR_IN_SRS_UNIT to be taken into account.
1807 : * </li>
1808 : * <li>
1809 : * REPROJECTION_APPROX_ERROR_IN_SRC_SRS_UNIT=err_threshold_in_src_SRS_units.
1810 : * (GDAL >= 2.2) Use an approximate transformer for the coordinate
1811 : * reprojection. Must be used together with
1812 : * REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT to be taken into account.
1813 : * </li>
1814 : * <li>
1815 : * REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT=err_threshold_in_dst_SRS_units.
1816 : * (GDAL >= 2.2) Use an approximate transformer for the coordinate
1817 : * reprojection. Must be used together with
1818 : * REPROJECTION_APPROX_ERROR_IN_SRC_SRS_UNIT to be taken into account.
1819 : * </li>
1820 : * <li>
1821 : * AREA_OF_INTEREST=west_lon_deg,south_lat_deg,east_lon_deg,north_lat_deg. (GDAL
1822 : * >= 3.0) Area of interest, used to compute the best coordinate operation
1823 : * between the source and target SRS. If not specified, the bounding box of the
1824 : * source raster will be used.
1825 : * </li>
1826 : * <li> GEOLOC_BACKMAP_OVERSAMPLE_FACTOR=[0.1,2]. (GDAL >= 3.5) Oversample
1827 : * factor used to derive the size of the "backmap" used for geolocation array
1828 : * transformers. Default value is 1.3.
1829 : * </li>
1830 : * <li> GEOLOC_USE_TEMP_DATASETS=YES/NO.
1831 : * (GDAL >= 3.5) Whether temporary GeoTIFF datasets should be used to store
1832 : * the backmap. The default is NO, that is to use in-memory arrays, unless the
1833 : * number of pixels of the geolocation array is greater than 16 megapixels.
1834 : * </li>
1835 : * <li>
1836 : * GEOLOC_ARRAY/SRC_GEOLOC_ARRAY=filename. (GDAL >= 3.5.2) Name of a GDAL
1837 : * dataset containing a geolocation array and associated metadata. This is an
1838 : * alternative to having geolocation information described in the GEOLOCATION
1839 : * metadata domain of the source dataset. The dataset specified may have a
1840 : * GEOLOCATION metadata domain containing appropriate metadata, however default
1841 : * values are assigned for all omitted items. X_BAND defaults to 1 and Y_BAND to
1842 : * 2, however the dataset must contain exactly 2 bands. PIXEL_OFFSET and
1843 : * LINE_OFFSET default to 0. PIXEL_STEP and LINE_STEP default to the ratio of
1844 : * the width/height of the source dataset divided by the with/height of the
1845 : * geolocation array. SRS defaults to the geolocation array dataset's spatial
1846 : * reference system if set, otherwise WGS84 is used.
1847 : * GEOREFERENCING_CONVENTION is selected from the main metadata domain if it
1848 : * is omitted from the GEOLOCATION domain, and if not available
1849 : * TOP_LEFT_CORNER is assigned as a default.
1850 : * If GEOLOC_ARRAY is set SRC_METHOD
1851 : * defaults to GEOLOC_ARRAY.
1852 : * </li>
1853 : * <li>DST_GEOLOC_ARRAY=filename. (GDAL >= 3.5.2) Name of a
1854 : * GDAL dataset that contains at least 2 bands with the X and Y geolocation
1855 : * bands. This is an alternative to having geolocation information described in
1856 : * the GEOLOCATION metadata domain of the destination dataset. See
1857 : * SRC_GEOLOC_ARRAY description for details, assumptions, and defaults. If this
1858 : * option is set, DST_METHOD=GEOLOC_ARRAY will be assumed if not set.
1859 : * </li>
1860 : * </ul>
1861 : *
1862 : * The use case for the *_APPROX_ERROR_* options is when defining an approximate
1863 : * transformer on top of the GenImgProjTransformer globally is not practical.
1864 : * Such a use case is when the source dataset has RPC with a RPC DEM. In such
1865 : * case we don't want to use the approximate transformer on the RPC
1866 : * transformation, as the RPC DEM generally involves non-linearities that the
1867 : * approximate transformer will not detect. In such case, we must a
1868 : * non-approximated GenImgProjTransformer, but it might be worthwhile to use
1869 : * approximate sub- transformers, for example on coordinate reprojection. For
1870 : * example if warping from a source dataset with RPC to a destination dataset
1871 : * with a UTM projection, since the inverse UTM transformation is rather costly.
1872 : * In which case, one can use the REPROJECTION_APPROX_ERROR_IN_SRC_SRS_UNIT and
1873 : * REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT options.
1874 : *
1875 : * @param hSrcDS source dataset, or NULL.
1876 : * @param hDstDS destination dataset (or NULL).
1877 : * @param papszOptions NULL-terminated list of string options (or NULL).
1878 : *
1879 : * @return handle suitable for use GDALGenImgProjTransform(), and to be
1880 : * deallocated with GDALDestroyGenImgProjTransformer() or NULL on failure.
1881 : */
1882 : /* clang-format on */
1883 :
1884 1413 : void *GDALCreateGenImgProjTransformer2(GDALDatasetH hSrcDS, GDALDatasetH hDstDS,
1885 : char **papszOptions)
1886 :
1887 : {
1888 1413 : CSLConstList papszMD = nullptr;
1889 : GDALRPCInfoV2 sRPCInfo;
1890 1413 : const char *pszMethod = CSLFetchNameValue(papszOptions, "SRC_METHOD");
1891 1413 : if (pszMethod == nullptr)
1892 1402 : pszMethod = CSLFetchNameValue(papszOptions, "METHOD");
1893 1413 : const char *pszSrcSRS = CSLFetchNameValue(papszOptions, "SRC_SRS");
1894 1413 : const char *pszDstSRS = CSLFetchNameValue(papszOptions, "DST_SRS");
1895 :
1896 1413 : const char *pszValue = CSLFetchNameValue(papszOptions, "MAX_GCP_ORDER");
1897 1413 : const int nOrder = pszValue ? atoi(pszValue) : 0;
1898 :
1899 1413 : pszValue = CSLFetchNameValue(papszOptions, "GCPS_OK");
1900 1413 : const bool bGCPUseOK = pszValue ? CPLTestBool(pszValue) : true;
1901 :
1902 1413 : pszValue = CSLFetchNameValue(papszOptions, "REFINE_MINIMUM_GCPS");
1903 1413 : const int nMinimumGcps = pszValue ? atoi(pszValue) : -1;
1904 :
1905 1413 : pszValue = CSLFetchNameValue(papszOptions, "REFINE_TOLERANCE");
1906 1413 : const bool bRefine = pszValue != nullptr;
1907 1413 : const double dfTolerance = pszValue ? CPLAtof(pszValue) : 0.0;
1908 :
1909 1413 : double dfWestLongitudeDeg = 0.0;
1910 1413 : double dfSouthLatitudeDeg = 0.0;
1911 1413 : double dfEastLongitudeDeg = 0.0;
1912 1413 : double dfNorthLatitudeDeg = 0.0;
1913 1413 : bool bHasAreaOfInterest = false;
1914 1413 : pszValue = CSLFetchNameValue(papszOptions, "AREA_OF_INTEREST");
1915 1413 : if (pszValue)
1916 : {
1917 0 : char **papszTokens = CSLTokenizeString2(pszValue, ", ", 0);
1918 0 : if (CSLCount(papszTokens) == 4)
1919 : {
1920 0 : dfWestLongitudeDeg = CPLAtof(papszTokens[0]);
1921 0 : dfSouthLatitudeDeg = CPLAtof(papszTokens[1]);
1922 0 : dfEastLongitudeDeg = CPLAtof(papszTokens[2]);
1923 0 : dfNorthLatitudeDeg = CPLAtof(papszTokens[3]);
1924 0 : bHasAreaOfInterest = true;
1925 : }
1926 0 : CSLDestroy(papszTokens);
1927 : }
1928 :
1929 1413 : const char *pszCO = CSLFetchNameValue(papszOptions, "COORDINATE_OPERATION");
1930 :
1931 : const auto SetAxisMapping =
1932 3160 : [papszOptions](OGRSpatialReference &oSRS, const char *pszPrefix)
1933 : {
1934 1054 : const char *pszMapping = CSLFetchNameValue(
1935 2108 : papszOptions, std::string(pszPrefix)
1936 1054 : .append("_DATA_AXIS_TO_SRS_AXIS_MAPPING")
1937 : .c_str());
1938 1054 : if (pszMapping)
1939 : {
1940 4 : CPLStringList aosTokens(CSLTokenizeString2(pszMapping, ",", 0));
1941 4 : std::vector<int> anMapping;
1942 6 : for (int i = 0; i < aosTokens.size(); ++i)
1943 4 : anMapping.push_back(atoi(aosTokens[i]));
1944 2 : oSRS.SetDataAxisToSRSAxisMapping(anMapping);
1945 : }
1946 : else
1947 : {
1948 1052 : const char *pszStrategy = CSLFetchNameValueDef(
1949 : papszOptions,
1950 2104 : std::string(pszPrefix).append("_AXIS_MAPPING_STRATEGY").c_str(),
1951 : "TRADITIONAL_GIS_ORDER");
1952 1052 : if (EQUAL(pszStrategy, "TRADITIONAL_GIS_ORDER"))
1953 1051 : oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1954 1 : else if (EQUAL(pszStrategy, "AUTHORITY_COMPLIANT"))
1955 1 : oSRS.SetAxisMappingStrategy(OAMS_AUTHORITY_COMPLIANT);
1956 : else
1957 : {
1958 0 : CPLError(CE_Warning, CPLE_AppDefined,
1959 : "Unrecognized value '%s' for %s", pszStrategy,
1960 0 : std::string(pszPrefix)
1961 0 : .append("_AXIS_MAPPING_STRATEGY")
1962 : .c_str());
1963 0 : return false;
1964 : }
1965 : }
1966 1054 : return true;
1967 1413 : };
1968 :
1969 2826 : OGRSpatialReference oSrcSRS;
1970 1413 : if (pszSrcSRS)
1971 : {
1972 234 : if (pszSrcSRS[0] != '\0' &&
1973 116 : oSrcSRS.SetFromUserInput(pszSrcSRS) != OGRERR_NONE)
1974 : {
1975 0 : CPLError(CE_Failure, CPLE_AppDefined,
1976 : "Failed to import coordinate system `%s'.", pszSrcSRS);
1977 0 : return nullptr;
1978 : }
1979 118 : if (!SetAxisMapping(oSrcSRS, "SRC_SRS"))
1980 0 : return nullptr;
1981 : }
1982 :
1983 2826 : OGRSpatialReference oDstSRS;
1984 1413 : if (pszDstSRS)
1985 : {
1986 1870 : if (pszDstSRS[0] != '\0' &&
1987 934 : oDstSRS.SetFromUserInput(pszDstSRS) != OGRERR_NONE)
1988 : {
1989 0 : CPLError(CE_Failure, CPLE_AppDefined,
1990 : "Failed to import coordinate system `%s'.", pszDstSRS);
1991 0 : return nullptr;
1992 : }
1993 936 : if (!SetAxisMapping(oDstSRS, "DST_SRS"))
1994 0 : return nullptr;
1995 : }
1996 :
1997 : const char *pszSrcGeolocArray =
1998 1413 : CSLFetchNameValueDef(papszOptions, "SRC_GEOLOC_ARRAY",
1999 : CSLFetchNameValue(papszOptions, "GEOLOC_ARRAY"));
2000 1413 : if (pszMethod == nullptr && pszSrcGeolocArray != nullptr)
2001 7 : pszMethod = "GEOLOC_ARRAY";
2002 :
2003 : /* -------------------------------------------------------------------- */
2004 : /* Initialize the transform info. */
2005 : /* -------------------------------------------------------------------- */
2006 : GDALGenImgProjTransformInfo *psInfo =
2007 1413 : GDALCreateGenImgProjTransformerInternal();
2008 :
2009 1413 : bool bCanUseSrcGeoTransform = false;
2010 :
2011 : /* -------------------------------------------------------------------- */
2012 : /* Get forward and inverse geotransform for the source image. */
2013 : /* -------------------------------------------------------------------- */
2014 1413 : if (hSrcDS == nullptr ||
2015 73 : (pszMethod != nullptr && EQUAL(pszMethod, "NO_GEOTRANSFORM")))
2016 : {
2017 89 : psInfo->adfSrcGeoTransform[0] = 0.0;
2018 89 : psInfo->adfSrcGeoTransform[1] = 1.0;
2019 89 : psInfo->adfSrcGeoTransform[2] = 0.0;
2020 89 : psInfo->adfSrcGeoTransform[3] = 0.0;
2021 89 : psInfo->adfSrcGeoTransform[4] = 0.0;
2022 89 : psInfo->adfSrcGeoTransform[5] = 1.0;
2023 89 : memcpy(psInfo->adfSrcInvGeoTransform, psInfo->adfSrcGeoTransform,
2024 : sizeof(double) * 6);
2025 : }
2026 2586 : else if ((pszMethod == nullptr || EQUAL(pszMethod, "GEOTRANSFORM")) &&
2027 1262 : GDALGetGeoTransform(hSrcDS, psInfo->adfSrcGeoTransform) == CE_None)
2028 : {
2029 1190 : if (!GDALInvGeoTransform(psInfo->adfSrcGeoTransform,
2030 1190 : psInfo->adfSrcInvGeoTransform))
2031 : {
2032 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
2033 0 : GDALDestroyGenImgProjTransformer(psInfo);
2034 0 : return nullptr;
2035 : }
2036 1190 : if (pszSrcSRS == nullptr)
2037 : {
2038 1117 : auto hSRS = GDALGetSpatialRef(hSrcDS);
2039 1117 : if (hSRS)
2040 922 : oSrcSRS = *(OGRSpatialReference::FromHandle(hSRS));
2041 : }
2042 1190 : if (!bHasAreaOfInterest && pszCO == nullptr && !oSrcSRS.IsEmpty())
2043 : {
2044 990 : GDALComputeAreaOfInterest(&oSrcSRS, psInfo->adfSrcGeoTransform,
2045 : GDALGetRasterXSize(hSrcDS),
2046 : GDALGetRasterYSize(hSrcDS),
2047 : dfWestLongitudeDeg, dfSouthLatitudeDeg,
2048 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
2049 : }
2050 1190 : bCanUseSrcGeoTransform = true;
2051 : }
2052 134 : else if (bGCPUseOK &&
2053 62 : (pszMethod == nullptr || EQUAL(pszMethod, "GCP_POLYNOMIAL")) &&
2054 268 : GDALGetGCPCount(hSrcDS) > 0 && nOrder >= 0)
2055 : {
2056 29 : if (pszSrcSRS == nullptr)
2057 : {
2058 28 : auto hSRS = GDALGetGCPSpatialRef(hSrcDS);
2059 28 : if (hSRS)
2060 28 : oSrcSRS = *(OGRSpatialReference::FromHandle(hSRS));
2061 : }
2062 :
2063 29 : const auto nGCPCount = GDALGetGCPCount(hSrcDS);
2064 29 : auto pasGCPList = GDALDuplicateGCPs(nGCPCount, GDALGetGCPs(hSrcDS));
2065 29 : GDALGCPAntimeridianUnwrap(nGCPCount, pasGCPList, oSrcSRS, papszOptions);
2066 :
2067 29 : if (bRefine)
2068 : {
2069 0 : psInfo->pSrcTransformArg = GDALCreateGCPRefineTransformer(
2070 : nGCPCount, pasGCPList, nOrder, FALSE, dfTolerance,
2071 : nMinimumGcps);
2072 : }
2073 : else
2074 : {
2075 29 : psInfo->pSrcTransformArg =
2076 29 : GDALCreateGCPTransformer(nGCPCount, pasGCPList, nOrder, FALSE);
2077 : }
2078 :
2079 29 : GDALDeinitGCPs(nGCPCount, pasGCPList);
2080 29 : CPLFree(pasGCPList);
2081 :
2082 29 : if (psInfo->pSrcTransformArg == nullptr)
2083 : {
2084 0 : GDALDestroyGenImgProjTransformer(psInfo);
2085 0 : return nullptr;
2086 : }
2087 29 : psInfo->pSrcTransformer = GDALGCPTransform;
2088 : }
2089 :
2090 116 : else if (bGCPUseOK && GDALGetGCPCount(hSrcDS) > 0 && nOrder <= 0 &&
2091 11 : (pszMethod == nullptr || EQUAL(pszMethod, "GCP_TPS")))
2092 : {
2093 11 : if (pszSrcSRS == nullptr)
2094 : {
2095 11 : auto hSRS = GDALGetGCPSpatialRef(hSrcDS);
2096 11 : if (hSRS)
2097 10 : oSrcSRS = *(OGRSpatialReference::FromHandle(hSRS));
2098 : }
2099 :
2100 11 : const auto nGCPCount = GDALGetGCPCount(hSrcDS);
2101 11 : auto pasGCPList = GDALDuplicateGCPs(nGCPCount, GDALGetGCPs(hSrcDS));
2102 11 : GDALGCPAntimeridianUnwrap(nGCPCount, pasGCPList, oSrcSRS, papszOptions);
2103 :
2104 11 : psInfo->pSrcTransformArg = GDALCreateTPSTransformerInt(
2105 : nGCPCount, pasGCPList, FALSE, papszOptions);
2106 :
2107 11 : GDALDeinitGCPs(nGCPCount, pasGCPList);
2108 11 : CPLFree(pasGCPList);
2109 :
2110 11 : if (psInfo->pSrcTransformArg == nullptr)
2111 : {
2112 2 : GDALDestroyGenImgProjTransformer(psInfo);
2113 2 : return nullptr;
2114 : }
2115 9 : psInfo->pSrcTransformer = GDALTPSTransform;
2116 : }
2117 :
2118 48 : else if ((pszMethod == nullptr || EQUAL(pszMethod, "RPC")) &&
2119 185 : (papszMD = GDALGetMetadata(hSrcDS, "RPC")) != nullptr &&
2120 43 : GDALExtractRPCInfoV2(papszMD, &sRPCInfo))
2121 : {
2122 43 : psInfo->pSrcTransformArg =
2123 43 : GDALCreateRPCTransformerV2(&sRPCInfo, FALSE, 0, papszOptions);
2124 43 : if (psInfo->pSrcTransformArg == nullptr)
2125 : {
2126 1 : GDALDestroyGenImgProjTransformer(psInfo);
2127 1 : return nullptr;
2128 : }
2129 42 : psInfo->pSrcTransformer = GDALRPCTransform;
2130 42 : if (pszSrcSRS == nullptr)
2131 : {
2132 42 : oSrcSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
2133 42 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2134 : }
2135 : }
2136 :
2137 102 : else if ((pszMethod == nullptr || EQUAL(pszMethod, "GEOLOC_ARRAY")) &&
2138 51 : ((papszMD = GDALGetMetadata(hSrcDS, "GEOLOCATION")) != nullptr ||
2139 : pszSrcGeolocArray != nullptr))
2140 : {
2141 44 : CPLStringList aosGeolocMD; // keep in this scope
2142 44 : if (pszSrcGeolocArray != nullptr)
2143 : {
2144 7 : if (papszMD != nullptr)
2145 : {
2146 0 : CPLError(
2147 : CE_Warning, CPLE_AppDefined,
2148 : "Both GEOLOCATION metadata domain on the source dataset "
2149 : "and [SRC_]GEOLOC_ARRAY transformer option are set. "
2150 : "Only using the later.");
2151 : }
2152 : aosGeolocMD =
2153 14 : GDALCreateGeolocationMetadata(hSrcDS, pszSrcGeolocArray,
2154 7 : /* bIsSource= */ true);
2155 7 : if (aosGeolocMD.empty())
2156 : {
2157 2 : GDALDestroyGenImgProjTransformer(psInfo);
2158 2 : return nullptr;
2159 : }
2160 5 : papszMD = aosGeolocMD.List();
2161 : }
2162 :
2163 42 : psInfo->pSrcTransformArg = GDALCreateGeoLocTransformerEx(
2164 : hSrcDS, papszMD, FALSE, nullptr, papszOptions);
2165 42 : if (psInfo->pSrcTransformArg == nullptr)
2166 : {
2167 1 : GDALDestroyGenImgProjTransformer(psInfo);
2168 1 : return nullptr;
2169 : }
2170 41 : psInfo->pSrcTransformer = GDALGeoLocTransform;
2171 41 : if (pszSrcSRS == nullptr)
2172 : {
2173 41 : pszSrcSRS = CSLFetchNameValue(papszMD, "SRS");
2174 41 : if (pszSrcSRS)
2175 : {
2176 39 : oSrcSRS.SetFromUserInput(pszSrcSRS);
2177 39 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2178 : }
2179 : }
2180 : }
2181 :
2182 7 : else if (pszMethod != nullptr)
2183 : {
2184 0 : CPLError(CE_Failure, CPLE_AppDefined,
2185 : "Unable to compute a %s based transformation between "
2186 : "pixel/line and georeferenced coordinates for %s.",
2187 : pszMethod, GDALGetDescription(hSrcDS));
2188 :
2189 0 : GDALDestroyGenImgProjTransformer(psInfo);
2190 0 : return nullptr;
2191 : }
2192 :
2193 : else
2194 : {
2195 7 : CPLError(CE_Failure, CPLE_AppDefined,
2196 : "Unable to compute a transformation between pixel/line "
2197 : "and georeferenced coordinates for %s. "
2198 : "There is no affine transformation and no GCPs. "
2199 : "Specify transformation option SRC_METHOD=NO_GEOTRANSFORM "
2200 : "to bypass this check.",
2201 : GDALGetDescription(hSrcDS));
2202 :
2203 7 : GDALDestroyGenImgProjTransformer(psInfo);
2204 7 : return nullptr;
2205 : }
2206 :
2207 : /* -------------------------------------------------------------------- */
2208 : /* Handle optional source approximation transformer. */
2209 : /* -------------------------------------------------------------------- */
2210 1400 : if (psInfo->pSrcTransformer)
2211 : {
2212 : const char *pszSrcApproxErrorFwd =
2213 121 : CSLFetchNameValue(papszOptions, "SRC_APPROX_ERROR_IN_SRS_UNIT");
2214 : const char *pszSrcApproxErrorReverse =
2215 121 : CSLFetchNameValue(papszOptions, "SRC_APPROX_ERROR_IN_PIXEL");
2216 121 : if (pszSrcApproxErrorFwd && pszSrcApproxErrorReverse)
2217 : {
2218 1 : void *pArg = GDALCreateApproxTransformer2(
2219 : psInfo->pSrcTransformer, psInfo->pSrcTransformArg,
2220 : CPLAtof(pszSrcApproxErrorFwd),
2221 : CPLAtof(pszSrcApproxErrorReverse));
2222 1 : if (pArg == nullptr)
2223 : {
2224 0 : GDALDestroyGenImgProjTransformer(psInfo);
2225 0 : return nullptr;
2226 : }
2227 1 : psInfo->pSrcTransformArg = pArg;
2228 1 : psInfo->pSrcTransformer = GDALApproxTransform;
2229 1 : GDALApproxTransformerOwnsSubtransformer(psInfo->pSrcTransformArg,
2230 : TRUE);
2231 : }
2232 : }
2233 :
2234 : /* -------------------------------------------------------------------- */
2235 : /* Get forward and inverse geotransform for destination image. */
2236 : /* If we have no destination use a unit transform. */
2237 : /* -------------------------------------------------------------------- */
2238 1400 : const char *pszDstMethod = CSLFetchNameValue(papszOptions, "DST_METHOD");
2239 : const char *pszDstGeolocArray =
2240 1400 : CSLFetchNameValue(papszOptions, "DST_GEOLOC_ARRAY");
2241 1400 : if (pszDstMethod == nullptr && pszDstGeolocArray != nullptr)
2242 2 : pszDstMethod = "GEOLOC_ARRAY";
2243 :
2244 1400 : if (!hDstDS ||
2245 5 : (pszDstMethod != nullptr && EQUAL(pszDstMethod, "NO_GEOTRANSFORM")))
2246 : {
2247 1051 : psInfo->adfDstGeoTransform[0] = 0.0;
2248 1051 : psInfo->adfDstGeoTransform[1] = 1.0;
2249 1051 : psInfo->adfDstGeoTransform[2] = 0.0;
2250 1051 : psInfo->adfDstGeoTransform[3] = 0.0;
2251 1051 : psInfo->adfDstGeoTransform[4] = 0.0;
2252 1051 : psInfo->adfDstGeoTransform[5] = 1.0;
2253 1051 : memcpy(psInfo->adfDstInvGeoTransform, psInfo->adfDstGeoTransform,
2254 : sizeof(double) * 6);
2255 : }
2256 696 : else if ((pszDstMethod == nullptr || EQUAL(pszDstMethod, "GEOTRANSFORM")) &&
2257 347 : GDALGetGeoTransform(hDstDS, psInfo->adfDstGeoTransform) == CE_None)
2258 : {
2259 340 : if (pszDstSRS == nullptr)
2260 : {
2261 258 : auto hSRS = GDALGetSpatialRef(hDstDS);
2262 258 : if (hSRS)
2263 175 : oDstSRS = *(OGRSpatialReference::FromHandle(hSRS));
2264 : }
2265 340 : if (!GDALInvGeoTransform(psInfo->adfDstGeoTransform,
2266 340 : psInfo->adfDstInvGeoTransform))
2267 : {
2268 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
2269 0 : GDALDestroyGenImgProjTransformer(psInfo);
2270 0 : return nullptr;
2271 : }
2272 : }
2273 9 : else if (bGCPUseOK &&
2274 2 : (pszDstMethod == nullptr ||
2275 2 : EQUAL(pszDstMethod, "GCP_POLYNOMIAL")) &&
2276 18 : GDALGetGCPCount(hDstDS) > 0 && nOrder >= 0)
2277 : {
2278 0 : if (pszDstSRS == nullptr)
2279 : {
2280 0 : auto hSRS = GDALGetGCPSpatialRef(hDstDS);
2281 0 : if (hSRS)
2282 0 : oDstSRS = *(OGRSpatialReference::FromHandle(hSRS));
2283 : }
2284 :
2285 0 : const auto nGCPCount = GDALGetGCPCount(hDstDS);
2286 0 : auto pasGCPList = GDALDuplicateGCPs(nGCPCount, GDALGetGCPs(hDstDS));
2287 0 : GDALGCPAntimeridianUnwrap(nGCPCount, pasGCPList, oDstSRS, papszOptions);
2288 :
2289 0 : if (bRefine)
2290 : {
2291 0 : psInfo->pDstTransformArg = GDALCreateGCPRefineTransformer(
2292 : nGCPCount, pasGCPList, nOrder, FALSE, dfTolerance,
2293 : nMinimumGcps);
2294 : }
2295 : else
2296 : {
2297 0 : psInfo->pDstTransformArg =
2298 0 : GDALCreateGCPTransformer(nGCPCount, pasGCPList, nOrder, FALSE);
2299 : }
2300 :
2301 0 : GDALDeinitGCPs(nGCPCount, pasGCPList);
2302 0 : CPLFree(pasGCPList);
2303 :
2304 0 : if (psInfo->pDstTransformArg == nullptr)
2305 : {
2306 0 : GDALDestroyGenImgProjTransformer(psInfo);
2307 0 : return nullptr;
2308 : }
2309 0 : psInfo->pDstTransformer = GDALGCPTransform;
2310 : }
2311 9 : else if (bGCPUseOK && GDALGetGCPCount(hDstDS) > 0 && nOrder <= 0 &&
2312 0 : (pszDstMethod == nullptr || EQUAL(pszDstMethod, "GCP_TPS")))
2313 : {
2314 0 : if (pszDstSRS == nullptr)
2315 : {
2316 0 : auto hSRS = GDALGetGCPSpatialRef(hDstDS);
2317 0 : if (hSRS)
2318 0 : oDstSRS = *(OGRSpatialReference::FromHandle(hSRS));
2319 : }
2320 :
2321 0 : const auto nGCPCount = GDALGetGCPCount(hDstDS);
2322 0 : auto pasGCPList = GDALDuplicateGCPs(nGCPCount, GDALGetGCPs(hDstDS));
2323 0 : GDALGCPAntimeridianUnwrap(nGCPCount, pasGCPList, oDstSRS, papszOptions);
2324 :
2325 0 : psInfo->pDstTransformArg = GDALCreateTPSTransformerInt(
2326 : nGCPCount, pasGCPList, FALSE, papszOptions);
2327 :
2328 0 : GDALDeinitGCPs(nGCPCount, pasGCPList);
2329 0 : CPLFree(pasGCPList);
2330 :
2331 0 : if (psInfo->pDstTransformArg == nullptr)
2332 : {
2333 0 : GDALDestroyGenImgProjTransformer(psInfo);
2334 0 : return nullptr;
2335 : }
2336 0 : psInfo->pDstTransformer = GDALTPSTransform;
2337 : }
2338 2 : else if ((pszDstMethod == nullptr || EQUAL(pszDstMethod, "RPC")) &&
2339 13 : (papszMD = GDALGetMetadata(hDstDS, "RPC")) != nullptr &&
2340 2 : GDALExtractRPCInfoV2(papszMD, &sRPCInfo))
2341 : {
2342 2 : psInfo->pDstTransformArg =
2343 2 : GDALCreateRPCTransformerV2(&sRPCInfo, FALSE, 0, papszOptions);
2344 2 : if (psInfo->pDstTransformArg == nullptr)
2345 : {
2346 0 : GDALDestroyGenImgProjTransformer(psInfo);
2347 0 : return nullptr;
2348 : }
2349 2 : psInfo->pDstTransformer = GDALRPCTransform;
2350 2 : if (pszDstSRS == nullptr)
2351 : {
2352 2 : oDstSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
2353 2 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2354 : }
2355 : }
2356 14 : else if ((pszDstMethod == nullptr || EQUAL(pszDstMethod, "GEOLOC_ARRAY")) &&
2357 7 : ((papszMD = GDALGetMetadata(hDstDS, "GEOLOCATION")) != nullptr ||
2358 : pszDstGeolocArray != nullptr))
2359 : {
2360 7 : CPLStringList aosGeolocMD; // keep in this scope
2361 7 : if (pszDstGeolocArray != nullptr)
2362 : {
2363 2 : if (papszMD != nullptr)
2364 : {
2365 0 : CPLError(
2366 : CE_Warning, CPLE_AppDefined,
2367 : "Both GEOLOCATION metadata domain on the target dataset "
2368 : "and DST_GEOLOC_ARRAY transformer option are set. "
2369 : "Only using the later.");
2370 : }
2371 : aosGeolocMD =
2372 4 : GDALCreateGeolocationMetadata(hDstDS, pszDstGeolocArray,
2373 2 : /* bIsSource= */ false);
2374 2 : if (aosGeolocMD.empty())
2375 : {
2376 1 : GDALDestroyGenImgProjTransformer(psInfo);
2377 1 : return nullptr;
2378 : }
2379 1 : papszMD = aosGeolocMD.List();
2380 : }
2381 :
2382 6 : psInfo->pDstTransformArg = GDALCreateGeoLocTransformerEx(
2383 : hDstDS, papszMD, FALSE, nullptr, papszOptions);
2384 6 : if (psInfo->pDstTransformArg == nullptr)
2385 : {
2386 1 : GDALDestroyGenImgProjTransformer(psInfo);
2387 1 : return nullptr;
2388 : }
2389 5 : psInfo->pDstTransformer = GDALGeoLocTransform;
2390 5 : if (pszDstSRS == nullptr)
2391 : {
2392 5 : pszDstSRS = CSLFetchNameValue(papszMD, "SRS");
2393 5 : if (pszDstSRS)
2394 : {
2395 5 : oDstSRS.SetFromUserInput(pszDstSRS);
2396 5 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2397 : }
2398 : }
2399 : }
2400 : else
2401 : {
2402 0 : CPLError(CE_Failure, CPLE_AppDefined,
2403 : "Unable to compute a transformation between pixel/line "
2404 : "and georeferenced coordinates for %s. "
2405 : "There is no affine transformation and no GCPs. "
2406 : "Specify transformation option DST_METHOD=NO_GEOTRANSFORM "
2407 : "to bypass this check.",
2408 : GDALGetDescription(hDstDS));
2409 :
2410 0 : GDALDestroyGenImgProjTransformer(psInfo);
2411 0 : return nullptr;
2412 : }
2413 :
2414 : /* -------------------------------------------------------------------- */
2415 : /* Handle optional destination approximation transformer. */
2416 : /* -------------------------------------------------------------------- */
2417 1398 : if (psInfo->pDstTransformer)
2418 : {
2419 : const char *pszDstApproxErrorFwd =
2420 7 : CSLFetchNameValue(papszOptions, "DST_APPROX_ERROR_IN_PIXEL");
2421 : const char *pszDstApproxErrorReverse =
2422 7 : CSLFetchNameValue(papszOptions, "DST_APPROX_ERROR_IN_SRS_UNIT");
2423 7 : if (pszDstApproxErrorFwd && pszDstApproxErrorReverse)
2424 : {
2425 0 : void *pArg = GDALCreateApproxTransformer2(
2426 : psInfo->pDstTransformer, psInfo->pDstTransformArg,
2427 : CPLAtof(pszDstApproxErrorFwd),
2428 : CPLAtof(pszDstApproxErrorReverse));
2429 0 : if (pArg == nullptr)
2430 : {
2431 0 : GDALDestroyGenImgProjTransformer(psInfo);
2432 0 : return nullptr;
2433 : }
2434 0 : psInfo->pDstTransformArg = pArg;
2435 0 : psInfo->pDstTransformer = GDALApproxTransform;
2436 0 : GDALApproxTransformerOwnsSubtransformer(psInfo->pDstTransformArg,
2437 : TRUE);
2438 : }
2439 : }
2440 :
2441 : /* -------------------------------------------------------------------- */
2442 : /* Setup reprojection. */
2443 : /* -------------------------------------------------------------------- */
2444 :
2445 1398 : if (CPLFetchBool(papszOptions, "STRIP_VERT_CS", false))
2446 : {
2447 0 : if (oSrcSRS.IsCompound())
2448 : {
2449 0 : oSrcSRS.StripVertical();
2450 : }
2451 0 : if (oDstSRS.IsCompound())
2452 : {
2453 0 : oDstSRS.StripVertical();
2454 : }
2455 : }
2456 :
2457 : const bool bMayInsertCenterLong =
2458 2390 : (bCanUseSrcGeoTransform && !oSrcSRS.IsEmpty() && hSrcDS &&
2459 992 : CPLFetchBool(papszOptions, "INSERT_CENTER_LONG", true));
2460 : const char *pszSrcCoordEpoch =
2461 1398 : CSLFetchNameValue(papszOptions, "SRC_COORDINATE_EPOCH");
2462 : const char *pszDstCoordEpoch =
2463 1398 : CSLFetchNameValue(papszOptions, "DST_COORDINATE_EPOCH");
2464 2551 : if ((!oSrcSRS.IsEmpty() && !oDstSRS.IsEmpty() &&
2465 1052 : (pszSrcCoordEpoch || pszDstCoordEpoch || !oSrcSRS.IsSame(&oDstSRS) ||
2466 2551 : (oSrcSRS.IsGeographic() && bMayInsertCenterLong))) ||
2467 : pszCO)
2468 : {
2469 732 : CPLStringList aosOptions;
2470 :
2471 732 : if (bMayInsertCenterLong)
2472 : {
2473 715 : InsertCenterLong(hSrcDS, &oSrcSRS, aosOptions);
2474 : }
2475 :
2476 732 : if (CPLFetchBool(papszOptions, "PROMOTE_TO_3D", false))
2477 : {
2478 19 : oSrcSRS.PromoteTo3D(nullptr);
2479 19 : oDstSRS.PromoteTo3D(nullptr);
2480 : }
2481 :
2482 732 : if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
2483 29 : dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))
2484 : {
2485 : aosOptions.SetNameValue(
2486 : "AREA_OF_INTEREST",
2487 : CPLSPrintf("%.16g,%.16g,%.16g,%.16g", dfWestLongitudeDeg,
2488 : dfSouthLatitudeDeg, dfEastLongitudeDeg,
2489 703 : dfNorthLatitudeDeg));
2490 : }
2491 732 : if (pszCO)
2492 : {
2493 7 : aosOptions.SetNameValue("COORDINATE_OPERATION", pszCO);
2494 : }
2495 :
2496 : const char *pszCoordEpoch =
2497 732 : CSLFetchNameValue(papszOptions, "COORDINATE_EPOCH");
2498 732 : if (pszCoordEpoch)
2499 : {
2500 1 : aosOptions.SetNameValue("COORDINATE_EPOCH", pszCoordEpoch);
2501 : }
2502 :
2503 732 : if (pszSrcCoordEpoch)
2504 : {
2505 0 : aosOptions.SetNameValue("SRC_COORDINATE_EPOCH", pszSrcCoordEpoch);
2506 0 : oSrcSRS.SetCoordinateEpoch(CPLAtof(pszSrcCoordEpoch));
2507 : }
2508 :
2509 732 : if (pszDstCoordEpoch)
2510 : {
2511 0 : aosOptions.SetNameValue("DST_COORDINATE_EPOCH", pszDstCoordEpoch);
2512 0 : oDstSRS.SetCoordinateEpoch(CPLAtof(pszDstCoordEpoch));
2513 : }
2514 :
2515 736 : psInfo->pReprojectArg = GDALCreateReprojectionTransformerEx(
2516 732 : !oSrcSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oSrcSRS)
2517 : : nullptr,
2518 732 : !oDstSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oDstSRS)
2519 : : nullptr,
2520 732 : aosOptions.List());
2521 :
2522 732 : if (pszCO)
2523 : {
2524 7 : psInfo->bHasCustomTransformationPipeline = true;
2525 : }
2526 :
2527 732 : if (psInfo->pReprojectArg == nullptr)
2528 : {
2529 1 : GDALDestroyGenImgProjTransformer(psInfo);
2530 1 : return nullptr;
2531 : }
2532 731 : psInfo->pReproject = GDALReprojectionTransform;
2533 :
2534 : /* --------------------------------------------------------------------
2535 : */
2536 : /* Handle optional reprojection approximation transformer. */
2537 : /* --------------------------------------------------------------------
2538 : */
2539 731 : const char *psApproxErrorFwd = CSLFetchNameValue(
2540 : papszOptions, "REPROJECTION_APPROX_ERROR_IN_DST_SRS_UNIT");
2541 731 : const char *psApproxErrorReverse = CSLFetchNameValue(
2542 : papszOptions, "REPROJECTION_APPROX_ERROR_IN_SRC_SRS_UNIT");
2543 731 : if (psApproxErrorFwd && psApproxErrorReverse)
2544 : {
2545 1 : void *pArg = GDALCreateApproxTransformer2(
2546 : psInfo->pReproject, psInfo->pReprojectArg,
2547 : CPLAtof(psApproxErrorFwd), CPLAtof(psApproxErrorReverse));
2548 1 : if (pArg == nullptr)
2549 : {
2550 0 : GDALDestroyGenImgProjTransformer(psInfo);
2551 0 : return nullptr;
2552 : }
2553 1 : psInfo->pReprojectArg = pArg;
2554 1 : psInfo->pReproject = GDALApproxTransform;
2555 1 : GDALApproxTransformerOwnsSubtransformer(psInfo->pReprojectArg,
2556 : TRUE);
2557 : }
2558 : }
2559 :
2560 1397 : return psInfo;
2561 : }
2562 :
2563 : /************************************************************************/
2564 : /* GDALRefreshGenImgProjTransformer() */
2565 : /************************************************************************/
2566 :
2567 968 : void GDALRefreshGenImgProjTransformer(void *hTransformArg)
2568 : {
2569 968 : GDALGenImgProjTransformInfo *psInfo =
2570 : static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
2571 :
2572 1609 : if (psInfo->pReprojectArg &&
2573 641 : psInfo->bCheckWithInvertPROJ != GetCurrentCheckWithInvertPROJ())
2574 : {
2575 66 : psInfo->bCheckWithInvertPROJ = !psInfo->bCheckWithInvertPROJ;
2576 :
2577 : CPLXMLNode *psXML =
2578 66 : GDALSerializeTransformer(psInfo->pReproject, psInfo->pReprojectArg);
2579 66 : GDALDestroyTransformer(psInfo->pReprojectArg);
2580 66 : GDALDeserializeTransformer(psXML, &psInfo->pReproject,
2581 : &psInfo->pReprojectArg);
2582 66 : CPLDestroyXMLNode(psXML);
2583 : }
2584 968 : }
2585 :
2586 : /************************************************************************/
2587 : /* GDALCreateGenImgProjTransformer3() */
2588 : /************************************************************************/
2589 :
2590 : /**
2591 : * Create image to image transformer.
2592 : *
2593 : * This function creates a transformation object that maps from pixel/line
2594 : * coordinates on one image to pixel/line coordinates on another image. The
2595 : * images may potentially be georeferenced in different coordinate systems,
2596 : * and may used GCPs to map between their pixel/line coordinates and
2597 : * georeferenced coordinates (as opposed to the default assumption that their
2598 : * geotransform should be used).
2599 : *
2600 : * This transformer potentially performs three concatenated transformations.
2601 : *
2602 : * The first stage is from source image pixel/line coordinates to source
2603 : * image georeferenced coordinates, and may be done using the geotransform,
2604 : * or if not defined using a polynomial model derived from GCPs. If GCPs
2605 : * are used this stage is accomplished using GDALGCPTransform().
2606 : *
2607 : * The second stage is to change projections from the source coordinate system
2608 : * to the destination coordinate system, assuming they differ. This is
2609 : * accomplished internally using GDALReprojectionTransform().
2610 : *
2611 : * The third stage is converting from destination image georeferenced
2612 : * coordinates to destination image coordinates. This is done using the
2613 : * destination image geotransform, or if not available, using a polynomial
2614 : * model derived from GCPs. If GCPs are used this stage is accomplished using
2615 : * GDALGCPTransform(). This stage is skipped if hDstDS is NULL when the
2616 : * transformation is created.
2617 : *
2618 : * @param pszSrcWKT source WKT (or NULL).
2619 : * @param padfSrcGeoTransform source geotransform (or NULL).
2620 : * @param pszDstWKT destination WKT (or NULL).
2621 : * @param padfDstGeoTransform destination geotransform (or NULL).
2622 : *
2623 : * @return handle suitable for use GDALGenImgProjTransform(), and to be
2624 : * deallocated with GDALDestroyGenImgProjTransformer() or NULL on failure.
2625 : */
2626 :
2627 0 : void *GDALCreateGenImgProjTransformer3(const char *pszSrcWKT,
2628 : const double *padfSrcGeoTransform,
2629 : const char *pszDstWKT,
2630 : const double *padfDstGeoTransform)
2631 :
2632 : {
2633 0 : OGRSpatialReference oSrcSRS;
2634 0 : if (pszSrcWKT)
2635 : {
2636 0 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2637 0 : if (pszSrcWKT[0] != '\0' &&
2638 0 : oSrcSRS.importFromWkt(pszSrcWKT) != OGRERR_NONE)
2639 : {
2640 0 : CPLError(CE_Failure, CPLE_AppDefined,
2641 : "Failed to import coordinate system `%s'.", pszSrcWKT);
2642 0 : return nullptr;
2643 : }
2644 : }
2645 :
2646 0 : OGRSpatialReference oDstSRS;
2647 0 : if (pszDstWKT)
2648 : {
2649 0 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2650 0 : if (pszDstWKT[0] != '\0' &&
2651 0 : oDstSRS.importFromWkt(pszDstWKT) != OGRERR_NONE)
2652 : {
2653 0 : CPLError(CE_Failure, CPLE_AppDefined,
2654 : "Failed to import coordinate system `%s'.", pszDstWKT);
2655 0 : return nullptr;
2656 : }
2657 : }
2658 0 : return GDALCreateGenImgProjTransformer4(
2659 : OGRSpatialReference::ToHandle(&oSrcSRS), padfSrcGeoTransform,
2660 0 : OGRSpatialReference::ToHandle(&oDstSRS), padfDstGeoTransform, nullptr);
2661 : }
2662 :
2663 : /************************************************************************/
2664 : /* GDALCreateGenImgProjTransformer4() */
2665 : /************************************************************************/
2666 :
2667 : /**
2668 : * Create image to image transformer.
2669 : *
2670 : * Similar to GDALCreateGenImgProjTransformer3(), except that it takes
2671 : * OGRSpatialReferenceH objects and options.
2672 : * The options are the ones supported by GDALCreateReprojectionTransformerEx()
2673 : *
2674 : * @since GDAL 3.0
2675 : */
2676 16 : void *GDALCreateGenImgProjTransformer4(OGRSpatialReferenceH hSrcSRS,
2677 : const double *padfSrcGeoTransform,
2678 : OGRSpatialReferenceH hDstSRS,
2679 : const double *padfDstGeoTransform,
2680 : const char *const *papszOptions)
2681 : {
2682 : /* -------------------------------------------------------------------- */
2683 : /* Initialize the transform info. */
2684 : /* -------------------------------------------------------------------- */
2685 : GDALGenImgProjTransformInfo *psInfo =
2686 16 : GDALCreateGenImgProjTransformerInternal();
2687 :
2688 : /* -------------------------------------------------------------------- */
2689 : /* Get forward and inverse geotransform for the source image. */
2690 : /* -------------------------------------------------------------------- */
2691 16 : if (padfSrcGeoTransform)
2692 : {
2693 16 : memcpy(psInfo->adfSrcGeoTransform, padfSrcGeoTransform,
2694 : sizeof(psInfo->adfSrcGeoTransform));
2695 16 : if (!GDALInvGeoTransform(psInfo->adfSrcGeoTransform,
2696 16 : psInfo->adfSrcInvGeoTransform))
2697 : {
2698 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
2699 1 : GDALDestroyGenImgProjTransformer(psInfo);
2700 1 : return nullptr;
2701 : }
2702 : }
2703 : else
2704 : {
2705 0 : psInfo->adfSrcGeoTransform[0] = 0.0;
2706 0 : psInfo->adfSrcGeoTransform[1] = 1.0;
2707 0 : psInfo->adfSrcGeoTransform[2] = 0.0;
2708 0 : psInfo->adfSrcGeoTransform[3] = 0.0;
2709 0 : psInfo->adfSrcGeoTransform[4] = 0.0;
2710 0 : psInfo->adfSrcGeoTransform[5] = 1.0;
2711 0 : memcpy(psInfo->adfSrcInvGeoTransform, psInfo->adfSrcGeoTransform,
2712 : sizeof(double) * 6);
2713 : }
2714 :
2715 : /* -------------------------------------------------------------------- */
2716 : /* Setup reprojection. */
2717 : /* -------------------------------------------------------------------- */
2718 15 : OGRSpatialReference *poSrcSRS = OGRSpatialReference::FromHandle(hSrcSRS);
2719 15 : OGRSpatialReference *poDstSRS = OGRSpatialReference::FromHandle(hDstSRS);
2720 30 : if (!poSrcSRS->IsEmpty() && !poDstSRS->IsEmpty() &&
2721 15 : !poSrcSRS->IsSame(poDstSRS))
2722 : {
2723 4 : psInfo->pReprojectArg =
2724 4 : GDALCreateReprojectionTransformerEx(hSrcSRS, hDstSRS, papszOptions);
2725 4 : if (psInfo->pReprojectArg == nullptr)
2726 : {
2727 1 : GDALDestroyGenImgProjTransformer(psInfo);
2728 1 : return nullptr;
2729 : }
2730 3 : psInfo->pReproject = GDALReprojectionTransform;
2731 : }
2732 :
2733 : /* -------------------------------------------------------------------- */
2734 : /* Get forward and inverse geotransform for destination image. */
2735 : /* If we have no destination matrix use a unit transform. */
2736 : /* -------------------------------------------------------------------- */
2737 14 : if (padfDstGeoTransform)
2738 : {
2739 14 : memcpy(psInfo->adfDstGeoTransform, padfDstGeoTransform,
2740 : sizeof(psInfo->adfDstGeoTransform));
2741 14 : if (!GDALInvGeoTransform(psInfo->adfDstGeoTransform,
2742 14 : psInfo->adfDstInvGeoTransform))
2743 : {
2744 1 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
2745 1 : GDALDestroyGenImgProjTransformer(psInfo);
2746 1 : return nullptr;
2747 : }
2748 : }
2749 : else
2750 : {
2751 0 : psInfo->adfDstGeoTransform[0] = 0.0;
2752 0 : psInfo->adfDstGeoTransform[1] = 1.0;
2753 0 : psInfo->adfDstGeoTransform[2] = 0.0;
2754 0 : psInfo->adfDstGeoTransform[3] = 0.0;
2755 0 : psInfo->adfDstGeoTransform[4] = 0.0;
2756 0 : psInfo->adfDstGeoTransform[5] = 1.0;
2757 0 : memcpy(psInfo->adfDstInvGeoTransform, psInfo->adfDstGeoTransform,
2758 : sizeof(double) * 6);
2759 : }
2760 :
2761 13 : return psInfo;
2762 : }
2763 :
2764 : /************************************************************************/
2765 : /* GDALSetGenImgProjTransformerDstGeoTransform() */
2766 : /************************************************************************/
2767 :
2768 : /**
2769 : * Set GenImgProj output geotransform.
2770 : *
2771 : * Normally the "destination geotransform", or transformation between
2772 : * georeferenced output coordinates and pixel/line coordinates on the
2773 : * destination file is extracted from the destination file by
2774 : * GDALCreateGenImgProjTransformer() and stored in the GenImgProj private
2775 : * info. However, sometimes it is inconvenient to have an output file
2776 : * handle with appropriate geotransform information when creating the
2777 : * transformation. For these cases, this function can be used to apply
2778 : * the destination geotransform.
2779 : *
2780 : * @param hTransformArg the handle to update.
2781 : * @param padfGeoTransform the destination geotransform to apply (six doubles).
2782 : */
2783 :
2784 828 : void GDALSetGenImgProjTransformerDstGeoTransform(void *hTransformArg,
2785 : const double *padfGeoTransform)
2786 :
2787 : {
2788 828 : VALIDATE_POINTER0(hTransformArg,
2789 : "GDALSetGenImgProjTransformerDstGeoTransform");
2790 :
2791 828 : GDALGenImgProjTransformInfo *psInfo =
2792 : static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
2793 :
2794 828 : memcpy(psInfo->adfDstGeoTransform, padfGeoTransform, sizeof(double) * 6);
2795 828 : if (!GDALInvGeoTransform(psInfo->adfDstGeoTransform,
2796 828 : psInfo->adfDstInvGeoTransform))
2797 : {
2798 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
2799 : }
2800 : }
2801 :
2802 : /************************************************************************/
2803 : /* GDALDestroyGenImgProjTransformer() */
2804 : /************************************************************************/
2805 :
2806 : /**
2807 : * GenImgProjTransformer deallocator.
2808 : *
2809 : * This function is used to deallocate the handle created with
2810 : * GDALCreateGenImgProjTransformer().
2811 : *
2812 : * @param hTransformArg the handle to deallocate.
2813 : */
2814 :
2815 1659 : void GDALDestroyGenImgProjTransformer(void *hTransformArg)
2816 :
2817 : {
2818 1659 : if (hTransformArg == nullptr)
2819 0 : return;
2820 :
2821 1659 : GDALGenImgProjTransformInfo *psInfo =
2822 : static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
2823 :
2824 1659 : if (psInfo->pSrcTransformArg != nullptr)
2825 144 : GDALDestroyTransformer(psInfo->pSrcTransformArg);
2826 :
2827 1659 : if (psInfo->pDstTransformArg != nullptr)
2828 7 : GDALDestroyTransformer(psInfo->pDstTransformArg);
2829 :
2830 1659 : if (psInfo->pReprojectArg != nullptr)
2831 842 : GDALDestroyTransformer(psInfo->pReprojectArg);
2832 :
2833 1659 : CPLFree(psInfo);
2834 : }
2835 :
2836 : /************************************************************************/
2837 : /* GDALGenImgProjTransform() */
2838 : /************************************************************************/
2839 :
2840 : /**
2841 : * Perform general image reprojection transformation.
2842 : *
2843 : * Actually performs the transformation setup in
2844 : * GDALCreateGenImgProjTransformer(). This function matches the signature
2845 : * required by the GDALTransformerFunc(), and more details on the arguments
2846 : * can be found in that topic.
2847 : */
2848 :
2849 : #ifdef DEBUG_APPROX_TRANSFORMER
2850 : int countGDALGenImgProjTransform = 0;
2851 : #endif
2852 :
2853 1939960 : int GDALGenImgProjTransform(void *pTransformArgIn, int bDstToSrc,
2854 : int nPointCount, double *padfX, double *padfY,
2855 : double *padfZ, int *panSuccess)
2856 : {
2857 1939960 : GDALGenImgProjTransformInfo *psInfo =
2858 : static_cast<GDALGenImgProjTransformInfo *>(pTransformArgIn);
2859 :
2860 : #ifdef DEBUG_APPROX_TRANSFORMER
2861 : CPLAssert(nPointCount > 0);
2862 : countGDALGenImgProjTransform += nPointCount;
2863 : #endif
2864 :
2865 21590600 : for (int i = 0; i < nPointCount; i++)
2866 : {
2867 19650700 : panSuccess[i] = (padfX[i] != HUGE_VAL && padfY[i] != HUGE_VAL);
2868 : }
2869 :
2870 : /* -------------------------------------------------------------------- */
2871 : /* Convert from src (dst) pixel/line to src (dst) */
2872 : /* georeferenced coordinates. */
2873 : /* -------------------------------------------------------------------- */
2874 1939960 : double *padfGeoTransform = nullptr;
2875 1939960 : void *pTransformArg = nullptr;
2876 1939960 : GDALTransformerFunc pTransformer = nullptr;
2877 1939960 : if (bDstToSrc)
2878 : {
2879 1559140 : padfGeoTransform = psInfo->adfDstGeoTransform;
2880 1559140 : pTransformArg = psInfo->pDstTransformArg;
2881 1559140 : pTransformer = psInfo->pDstTransformer;
2882 : }
2883 : else
2884 : {
2885 380817 : padfGeoTransform = psInfo->adfSrcGeoTransform;
2886 380817 : pTransformArg = psInfo->pSrcTransformArg;
2887 380817 : pTransformer = psInfo->pSrcTransformer;
2888 : }
2889 :
2890 1939960 : int ret = TRUE;
2891 1939960 : if (pTransformArg != nullptr)
2892 : {
2893 41607 : if (!pTransformer(pTransformArg, FALSE, nPointCount, padfX, padfY,
2894 : padfZ, panSuccess))
2895 1772 : ret = FALSE;
2896 : }
2897 : else
2898 : {
2899 21470400 : for (int i = 0; i < nPointCount; i++)
2900 : {
2901 19572000 : if (!panSuccess[i])
2902 1758 : continue;
2903 :
2904 19570300 : const double dfNewX = padfGeoTransform[0] +
2905 19570300 : padfX[i] * padfGeoTransform[1] +
2906 19570300 : padfY[i] * padfGeoTransform[2];
2907 19570300 : const double dfNewY = padfGeoTransform[3] +
2908 19570300 : padfX[i] * padfGeoTransform[4] +
2909 19570300 : padfY[i] * padfGeoTransform[5];
2910 :
2911 19570300 : padfX[i] = dfNewX;
2912 19570300 : padfY[i] = dfNewY;
2913 : }
2914 : }
2915 :
2916 : /* -------------------------------------------------------------------- */
2917 : /* Reproject if needed. */
2918 : /* -------------------------------------------------------------------- */
2919 1939960 : if (psInfo->pReprojectArg)
2920 : {
2921 1609160 : if (!psInfo->pReproject(psInfo->pReprojectArg, bDstToSrc, nPointCount,
2922 : padfX, padfY, padfZ, panSuccess))
2923 18077 : ret = FALSE;
2924 : }
2925 :
2926 : /* -------------------------------------------------------------------- */
2927 : /* Convert dst (src) georef coordinates back to pixel/line. */
2928 : /* -------------------------------------------------------------------- */
2929 1939950 : if (bDstToSrc)
2930 : {
2931 1559140 : padfGeoTransform = psInfo->adfSrcInvGeoTransform;
2932 1559140 : pTransformArg = psInfo->pSrcTransformArg;
2933 1559140 : pTransformer = psInfo->pSrcTransformer;
2934 : }
2935 : else
2936 : {
2937 380817 : padfGeoTransform = psInfo->adfDstInvGeoTransform;
2938 380817 : pTransformArg = psInfo->pDstTransformArg;
2939 380817 : pTransformer = psInfo->pDstTransformer;
2940 : }
2941 :
2942 1939950 : if (pTransformArg != nullptr)
2943 : {
2944 51473 : if (!pTransformer(pTransformArg, TRUE, nPointCount, padfX, padfY, padfZ,
2945 : panSuccess))
2946 1640 : ret = FALSE;
2947 : }
2948 : else
2949 : {
2950 20686300 : for (int i = 0; i < nPointCount; i++)
2951 : {
2952 18797800 : if (!panSuccess[i])
2953 3523210 : continue;
2954 :
2955 15274600 : const double dfNewX = padfGeoTransform[0] +
2956 15274600 : padfX[i] * padfGeoTransform[1] +
2957 15274600 : padfY[i] * padfGeoTransform[2];
2958 15274600 : const double dfNewY = padfGeoTransform[3] +
2959 15274600 : padfX[i] * padfGeoTransform[4] +
2960 15274600 : padfY[i] * padfGeoTransform[5];
2961 :
2962 15274600 : padfX[i] = dfNewX;
2963 15274600 : padfY[i] = dfNewY;
2964 : }
2965 : }
2966 :
2967 1939950 : return ret;
2968 : }
2969 :
2970 : /************************************************************************/
2971 : /* GDALTransformLonLatToDestGenImgProjTransformer() */
2972 : /************************************************************************/
2973 :
2974 2536 : int GDALTransformLonLatToDestGenImgProjTransformer(void *hTransformArg,
2975 : double *pdfX, double *pdfY)
2976 : {
2977 2536 : GDALGenImgProjTransformInfo *psInfo =
2978 : static_cast<GDALGenImgProjTransformInfo *>(hTransformArg);
2979 :
2980 2536 : if (psInfo->pReprojectArg == nullptr ||
2981 1418 : psInfo->pReproject != GDALReprojectionTransform)
2982 1122 : return false;
2983 :
2984 1414 : GDALReprojectionTransformInfo *psReprojInfo =
2985 : static_cast<GDALReprojectionTransformInfo *>(psInfo->pReprojectArg);
2986 2828 : if (psReprojInfo->poForwardTransform == nullptr ||
2987 1414 : psReprojInfo->poForwardTransform->GetSourceCS() == nullptr)
2988 2 : return false;
2989 :
2990 1412 : double z = 0;
2991 1412 : int success = true;
2992 1412 : auto poSourceCRS = psReprojInfo->poForwardTransform->GetSourceCS();
2993 2504 : if (poSourceCRS->IsGeographic() &&
2994 1092 : std::fabs(poSourceCRS->GetAngularUnits() -
2995 1092 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-9)
2996 : {
2997 : // Optimization to avoid creating a OGRCoordinateTransformation
2998 1090 : OGRAxisOrientation eSourceFirstAxisOrient = OAO_Other;
2999 1090 : poSourceCRS->GetAxis(nullptr, 0, &eSourceFirstAxisOrient);
3000 1090 : const auto &mapping = poSourceCRS->GetDataAxisToSRSAxisMapping();
3001 2180 : if ((mapping[0] == 2 && eSourceFirstAxisOrient == OAO_East) ||
3002 1090 : (mapping[0] == 1 && eSourceFirstAxisOrient != OAO_East))
3003 : {
3004 6 : std::swap(*pdfX, *pdfY);
3005 : }
3006 : }
3007 : else
3008 : {
3009 : auto poLongLat =
3010 322 : std::unique_ptr<OGRSpatialReference>(poSourceCRS->CloneGeogCS());
3011 322 : if (poLongLat == nullptr)
3012 0 : return false;
3013 322 : poLongLat->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3014 :
3015 : const bool bCurrentCheckWithInvertProj =
3016 322 : GetCurrentCheckWithInvertPROJ();
3017 322 : if (!bCurrentCheckWithInvertProj)
3018 322 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", "YES");
3019 : auto poCT = std::unique_ptr<OGRCoordinateTransformation>(
3020 322 : OGRCreateCoordinateTransformation(poLongLat.get(), poSourceCRS));
3021 322 : if (!bCurrentCheckWithInvertProj)
3022 322 : CPLSetThreadLocalConfigOption("CHECK_WITH_INVERT_PROJ", nullptr);
3023 322 : if (poCT == nullptr)
3024 0 : return false;
3025 :
3026 322 : poCT->SetEmitErrors(false);
3027 322 : if (!poCT->Transform(1, pdfX, pdfY))
3028 2 : return false;
3029 :
3030 320 : if (!psInfo->pReproject(psInfo->pReprojectArg, false, 1, pdfX, pdfY, &z,
3031 502 : &success) ||
3032 182 : !success)
3033 : {
3034 138 : return false;
3035 : }
3036 : }
3037 :
3038 1272 : double *padfGeoTransform = psInfo->adfDstInvGeoTransform;
3039 1272 : void *pTransformArg = psInfo->pDstTransformArg;
3040 1272 : GDALTransformerFunc pTransformer = psInfo->pDstTransformer;
3041 1272 : if (pTransformArg != nullptr)
3042 : {
3043 4 : if (!pTransformer(pTransformArg, TRUE, 1, pdfX, pdfY, &z, &success) ||
3044 0 : !success)
3045 : {
3046 4 : return false;
3047 : }
3048 : }
3049 : else
3050 : {
3051 1268 : const double dfNewX = padfGeoTransform[0] +
3052 1268 : pdfX[0] * padfGeoTransform[1] +
3053 1268 : pdfY[0] * padfGeoTransform[2];
3054 1268 : const double dfNewY = padfGeoTransform[3] +
3055 1268 : pdfX[0] * padfGeoTransform[4] +
3056 1268 : pdfY[0] * padfGeoTransform[5];
3057 :
3058 1268 : pdfX[0] = dfNewX;
3059 1268 : pdfY[0] = dfNewY;
3060 : }
3061 :
3062 1268 : return true;
3063 : }
3064 :
3065 : /************************************************************************/
3066 : /* GDALSerializeGenImgProjTransformer() */
3067 : /************************************************************************/
3068 :
3069 77 : static CPLXMLNode *GDALSerializeGenImgProjTransformer(void *pTransformArg)
3070 :
3071 : {
3072 77 : GDALGenImgProjTransformInfo *psInfo =
3073 : static_cast<GDALGenImgProjTransformInfo *>(pTransformArg);
3074 :
3075 : CPLXMLNode *psTree =
3076 77 : CPLCreateXMLNode(nullptr, CXT_Element, "GenImgProjTransformer");
3077 :
3078 77 : char szWork[200] = {};
3079 :
3080 : /* -------------------------------------------------------------------- */
3081 : /* Handle source transformation. */
3082 : /* -------------------------------------------------------------------- */
3083 77 : if (psInfo->pSrcTransformArg != nullptr)
3084 : {
3085 13 : CPLXMLNode *psTransformer = GDALSerializeTransformer(
3086 : psInfo->pSrcTransformer, psInfo->pSrcTransformArg);
3087 13 : if (psTransformer != nullptr)
3088 : {
3089 : CPLXMLNode *psTransformerContainer =
3090 13 : CPLCreateXMLNode(psTree, CXT_Element,
3091 : CPLSPrintf("Src%s", psTransformer->pszValue));
3092 :
3093 13 : CPLAddXMLChild(psTransformerContainer, psTransformer);
3094 : }
3095 : }
3096 :
3097 : /* -------------------------------------------------------------------- */
3098 : /* Handle source geotransforms. */
3099 : /* -------------------------------------------------------------------- */
3100 : else
3101 : {
3102 64 : CPLsnprintf(
3103 : szWork, sizeof(szWork), "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g",
3104 : psInfo->adfSrcGeoTransform[0], psInfo->adfSrcGeoTransform[1],
3105 : psInfo->adfSrcGeoTransform[2], psInfo->adfSrcGeoTransform[3],
3106 : psInfo->adfSrcGeoTransform[4], psInfo->adfSrcGeoTransform[5]);
3107 64 : CPLCreateXMLElementAndValue(psTree, "SrcGeoTransform", szWork);
3108 :
3109 64 : CPLsnprintf(
3110 : szWork, sizeof(szWork), "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g",
3111 : psInfo->adfSrcInvGeoTransform[0], psInfo->adfSrcInvGeoTransform[1],
3112 : psInfo->adfSrcInvGeoTransform[2], psInfo->adfSrcInvGeoTransform[3],
3113 : psInfo->adfSrcInvGeoTransform[4], psInfo->adfSrcInvGeoTransform[5]);
3114 64 : CPLCreateXMLElementAndValue(psTree, "SrcInvGeoTransform", szWork);
3115 : }
3116 :
3117 : /* -------------------------------------------------------------------- */
3118 : /* Handle dest transformation. */
3119 : /* -------------------------------------------------------------------- */
3120 77 : if (psInfo->pDstTransformArg != nullptr)
3121 : {
3122 0 : CPLXMLNode *psTransformer = GDALSerializeTransformer(
3123 : psInfo->pDstTransformer, psInfo->pDstTransformArg);
3124 0 : if (psTransformer != nullptr)
3125 : {
3126 : CPLXMLNode *psTransformerContainer =
3127 0 : CPLCreateXMLNode(psTree, CXT_Element,
3128 : CPLSPrintf("Dst%s", psTransformer->pszValue));
3129 :
3130 0 : CPLAddXMLChild(psTransformerContainer, psTransformer);
3131 : }
3132 : }
3133 :
3134 : /* -------------------------------------------------------------------- */
3135 : /* Handle destination geotransforms. */
3136 : /* -------------------------------------------------------------------- */
3137 : else
3138 : {
3139 77 : CPLsnprintf(
3140 : szWork, sizeof(szWork), "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g",
3141 : psInfo->adfDstGeoTransform[0], psInfo->adfDstGeoTransform[1],
3142 : psInfo->adfDstGeoTransform[2], psInfo->adfDstGeoTransform[3],
3143 : psInfo->adfDstGeoTransform[4], psInfo->adfDstGeoTransform[5]);
3144 77 : CPLCreateXMLElementAndValue(psTree, "DstGeoTransform", szWork);
3145 :
3146 77 : CPLsnprintf(
3147 : szWork, sizeof(szWork), "%.17g,%.17g,%.17g,%.17g,%.17g,%.17g",
3148 : psInfo->adfDstInvGeoTransform[0], psInfo->adfDstInvGeoTransform[1],
3149 : psInfo->adfDstInvGeoTransform[2], psInfo->adfDstInvGeoTransform[3],
3150 : psInfo->adfDstInvGeoTransform[4], psInfo->adfDstInvGeoTransform[5]);
3151 77 : CPLCreateXMLElementAndValue(psTree, "DstInvGeoTransform", szWork);
3152 : }
3153 :
3154 : /* -------------------------------------------------------------------- */
3155 : /* Do we have a reprojection transformer? */
3156 : /* -------------------------------------------------------------------- */
3157 77 : if (psInfo->pReprojectArg != nullptr)
3158 : {
3159 :
3160 : CPLXMLNode *psTransformerContainer =
3161 50 : CPLCreateXMLNode(psTree, CXT_Element, "ReprojectTransformer");
3162 :
3163 : CPLXMLNode *psTransformer =
3164 50 : GDALSerializeTransformer(psInfo->pReproject, psInfo->pReprojectArg);
3165 50 : if (psTransformer != nullptr)
3166 50 : CPLAddXMLChild(psTransformerContainer, psTransformer);
3167 : }
3168 :
3169 77 : return psTree;
3170 : }
3171 :
3172 : /************************************************************************/
3173 : /* GDALDeserializeGeoTransform() */
3174 : /************************************************************************/
3175 :
3176 735 : static void GDALDeserializeGeoTransform(const char *pszGT,
3177 : double adfGeoTransform[6])
3178 : {
3179 735 : CPLsscanf(pszGT, "%lf,%lf,%lf,%lf,%lf,%lf", adfGeoTransform + 0,
3180 : adfGeoTransform + 1, adfGeoTransform + 2, adfGeoTransform + 3,
3181 : adfGeoTransform + 4, adfGeoTransform + 5);
3182 735 : }
3183 :
3184 : /************************************************************************/
3185 : /* GDALDeserializeGenImgProjTransformer() */
3186 : /************************************************************************/
3187 :
3188 198 : void *GDALDeserializeGenImgProjTransformer(CPLXMLNode *psTree)
3189 :
3190 : {
3191 : /* -------------------------------------------------------------------- */
3192 : /* Initialize the transform info. */
3193 : /* -------------------------------------------------------------------- */
3194 : GDALGenImgProjTransformInfo *psInfo =
3195 198 : GDALCreateGenImgProjTransformerInternal();
3196 :
3197 : /* -------------------------------------------------------------------- */
3198 : /* SrcGeotransform */
3199 : /* -------------------------------------------------------------------- */
3200 198 : if (CPLGetXMLNode(psTree, "SrcGeoTransform") != nullptr)
3201 : {
3202 182 : GDALDeserializeGeoTransform(
3203 : CPLGetXMLValue(psTree, "SrcGeoTransform", ""),
3204 182 : psInfo->adfSrcGeoTransform);
3205 :
3206 182 : if (CPLGetXMLNode(psTree, "SrcInvGeoTransform") != nullptr)
3207 : {
3208 182 : GDALDeserializeGeoTransform(
3209 : CPLGetXMLValue(psTree, "SrcInvGeoTransform", ""),
3210 182 : psInfo->adfSrcInvGeoTransform);
3211 : }
3212 : else
3213 : {
3214 0 : if (!GDALInvGeoTransform(psInfo->adfSrcGeoTransform,
3215 0 : psInfo->adfSrcInvGeoTransform))
3216 : {
3217 0 : CPLError(CE_Failure, CPLE_AppDefined,
3218 : "Cannot invert geotransform");
3219 : }
3220 : }
3221 : }
3222 :
3223 : /* -------------------------------------------------------------------- */
3224 : /* Src Transform */
3225 : /* -------------------------------------------------------------------- */
3226 : else
3227 : {
3228 16 : for (CPLXMLNode *psIter = psTree->psChild; psIter != nullptr;
3229 0 : psIter = psIter->psNext)
3230 : {
3231 16 : if (psIter->eType == CXT_Element &&
3232 16 : STARTS_WITH_CI(psIter->pszValue, "Src"))
3233 : {
3234 16 : GDALDeserializeTransformer(psIter->psChild,
3235 : &psInfo->pSrcTransformer,
3236 : &psInfo->pSrcTransformArg);
3237 16 : break;
3238 : }
3239 : }
3240 : }
3241 :
3242 : /* -------------------------------------------------------------------- */
3243 : /* DstGeotransform */
3244 : /* -------------------------------------------------------------------- */
3245 198 : if (CPLGetXMLNode(psTree, "DstGeoTransform") != nullptr)
3246 : {
3247 198 : GDALDeserializeGeoTransform(
3248 : CPLGetXMLValue(psTree, "DstGeoTransform", ""),
3249 198 : psInfo->adfDstGeoTransform);
3250 :
3251 198 : if (CPLGetXMLNode(psTree, "DstInvGeoTransform") != nullptr)
3252 : {
3253 173 : GDALDeserializeGeoTransform(
3254 : CPLGetXMLValue(psTree, "DstInvGeoTransform", ""),
3255 173 : psInfo->adfDstInvGeoTransform);
3256 : }
3257 : else
3258 : {
3259 25 : if (!GDALInvGeoTransform(psInfo->adfDstGeoTransform,
3260 25 : psInfo->adfDstInvGeoTransform))
3261 : {
3262 0 : CPLError(CE_Failure, CPLE_AppDefined,
3263 : "Cannot invert geotransform");
3264 : }
3265 : }
3266 : }
3267 :
3268 : /* -------------------------------------------------------------------- */
3269 : /* Dst Transform */
3270 : /* -------------------------------------------------------------------- */
3271 : else
3272 : {
3273 0 : for (CPLXMLNode *psIter = psTree->psChild; psIter != nullptr;
3274 0 : psIter = psIter->psNext)
3275 : {
3276 0 : if (psIter->eType == CXT_Element &&
3277 0 : STARTS_WITH_CI(psIter->pszValue, "Dst"))
3278 : {
3279 0 : GDALDeserializeTransformer(psIter->psChild,
3280 : &psInfo->pDstTransformer,
3281 : &psInfo->pDstTransformArg);
3282 0 : break;
3283 : }
3284 : }
3285 : }
3286 :
3287 : /* -------------------------------------------------------------------- */
3288 : /* Reproject transformer */
3289 : /* -------------------------------------------------------------------- */
3290 198 : CPLXMLNode *psSubtree = CPLGetXMLNode(psTree, "ReprojectTransformer");
3291 198 : if (psSubtree != nullptr && psSubtree->psChild != nullptr)
3292 : {
3293 97 : GDALDeserializeTransformer(psSubtree->psChild, &psInfo->pReproject,
3294 : &psInfo->pReprojectArg);
3295 : }
3296 :
3297 198 : return psInfo;
3298 : }
3299 :
3300 : /************************************************************************/
3301 : /* GDALCreateReprojectionTransformer() */
3302 : /************************************************************************/
3303 :
3304 : /**
3305 : * Create reprojection transformer.
3306 : *
3307 : * Creates a callback data structure suitable for use with
3308 : * GDALReprojectionTransformation() to represent a transformation from
3309 : * one geographic or projected coordinate system to another. On input
3310 : * the coordinate systems are described in OpenGIS WKT format.
3311 : *
3312 : * Internally the OGRCoordinateTransformation object is used to implement
3313 : * the reprojection.
3314 : *
3315 : * @param pszSrcWKT the coordinate system for the source coordinate system.
3316 : * @param pszDstWKT the coordinate system for the destination coordinate
3317 : * system.
3318 : *
3319 : * @return Handle for use with GDALReprojectionTransform(), or NULL if the
3320 : * system fails to initialize the reprojection.
3321 : **/
3322 :
3323 0 : void *GDALCreateReprojectionTransformer(const char *pszSrcWKT,
3324 : const char *pszDstWKT)
3325 :
3326 : {
3327 : /* -------------------------------------------------------------------- */
3328 : /* Ingest the SRS definitions. */
3329 : /* -------------------------------------------------------------------- */
3330 0 : OGRSpatialReference oSrcSRS;
3331 0 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3332 0 : if (oSrcSRS.importFromWkt(pszSrcWKT) != OGRERR_NONE)
3333 : {
3334 0 : CPLError(CE_Failure, CPLE_AppDefined,
3335 : "Failed to import coordinate system `%s'.", pszSrcWKT);
3336 0 : return nullptr;
3337 : }
3338 :
3339 0 : OGRSpatialReference oDstSRS;
3340 0 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3341 0 : if (oDstSRS.importFromWkt(pszDstWKT) != OGRERR_NONE)
3342 : {
3343 0 : CPLError(CE_Failure, CPLE_AppDefined,
3344 : "Failed to import coordinate system `%s'.", pszSrcWKT);
3345 0 : return nullptr;
3346 : }
3347 :
3348 0 : return GDALCreateReprojectionTransformerEx(
3349 : OGRSpatialReference::ToHandle(&oSrcSRS),
3350 0 : OGRSpatialReference::ToHandle(&oDstSRS), nullptr);
3351 : }
3352 :
3353 : /************************************************************************/
3354 : /* GDALCreateReprojectionTransformerEx() */
3355 : /************************************************************************/
3356 :
3357 : /**
3358 : * Create reprojection transformer.
3359 : *
3360 : * Creates a callback data structure suitable for use with
3361 : * GDALReprojectionTransformation() to represent a transformation from
3362 : * one geographic or projected coordinate system to another.
3363 : *
3364 : * Internally the OGRCoordinateTransformation object is used to implement
3365 : * the reprojection.
3366 : *
3367 : * @param hSrcSRS the coordinate system for the source coordinate system.
3368 : * @param hDstSRS the coordinate system for the destination coordinate
3369 : * system.
3370 : * @param papszOptions NULL-terminated list of options, or NULL. Currently
3371 : * supported options are:
3372 : * <ul>
3373 : * <li>AREA_OF_INTEREST=west_long,south_lat,east_long,north_lat: Values in
3374 : * degrees. longitudes in [-180,180], latitudes in [-90,90].</li>
3375 : * <li>COORDINATE_OPERATION=string: PROJ or WKT string representing a
3376 : * coordinate operation, overriding the default computed transformation.</li>
3377 : * <li>COORDINATE_EPOCH=decimal_year: Coordinate epoch, expressed as a
3378 : * decimal year. Useful for time-dependent coordinate operations.</li>
3379 : * <li> SRC_COORDINATE_EPOCH: (GDAL >= 3.4) Coordinate epoch of source CRS,
3380 : * expressed as a decimal year. Useful for time-dependent coordinate
3381 : *operations.</li> <li> DST_COORDINATE_EPOCH: (GDAL >= 3.4) Coordinate epoch
3382 : *of target CRS, expressed as a decimal year. Useful for time-dependent
3383 : *coordinate operations.</li>
3384 : * </ul>
3385 : *
3386 : * @return Handle for use with GDALReprojectionTransform(), or NULL if the
3387 : * system fails to initialize the reprojection.
3388 : *
3389 : * @since GDAL 3.0
3390 : **/
3391 :
3392 910 : void *GDALCreateReprojectionTransformerEx(OGRSpatialReferenceH hSrcSRS,
3393 : OGRSpatialReferenceH hDstSRS,
3394 : const char *const *papszOptions)
3395 : {
3396 910 : OGRSpatialReference *poSrcSRS = OGRSpatialReference::FromHandle(hSrcSRS);
3397 910 : OGRSpatialReference *poDstSRS = OGRSpatialReference::FromHandle(hDstSRS);
3398 :
3399 : /* -------------------------------------------------------------------- */
3400 : /* Build the forward coordinate transformation. */
3401 : /* -------------------------------------------------------------------- */
3402 910 : double dfWestLongitudeDeg = 0.0;
3403 910 : double dfSouthLatitudeDeg = 0.0;
3404 910 : double dfEastLongitudeDeg = 0.0;
3405 910 : double dfNorthLatitudeDeg = 0.0;
3406 910 : const char *pszBBOX = CSLFetchNameValue(papszOptions, "AREA_OF_INTEREST");
3407 910 : if (pszBBOX)
3408 : {
3409 844 : char **papszTokens = CSLTokenizeString2(pszBBOX, ",", 0);
3410 844 : if (CSLCount(papszTokens) == 4)
3411 : {
3412 844 : dfWestLongitudeDeg = CPLAtof(papszTokens[0]);
3413 844 : dfSouthLatitudeDeg = CPLAtof(papszTokens[1]);
3414 844 : dfEastLongitudeDeg = CPLAtof(papszTokens[2]);
3415 844 : dfNorthLatitudeDeg = CPLAtof(papszTokens[3]);
3416 : }
3417 844 : CSLDestroy(papszTokens);
3418 : }
3419 910 : const char *pszCO = CSLFetchNameValue(papszOptions, "COORDINATE_OPERATION");
3420 :
3421 1820 : OGRCoordinateTransformationOptions optionsFwd;
3422 910 : if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
3423 : dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))
3424 : {
3425 844 : optionsFwd.SetAreaOfInterest(dfWestLongitudeDeg, dfSouthLatitudeDeg,
3426 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
3427 : }
3428 910 : if (pszCO)
3429 : {
3430 7 : optionsFwd.SetCoordinateOperation(pszCO, false);
3431 : }
3432 :
3433 910 : const char *pszCENTER_LONG = CSLFetchNameValue(papszOptions, "CENTER_LONG");
3434 910 : if (pszCENTER_LONG)
3435 : {
3436 605 : optionsFwd.SetSourceCenterLong(CPLAtof(pszCENTER_LONG));
3437 : }
3438 :
3439 : OGRCoordinateTransformation *poForwardTransform =
3440 910 : OGRCreateCoordinateTransformation(poSrcSRS, poDstSRS, optionsFwd);
3441 :
3442 910 : if (poForwardTransform == nullptr)
3443 : // OGRCreateCoordinateTransformation() will report errors on its own.
3444 2 : return nullptr;
3445 :
3446 908 : poForwardTransform->SetEmitErrors(false);
3447 :
3448 : /* -------------------------------------------------------------------- */
3449 : /* Create a structure to hold the transform info, and also */
3450 : /* build reverse transform. We assume that if the forward */
3451 : /* transform can be created, then so can the reverse one. */
3452 : /* -------------------------------------------------------------------- */
3453 908 : GDALReprojectionTransformInfo *psInfo = new GDALReprojectionTransformInfo();
3454 :
3455 908 : psInfo->papszOptions = CSLDuplicate(papszOptions);
3456 908 : psInfo->poForwardTransform = poForwardTransform;
3457 908 : psInfo->dfTime = CPLAtof(CSLFetchNameValueDef(
3458 : papszOptions, "COORDINATE_EPOCH",
3459 : CSLFetchNameValueDef(
3460 : papszOptions, "DST_COORDINATE_EPOCH",
3461 : CSLFetchNameValueDef(papszOptions, "SRC_COORDINATE_EPOCH", "0"))));
3462 908 : psInfo->poReverseTransform = poForwardTransform->GetInverse();
3463 :
3464 908 : if (psInfo->poReverseTransform)
3465 908 : psInfo->poReverseTransform->SetEmitErrors(false);
3466 :
3467 908 : memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
3468 : strlen(GDAL_GTI2_SIGNATURE));
3469 908 : psInfo->sTI.pszClassName = "GDALReprojectionTransformer";
3470 908 : psInfo->sTI.pfnTransform = GDALReprojectionTransform;
3471 908 : psInfo->sTI.pfnCleanup = GDALDestroyReprojectionTransformer;
3472 908 : psInfo->sTI.pfnSerialize = GDALSerializeReprojectionTransformer;
3473 :
3474 908 : return psInfo;
3475 : }
3476 :
3477 : /************************************************************************/
3478 : /* GDALDestroyReprojectionTransformer() */
3479 : /************************************************************************/
3480 :
3481 : /**
3482 : * Destroy reprojection transformation.
3483 : *
3484 : * @param pTransformArg the transformation handle returned by
3485 : * GDALCreateReprojectionTransformer().
3486 : */
3487 :
3488 908 : void GDALDestroyReprojectionTransformer(void *pTransformArg)
3489 :
3490 : {
3491 908 : if (pTransformArg == nullptr)
3492 0 : return;
3493 :
3494 908 : GDALReprojectionTransformInfo *psInfo =
3495 : static_cast<GDALReprojectionTransformInfo *>(pTransformArg);
3496 :
3497 908 : if (psInfo->poForwardTransform)
3498 908 : OGRCoordinateTransformation::DestroyCT(psInfo->poForwardTransform);
3499 :
3500 908 : if (psInfo->poReverseTransform)
3501 908 : OGRCoordinateTransformation::DestroyCT(psInfo->poReverseTransform);
3502 :
3503 908 : CSLDestroy(psInfo->papszOptions);
3504 :
3505 908 : delete psInfo;
3506 : }
3507 :
3508 : /************************************************************************/
3509 : /* GDALReprojectionTransform() */
3510 : /************************************************************************/
3511 :
3512 : /**
3513 : * Perform reprojection transformation.
3514 : *
3515 : * Actually performs the reprojection transformation described in
3516 : * GDALCreateReprojectionTransformer(). This function matches the
3517 : * GDALTransformerFunc() signature. Details of the arguments are described
3518 : * there.
3519 : */
3520 :
3521 1609480 : int GDALReprojectionTransform(void *pTransformArg, int bDstToSrc,
3522 : int nPointCount, double *padfX, double *padfY,
3523 : double *padfZ, int *panSuccess)
3524 :
3525 : {
3526 1609480 : GDALReprojectionTransformInfo *psInfo =
3527 : static_cast<GDALReprojectionTransformInfo *>(pTransformArg);
3528 : int bSuccess;
3529 :
3530 1609480 : std::vector<double> adfTime;
3531 1609480 : double *padfT = nullptr;
3532 1609480 : if (psInfo->dfTime != 0.0 && nPointCount > 0)
3533 : {
3534 1 : adfTime.resize(nPointCount, psInfo->dfTime);
3535 1 : padfT = &adfTime[0];
3536 : }
3537 :
3538 1609480 : if (bDstToSrc)
3539 : {
3540 1371380 : if (psInfo->poReverseTransform == nullptr)
3541 : {
3542 0 : CPLError(
3543 : CE_Failure, CPLE_AppDefined,
3544 : "Inverse coordinate transformation cannot be instantiated");
3545 0 : if (panSuccess)
3546 : {
3547 0 : for (int i = 0; i < nPointCount; i++)
3548 0 : panSuccess[i] = FALSE;
3549 : }
3550 0 : bSuccess = false;
3551 : }
3552 : else
3553 : {
3554 1371380 : bSuccess = psInfo->poReverseTransform->Transform(
3555 1371380 : nPointCount, padfX, padfY, padfZ, padfT, panSuccess);
3556 : }
3557 : }
3558 : else
3559 238099 : bSuccess = psInfo->poForwardTransform->Transform(
3560 238099 : nPointCount, padfX, padfY, padfZ, padfT, panSuccess);
3561 :
3562 3218960 : return bSuccess;
3563 : }
3564 :
3565 : /************************************************************************/
3566 : /* GDALSerializeReprojectionTransformer() */
3567 : /************************************************************************/
3568 :
3569 127 : static CPLXMLNode *GDALSerializeReprojectionTransformer(void *pTransformArg)
3570 :
3571 : {
3572 : CPLXMLNode *psTree;
3573 127 : GDALReprojectionTransformInfo *psInfo =
3574 : static_cast<GDALReprojectionTransformInfo *>(pTransformArg);
3575 :
3576 127 : psTree = CPLCreateXMLNode(nullptr, CXT_Element, "ReprojectionTransformer");
3577 :
3578 : /* -------------------------------------------------------------------- */
3579 : /* Handle SourceCS. */
3580 : /* -------------------------------------------------------------------- */
3581 254 : const auto ExportToWkt = [](const OGRSpatialReference *poSRS)
3582 : {
3583 : // Try first in WKT1 for backward compat
3584 : {
3585 254 : char *pszWKT = nullptr;
3586 254 : const char *const apszOptions[] = {"FORMAT=WKT1", nullptr};
3587 254 : CPLErrorHandlerPusher oHandler(CPLQuietErrorHandler);
3588 254 : CPLErrorStateBackuper oBackuper;
3589 254 : if (poSRS->exportToWkt(&pszWKT, apszOptions) == OGRERR_NONE)
3590 : {
3591 506 : std::string osRet(pszWKT);
3592 253 : CPLFree(pszWKT);
3593 253 : return osRet;
3594 : }
3595 1 : CPLFree(pszWKT);
3596 : }
3597 :
3598 1 : char *pszWKT = nullptr;
3599 1 : const char *const apszOptions[] = {"FORMAT=WKT2_2019", nullptr};
3600 1 : if (poSRS->exportToWkt(&pszWKT, apszOptions) == OGRERR_NONE)
3601 : {
3602 2 : std::string osRet(pszWKT);
3603 1 : CPLFree(pszWKT);
3604 1 : return osRet;
3605 : }
3606 0 : CPLFree(pszWKT);
3607 0 : return std::string();
3608 : };
3609 :
3610 127 : auto poSRS = psInfo->poForwardTransform->GetSourceCS();
3611 127 : if (poSRS)
3612 : {
3613 254 : const auto osWKT = ExportToWkt(poSRS);
3614 127 : CPLCreateXMLElementAndValue(psTree, "SourceSRS", osWKT.c_str());
3615 : }
3616 :
3617 : /* -------------------------------------------------------------------- */
3618 : /* Handle DestinationCS. */
3619 : /* -------------------------------------------------------------------- */
3620 127 : poSRS = psInfo->poForwardTransform->GetTargetCS();
3621 127 : if (poSRS)
3622 : {
3623 254 : const auto osWKT = ExportToWkt(poSRS);
3624 127 : CPLCreateXMLElementAndValue(psTree, "TargetSRS", osWKT.c_str());
3625 : }
3626 :
3627 : /* -------------------------------------------------------------------- */
3628 : /* Serialize options. */
3629 : /* -------------------------------------------------------------------- */
3630 127 : if (psInfo->papszOptions)
3631 : {
3632 : CPLXMLNode *psOptions =
3633 114 : CPLCreateXMLNode(psTree, CXT_Element, "Options");
3634 276 : for (auto iter = psInfo->papszOptions; *iter != nullptr; ++iter)
3635 : {
3636 162 : char *pszKey = nullptr;
3637 162 : const char *pszValue = CPLParseNameValue(*iter, &pszKey);
3638 162 : if (pszKey && pszValue)
3639 : {
3640 : auto elt =
3641 162 : CPLCreateXMLElementAndValue(psOptions, "Option", pszValue);
3642 162 : CPLAddXMLAttributeAndValue(elt, "key", pszKey);
3643 : }
3644 162 : CPLFree(pszKey);
3645 : }
3646 : }
3647 :
3648 127 : return psTree;
3649 : }
3650 :
3651 : /************************************************************************/
3652 : /* GDALDeserializeReprojectionTransformer() */
3653 : /************************************************************************/
3654 :
3655 174 : static void *GDALDeserializeReprojectionTransformer(CPLXMLNode *psTree)
3656 :
3657 : {
3658 174 : const char *pszSourceSRS = CPLGetXMLValue(psTree, "SourceSRS", nullptr);
3659 174 : const char *pszTargetSRS = CPLGetXMLValue(psTree, "TargetSRS", nullptr);
3660 174 : char *pszSourceWKT = nullptr, *pszTargetWKT = nullptr;
3661 174 : void *pResult = nullptr;
3662 :
3663 348 : OGRSpatialReference oSrcSRS;
3664 348 : OGRSpatialReference oDstSRS;
3665 :
3666 174 : oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3667 174 : oDstSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3668 174 : if (pszSourceSRS != nullptr)
3669 : {
3670 174 : oSrcSRS.SetFromUserInput(pszSourceSRS);
3671 : }
3672 :
3673 174 : if (pszTargetSRS != nullptr)
3674 : {
3675 174 : oDstSRS.SetFromUserInput(pszTargetSRS);
3676 : }
3677 :
3678 174 : CPLStringList aosList;
3679 174 : const CPLXMLNode *psOptions = CPLGetXMLNode(psTree, "Options");
3680 174 : if (psOptions)
3681 : {
3682 344 : for (auto iter = psOptions->psChild; iter; iter = iter->psNext)
3683 : {
3684 206 : if (iter->eType == CXT_Element &&
3685 206 : strcmp(iter->pszValue, "Option") == 0)
3686 : {
3687 206 : const char *pszKey = CPLGetXMLValue(iter, "key", nullptr);
3688 206 : const char *pszValue = CPLGetXMLValue(iter, nullptr, nullptr);
3689 206 : if (pszKey && pszValue)
3690 : {
3691 206 : aosList.SetNameValue(pszKey, pszValue);
3692 : }
3693 : }
3694 : }
3695 : }
3696 :
3697 174 : pResult = GDALCreateReprojectionTransformerEx(
3698 174 : !oSrcSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oSrcSRS) : nullptr,
3699 174 : !oDstSRS.IsEmpty() ? OGRSpatialReference::ToHandle(&oDstSRS) : nullptr,
3700 174 : aosList.List());
3701 :
3702 174 : CPLFree(pszSourceWKT);
3703 174 : CPLFree(pszTargetWKT);
3704 :
3705 348 : return pResult;
3706 : }
3707 :
3708 : /************************************************************************/
3709 : /* ==================================================================== */
3710 : /* Approximate transformer. */
3711 : /* ==================================================================== */
3712 : /************************************************************************/
3713 :
3714 : typedef struct
3715 : {
3716 : GDALTransformerInfo sTI;
3717 :
3718 : GDALTransformerFunc pfnBaseTransformer;
3719 : void *pBaseCBData;
3720 : double dfMaxErrorForward;
3721 : double dfMaxErrorReverse;
3722 :
3723 : int bOwnSubtransformer;
3724 : } ApproxTransformInfo;
3725 :
3726 : /************************************************************************/
3727 : /* GDALCreateSimilarApproxTransformer() */
3728 : /************************************************************************/
3729 :
3730 17 : static void *GDALCreateSimilarApproxTransformer(void *hTransformArg,
3731 : double dfSrcRatioX,
3732 : double dfSrcRatioY)
3733 : {
3734 17 : VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarApproxTransformer",
3735 : nullptr);
3736 :
3737 17 : ApproxTransformInfo *psInfo =
3738 : static_cast<ApproxTransformInfo *>(hTransformArg);
3739 :
3740 : ApproxTransformInfo *psClonedInfo = static_cast<ApproxTransformInfo *>(
3741 17 : CPLMalloc(sizeof(ApproxTransformInfo)));
3742 :
3743 17 : memcpy(psClonedInfo, psInfo, sizeof(ApproxTransformInfo));
3744 17 : if (psClonedInfo->pBaseCBData)
3745 : {
3746 17 : psClonedInfo->pBaseCBData = GDALCreateSimilarTransformer(
3747 : psInfo->pBaseCBData, dfSrcRatioX, dfSrcRatioY);
3748 17 : if (psClonedInfo->pBaseCBData == nullptr)
3749 : {
3750 0 : CPLFree(psClonedInfo);
3751 0 : return nullptr;
3752 : }
3753 : }
3754 17 : psClonedInfo->bOwnSubtransformer = TRUE;
3755 :
3756 17 : return psClonedInfo;
3757 : }
3758 :
3759 : /************************************************************************/
3760 : /* GDALSerializeApproxTransformer() */
3761 : /************************************************************************/
3762 :
3763 57 : static CPLXMLNode *GDALSerializeApproxTransformer(void *pTransformArg)
3764 :
3765 : {
3766 : CPLXMLNode *psTree;
3767 57 : ApproxTransformInfo *psInfo =
3768 : static_cast<ApproxTransformInfo *>(pTransformArg);
3769 :
3770 57 : psTree = CPLCreateXMLNode(nullptr, CXT_Element, "ApproxTransformer");
3771 :
3772 : /* -------------------------------------------------------------------- */
3773 : /* Attach max error. */
3774 : /* -------------------------------------------------------------------- */
3775 57 : if (psInfo->dfMaxErrorForward == psInfo->dfMaxErrorReverse)
3776 : {
3777 55 : CPLCreateXMLElementAndValue(
3778 : psTree, "MaxError",
3779 110 : CPLString().Printf("%g", psInfo->dfMaxErrorForward));
3780 : }
3781 : else
3782 : {
3783 2 : CPLCreateXMLElementAndValue(
3784 : psTree, "MaxErrorForward",
3785 4 : CPLString().Printf("%g", psInfo->dfMaxErrorForward));
3786 2 : CPLCreateXMLElementAndValue(
3787 : psTree, "MaxErrorReverse",
3788 4 : CPLString().Printf("%g", psInfo->dfMaxErrorReverse));
3789 : }
3790 :
3791 : /* -------------------------------------------------------------------- */
3792 : /* Capture underlying transformer. */
3793 : /* -------------------------------------------------------------------- */
3794 : CPLXMLNode *psTransformerContainer =
3795 57 : CPLCreateXMLNode(psTree, CXT_Element, "BaseTransformer");
3796 :
3797 57 : CPLXMLNode *psTransformer = GDALSerializeTransformer(
3798 : psInfo->pfnBaseTransformer, psInfo->pBaseCBData);
3799 57 : if (psTransformer != nullptr)
3800 57 : CPLAddXMLChild(psTransformerContainer, psTransformer);
3801 :
3802 57 : return psTree;
3803 : }
3804 :
3805 : /************************************************************************/
3806 : /* GDALCreateApproxTransformer() */
3807 : /************************************************************************/
3808 :
3809 : /**
3810 : * Create an approximating transformer.
3811 : *
3812 : * This function creates a context for an approximated transformer. Basically
3813 : * a high precision transformer is supplied as input and internally linear
3814 : * approximations are computed to generate results to within a defined
3815 : * precision.
3816 : *
3817 : * The approximation is actually done at the point where GDALApproxTransform()
3818 : * calls are made, and depend on the assumption that they are roughly linear.
3819 : * The first and last point passed in must be the extreme values and the
3820 : * intermediate values should describe a curve between the end points. The
3821 : * approximator transforms and centers using the approximate transformer, and
3822 : * then compares the true middle transformed value to a linear approximation
3823 : * based on the end points. If the error is within the supplied threshold then
3824 : * the end points are used to linearly approximate all the values otherwise the
3825 : * input points are split into two smaller sets, and the function is recursively
3826 : * called until a sufficiently small set of points is found that the linear
3827 : * approximation is OK, or that all the points are exactly computed.
3828 : *
3829 : * This function is very suitable for approximating transformation results
3830 : * from output pixel/line space to input coordinates for warpers that operate
3831 : * on one input scanline at a time. Care should be taken using it in other
3832 : * circumstances as little internal validation is done in order to keep things
3833 : * fast.
3834 : *
3835 : * @param pfnBaseTransformer the high precision transformer which should be
3836 : * approximated.
3837 : * @param pBaseTransformArg the callback argument for the high precision
3838 : * transformer.
3839 : * @param dfMaxError the maximum cartesian error in the "output" space that
3840 : * is to be accepted in the linear approximation, evaluated as a Manhattan
3841 : * distance.
3842 : *
3843 : * @return callback pointer suitable for use with GDALApproxTransform(). It
3844 : * should be deallocated with GDALDestroyApproxTransformer().
3845 : */
3846 :
3847 914 : void *GDALCreateApproxTransformer(GDALTransformerFunc pfnBaseTransformer,
3848 : void *pBaseTransformArg, double dfMaxError)
3849 :
3850 : {
3851 914 : return GDALCreateApproxTransformer2(pfnBaseTransformer, pBaseTransformArg,
3852 914 : dfMaxError, dfMaxError);
3853 : }
3854 :
3855 : static void *
3856 1066 : GDALCreateApproxTransformer2(GDALTransformerFunc pfnBaseTransformer,
3857 : void *pBaseTransformArg, double dfMaxErrorForward,
3858 : double dfMaxErrorReverse)
3859 :
3860 : {
3861 : ApproxTransformInfo *psATInfo = static_cast<ApproxTransformInfo *>(
3862 1066 : CPLMalloc(sizeof(ApproxTransformInfo)));
3863 1066 : psATInfo->pfnBaseTransformer = pfnBaseTransformer;
3864 1066 : psATInfo->pBaseCBData = pBaseTransformArg;
3865 1066 : psATInfo->dfMaxErrorForward = dfMaxErrorForward;
3866 1066 : psATInfo->dfMaxErrorReverse = dfMaxErrorReverse;
3867 1066 : psATInfo->bOwnSubtransformer = FALSE;
3868 :
3869 1066 : memcpy(psATInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
3870 : strlen(GDAL_GTI2_SIGNATURE));
3871 1066 : psATInfo->sTI.pszClassName = GDAL_APPROX_TRANSFORMER_CLASS_NAME;
3872 1066 : psATInfo->sTI.pfnTransform = GDALApproxTransform;
3873 1066 : psATInfo->sTI.pfnCleanup = GDALDestroyApproxTransformer;
3874 1066 : psATInfo->sTI.pfnSerialize = GDALSerializeApproxTransformer;
3875 1066 : psATInfo->sTI.pfnCreateSimilar = GDALCreateSimilarApproxTransformer;
3876 :
3877 1066 : return psATInfo;
3878 : }
3879 :
3880 : /************************************************************************/
3881 : /* GDALApproxTransformerOwnsSubtransformer() */
3882 : /************************************************************************/
3883 :
3884 : /** Set bOwnSubtransformer flag */
3885 1063 : void GDALApproxTransformerOwnsSubtransformer(void *pCBData, int bOwnFlag)
3886 :
3887 : {
3888 1063 : ApproxTransformInfo *psATInfo = static_cast<ApproxTransformInfo *>(pCBData);
3889 :
3890 1063 : psATInfo->bOwnSubtransformer = bOwnFlag;
3891 1063 : }
3892 :
3893 : /************************************************************************/
3894 : /* GDALDestroyApproxTransformer() */
3895 : /************************************************************************/
3896 :
3897 : /**
3898 : * Cleanup approximate transformer.
3899 : *
3900 : * Deallocates the resources allocated by GDALCreateApproxTransformer().
3901 : *
3902 : * @param pCBData callback data originally returned by
3903 : * GDALCreateApproxTransformer().
3904 : */
3905 :
3906 1083 : void GDALDestroyApproxTransformer(void *pCBData)
3907 :
3908 : {
3909 1083 : if (pCBData == nullptr)
3910 0 : return;
3911 :
3912 1083 : ApproxTransformInfo *psATInfo = static_cast<ApproxTransformInfo *>(pCBData);
3913 :
3914 1083 : if (psATInfo->bOwnSubtransformer)
3915 1080 : GDALDestroyTransformer(psATInfo->pBaseCBData);
3916 :
3917 1083 : CPLFree(pCBData);
3918 : }
3919 :
3920 : /************************************************************************/
3921 : /* GDALRefreshApproxTransformer() */
3922 : /************************************************************************/
3923 :
3924 44 : void GDALRefreshApproxTransformer(void *hTransformArg)
3925 : {
3926 44 : ApproxTransformInfo *psInfo =
3927 : static_cast<ApproxTransformInfo *>(hTransformArg);
3928 :
3929 44 : if (GDALIsTransformer(psInfo->pBaseCBData,
3930 : GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
3931 : {
3932 44 : GDALRefreshGenImgProjTransformer(psInfo->pBaseCBData);
3933 : }
3934 44 : }
3935 :
3936 : /************************************************************************/
3937 : /* GDALApproxTransformInternal() */
3938 : /************************************************************************/
3939 :
3940 1111520 : static int GDALApproxTransformInternal(void *pCBData, int bDstToSrc,
3941 : int nPoints, double *x, double *y,
3942 : double *z, int *panSuccess,
3943 : // SME = Start, Middle, End.
3944 : const double xSMETransformed[3],
3945 : const double ySMETransformed[3],
3946 : const double zSMETransformed[3])
3947 : {
3948 1111520 : ApproxTransformInfo *psATInfo = static_cast<ApproxTransformInfo *>(pCBData);
3949 1111520 : const int nMiddle = (nPoints - 1) / 2;
3950 :
3951 : #ifdef notdef_sanify_check
3952 : {
3953 : double x2[3] = {x[0], x[nMiddle], x[nPoints - 1]};
3954 : double y2[3] = {y[0], y[nMiddle], y[nPoints - 1]};
3955 : double z2[3] = {z[0], z[nMiddle], z[nPoints - 1]};
3956 : int anSuccess2[3] = {};
3957 :
3958 : const int bSuccess = psATInfo->pfnBaseTransformer(
3959 : psATInfo->pBaseCBData, bDstToSrc, 3, x2, y2, z2, anSuccess2);
3960 : CPLAssert(bSuccess);
3961 : CPLAssert(anSuccess2[0]);
3962 : CPLAssert(anSuccess2[1]);
3963 : CPLAssert(anSuccess2[2]);
3964 : CPLAssert(x2[0] == xSMETransformed[0]);
3965 : CPLAssert(y2[0] == ySMETransformed[0]);
3966 : CPLAssert(z2[0] == zSMETransformed[0]);
3967 : CPLAssert(x2[1] == xSMETransformed[1]);
3968 : CPLAssert(y2[1] == ySMETransformed[1]);
3969 : CPLAssert(z2[1] == zSMETransformed[1]);
3970 : CPLAssert(x2[2] == xSMETransformed[2]);
3971 : CPLAssert(y2[2] == ySMETransformed[2]);
3972 : CPLAssert(z2[2] == zSMETransformed[2]);
3973 : }
3974 : #endif
3975 :
3976 : #ifdef DEBUG_APPROX_TRANSFORMER
3977 : fprintf(stderr, "start (%.3f,%.3f) -> (%.3f,%.3f)\n", /*ok*/
3978 : x[0], y[0], xSMETransformed[0], ySMETransformed[0]);
3979 : fprintf(stderr, "middle (%.3f,%.3f) -> (%.3f,%.3f)\n", /*ok*/
3980 : x[nMiddle], y[nMiddle], xSMETransformed[1], ySMETransformed[1]);
3981 : fprintf(stderr, "end (%.3f,%.3f) -> (%.3f,%.3f)\n", /*ok*/
3982 : x[nPoints - 1], y[nPoints - 1], xSMETransformed[2],
3983 : ySMETransformed[2]);
3984 : #endif
3985 :
3986 : /* -------------------------------------------------------------------- */
3987 : /* Is the error at the middle acceptable relative to an */
3988 : /* interpolation of the middle position? */
3989 : /* -------------------------------------------------------------------- */
3990 1111520 : const double dfDeltaX =
3991 1111520 : (xSMETransformed[2] - xSMETransformed[0]) / (x[nPoints - 1] - x[0]);
3992 1111520 : const double dfDeltaY =
3993 1111520 : (ySMETransformed[2] - ySMETransformed[0]) / (x[nPoints - 1] - x[0]);
3994 1111520 : const double dfDeltaZ =
3995 1111520 : (zSMETransformed[2] - zSMETransformed[0]) / (x[nPoints - 1] - x[0]);
3996 :
3997 1111520 : const double dfError =
3998 1111520 : fabs((xSMETransformed[0] + dfDeltaX * (x[nMiddle] - x[0])) -
3999 1111520 : xSMETransformed[1]) +
4000 1111520 : fabs((ySMETransformed[0] + dfDeltaY * (x[nMiddle] - x[0])) -
4001 1111520 : ySMETransformed[1]);
4002 :
4003 1111520 : const double dfMaxError =
4004 1111520 : (bDstToSrc) ? psATInfo->dfMaxErrorReverse : psATInfo->dfMaxErrorForward;
4005 1111520 : if (dfError > dfMaxError)
4006 : {
4007 : #if DEBUG_VERBOSE
4008 : CPLDebug("GDAL",
4009 : "ApproxTransformer - "
4010 : "error %g over threshold %g, subdivide %d points.",
4011 : dfError, dfMaxError, nPoints);
4012 : #endif
4013 :
4014 589967 : double xMiddle[3] = {x[(nMiddle - 1) / 2], x[nMiddle - 1],
4015 589967 : x[nMiddle + (nPoints - nMiddle - 1) / 2]};
4016 589967 : double yMiddle[3] = {y[(nMiddle - 1) / 2], y[nMiddle - 1],
4017 589967 : y[nMiddle + (nPoints - nMiddle - 1) / 2]};
4018 589967 : double zMiddle[3] = {z[(nMiddle - 1) / 2], z[nMiddle - 1],
4019 589967 : z[nMiddle + (nPoints - nMiddle - 1) / 2]};
4020 :
4021 589967 : const bool bUseBaseTransformForHalf1 =
4022 462084 : nMiddle <= 5 || y[0] != y[nMiddle - 1] ||
4023 1514140 : y[0] != y[(nMiddle - 1) / 2] || x[0] == x[nMiddle - 1] ||
4024 462084 : x[0] == x[(nMiddle - 1) / 2];
4025 589967 : const bool bUseBaseTransformForHalf2 =
4026 476072 : nPoints - nMiddle <= 5 || y[nMiddle] != y[nPoints - 1] ||
4027 476072 : y[nMiddle] != y[nMiddle + (nPoints - nMiddle - 1) / 2] ||
4028 1542110 : x[nMiddle] == x[nPoints - 1] ||
4029 476072 : x[nMiddle] == x[nMiddle + (nPoints - nMiddle - 1) / 2];
4030 :
4031 589967 : int anSuccess2[3] = {};
4032 589967 : int bSuccess = FALSE;
4033 589967 : if (!bUseBaseTransformForHalf1 && !bUseBaseTransformForHalf2)
4034 462084 : bSuccess = psATInfo->pfnBaseTransformer(
4035 : psATInfo->pBaseCBData, bDstToSrc, 3, xMiddle, yMiddle, zMiddle,
4036 : anSuccess2);
4037 127883 : else if (!bUseBaseTransformForHalf1)
4038 : {
4039 0 : bSuccess = psATInfo->pfnBaseTransformer(
4040 : psATInfo->pBaseCBData, bDstToSrc, 2, xMiddle, yMiddle, zMiddle,
4041 : anSuccess2);
4042 0 : anSuccess2[2] = TRUE;
4043 : }
4044 127883 : else if (!bUseBaseTransformForHalf2)
4045 : {
4046 13988 : bSuccess = psATInfo->pfnBaseTransformer(
4047 : psATInfo->pBaseCBData, bDstToSrc, 1, xMiddle + 2, yMiddle + 2,
4048 : zMiddle + 2, anSuccess2 + 2);
4049 13988 : anSuccess2[0] = TRUE;
4050 13988 : anSuccess2[1] = TRUE;
4051 : }
4052 :
4053 589967 : if (!bSuccess || !anSuccess2[0] || !anSuccess2[1] || !anSuccess2[2])
4054 : {
4055 113911 : bSuccess = psATInfo->pfnBaseTransformer(
4056 : psATInfo->pBaseCBData, bDstToSrc, nMiddle - 1, x + 1, y + 1,
4057 : z + 1, panSuccess + 1);
4058 227822 : bSuccess &= psATInfo->pfnBaseTransformer(
4059 113911 : psATInfo->pBaseCBData, bDstToSrc, nPoints - nMiddle - 2,
4060 113911 : x + nMiddle + 1, y + nMiddle + 1, z + nMiddle + 1,
4061 113911 : panSuccess + nMiddle + 1);
4062 :
4063 113911 : x[0] = xSMETransformed[0];
4064 113911 : y[0] = ySMETransformed[0];
4065 113911 : z[0] = zSMETransformed[0];
4066 113911 : panSuccess[0] = TRUE;
4067 113911 : x[nMiddle] = xSMETransformed[1];
4068 113911 : y[nMiddle] = ySMETransformed[1];
4069 113911 : z[nMiddle] = zSMETransformed[1];
4070 113911 : panSuccess[nMiddle] = TRUE;
4071 113911 : x[nPoints - 1] = xSMETransformed[2];
4072 113911 : y[nPoints - 1] = ySMETransformed[2];
4073 113911 : z[nPoints - 1] = zSMETransformed[2];
4074 113911 : panSuccess[nPoints - 1] = TRUE;
4075 113911 : return bSuccess;
4076 : }
4077 :
4078 476056 : double x2[3] = {};
4079 476056 : double y2[3] = {};
4080 476056 : double z2[3] = {};
4081 476056 : if (!bUseBaseTransformForHalf1)
4082 : {
4083 462068 : x2[0] = xSMETransformed[0];
4084 462068 : y2[0] = ySMETransformed[0];
4085 462068 : z2[0] = zSMETransformed[0];
4086 462068 : x2[1] = xMiddle[0];
4087 462068 : y2[1] = yMiddle[0];
4088 462068 : z2[1] = zMiddle[0];
4089 462068 : x2[2] = xMiddle[1];
4090 462068 : y2[2] = yMiddle[1];
4091 462068 : z2[2] = zMiddle[1];
4092 :
4093 462068 : bSuccess = GDALApproxTransformInternal(
4094 : psATInfo, bDstToSrc, nMiddle, x, y, z, panSuccess, x2, y2, z2);
4095 : }
4096 : else
4097 : {
4098 13988 : bSuccess = psATInfo->pfnBaseTransformer(
4099 : psATInfo->pBaseCBData, bDstToSrc, nMiddle - 1, x + 1, y + 1,
4100 : z + 1, panSuccess + 1);
4101 13988 : x[0] = xSMETransformed[0];
4102 13988 : y[0] = ySMETransformed[0];
4103 13988 : z[0] = zSMETransformed[0];
4104 13988 : panSuccess[0] = TRUE;
4105 : }
4106 :
4107 476056 : if (!bSuccess)
4108 24 : return FALSE;
4109 :
4110 476032 : if (!bUseBaseTransformForHalf2)
4111 : {
4112 476032 : x2[0] = xSMETransformed[1];
4113 476032 : y2[0] = ySMETransformed[1];
4114 476032 : z2[0] = zSMETransformed[1];
4115 476032 : x2[1] = xMiddle[2];
4116 476032 : y2[1] = yMiddle[2];
4117 476032 : z2[1] = zMiddle[2];
4118 476032 : x2[2] = xSMETransformed[2];
4119 476032 : y2[2] = ySMETransformed[2];
4120 476032 : z2[2] = zSMETransformed[2];
4121 :
4122 476032 : bSuccess = GDALApproxTransformInternal(
4123 476032 : psATInfo, bDstToSrc, nPoints - nMiddle, x + nMiddle,
4124 476032 : y + nMiddle, z + nMiddle, panSuccess + nMiddle, x2, y2, z2);
4125 : }
4126 : else
4127 : {
4128 0 : bSuccess = psATInfo->pfnBaseTransformer(
4129 0 : psATInfo->pBaseCBData, bDstToSrc, nPoints - nMiddle - 2,
4130 0 : x + nMiddle + 1, y + nMiddle + 1, z + nMiddle + 1,
4131 0 : panSuccess + nMiddle + 1);
4132 :
4133 0 : x[nMiddle] = xSMETransformed[1];
4134 0 : y[nMiddle] = ySMETransformed[1];
4135 0 : z[nMiddle] = zSMETransformed[1];
4136 0 : panSuccess[nMiddle] = TRUE;
4137 0 : x[nPoints - 1] = xSMETransformed[2];
4138 0 : y[nPoints - 1] = ySMETransformed[2];
4139 0 : z[nPoints - 1] = zSMETransformed[2];
4140 0 : panSuccess[nPoints - 1] = TRUE;
4141 : }
4142 :
4143 476032 : if (!bSuccess)
4144 2 : return FALSE;
4145 :
4146 476030 : return TRUE;
4147 : }
4148 :
4149 : /* -------------------------------------------------------------------- */
4150 : /* Error is OK since this is just used to compute output bounds */
4151 : /* of newly created file for gdalwarper. So just use affine */
4152 : /* approximation of the reverse transform. Eventually we */
4153 : /* should implement iterative searching to find a result within */
4154 : /* our error threshold. */
4155 : /* NOTE: the above comment is not true: gdalwarp uses approximator */
4156 : /* also to compute the source pixel of each target pixel. */
4157 : /* -------------------------------------------------------------------- */
4158 97353100 : for (int i = nPoints - 1; i >= 0; i--)
4159 : {
4160 : #ifdef check_error
4161 : double xtemp = x[i];
4162 : double ytemp = y[i];
4163 : double ztemp = z[i];
4164 : double x_ori = xtemp;
4165 : double y_ori = ytemp;
4166 : int btemp = FALSE;
4167 : psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc, 1,
4168 : &xtemp, &ytemp, &ztemp, &btemp);
4169 : #endif
4170 96831500 : const double dfDist = (x[i] - x[0]);
4171 96831500 : x[i] = xSMETransformed[0] + dfDeltaX * dfDist;
4172 96831500 : y[i] = ySMETransformed[0] + dfDeltaY * dfDist;
4173 96831500 : z[i] = zSMETransformed[0] + dfDeltaZ * dfDist;
4174 : #ifdef check_error
4175 : const double dfError2 = fabs(x[i] - xtemp) + fabs(y[i] - ytemp);
4176 : if (dfError2 > 4 /*10 * dfMaxError*/)
4177 : {
4178 : /*ok*/ printf("Error = %f on (%f, %f)\n", dfError2, x_ori, y_ori);
4179 : }
4180 : #endif
4181 96831500 : panSuccess[i] = TRUE;
4182 : }
4183 :
4184 521558 : return TRUE;
4185 : }
4186 :
4187 : /************************************************************************/
4188 : /* GDALApproxTransform() */
4189 : /************************************************************************/
4190 :
4191 : /**
4192 : * Perform approximate transformation.
4193 : *
4194 : * Actually performs the approximate transformation described in
4195 : * GDALCreateApproxTransformer(). This function matches the
4196 : * GDALTransformerFunc() signature. Details of the arguments are described
4197 : * there.
4198 : */
4199 :
4200 515519 : int GDALApproxTransform(void *pCBData, int bDstToSrc, int nPoints, double *x,
4201 : double *y, double *z, int *panSuccess)
4202 :
4203 : {
4204 515519 : ApproxTransformInfo *psATInfo = static_cast<ApproxTransformInfo *>(pCBData);
4205 515519 : double x2[3] = {};
4206 515519 : double y2[3] = {};
4207 515519 : double z2[3] = {};
4208 515519 : int anSuccess2[3] = {};
4209 : int bSuccess;
4210 :
4211 515519 : const int nMiddle = (nPoints - 1) / 2;
4212 :
4213 : /* -------------------------------------------------------------------- */
4214 : /* Bail if our preconditions are not met, or if error is not */
4215 : /* acceptable. */
4216 : /* -------------------------------------------------------------------- */
4217 515519 : int bRet = FALSE;
4218 515519 : if (y[0] != y[nPoints - 1] || y[0] != y[nMiddle] ||
4219 507298 : x[0] == x[nPoints - 1] || x[0] == x[nMiddle] ||
4220 178962 : (psATInfo->dfMaxErrorForward == 0.0 &&
4221 178962 : psATInfo->dfMaxErrorReverse == 0.0) ||
4222 : nPoints <= 5)
4223 : {
4224 336971 : bRet = psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc,
4225 : nPoints, x, y, z, panSuccess);
4226 336972 : goto end;
4227 : }
4228 :
4229 : /* -------------------------------------------------------------------- */
4230 : /* Transform first, last and middle point. */
4231 : /* -------------------------------------------------------------------- */
4232 178548 : x2[0] = x[0];
4233 178548 : y2[0] = y[0];
4234 178548 : z2[0] = z[0];
4235 178548 : x2[1] = x[nMiddle];
4236 178548 : y2[1] = y[nMiddle];
4237 178548 : z2[1] = z[nMiddle];
4238 178548 : x2[2] = x[nPoints - 1];
4239 178548 : y2[2] = y[nPoints - 1];
4240 178548 : z2[2] = z[nPoints - 1];
4241 :
4242 178548 : bSuccess = psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc, 3,
4243 : x2, y2, z2, anSuccess2);
4244 178548 : if (!bSuccess || !anSuccess2[0] || !anSuccess2[1] || !anSuccess2[2])
4245 : {
4246 5123 : bRet = psATInfo->pfnBaseTransformer(psATInfo->pBaseCBData, bDstToSrc,
4247 : nPoints, x, y, z, panSuccess);
4248 5123 : goto end;
4249 : }
4250 :
4251 173425 : bRet = GDALApproxTransformInternal(pCBData, bDstToSrc, nPoints, x, y, z,
4252 : panSuccess, x2, y2, z2);
4253 :
4254 515520 : end:
4255 : #ifdef DEBUG_APPROX_TRANSFORMER
4256 : for (int i = 0; i < nPoints; i++)
4257 : fprintf(stderr, "[%d] (%.10f,%.10f) %d\n", /*ok*/
4258 : i, x[i], y[i], panSuccess[i]);
4259 : #endif
4260 :
4261 515520 : return bRet;
4262 : }
4263 :
4264 : /************************************************************************/
4265 : /* GDALDeserializeApproxTransformer() */
4266 : /************************************************************************/
4267 :
4268 150 : static void *GDALDeserializeApproxTransformer(CPLXMLNode *psTree)
4269 :
4270 : {
4271 150 : double dfMaxErrorForward = 0.25;
4272 150 : double dfMaxErrorReverse = 0.25;
4273 150 : const char *pszMaxError = CPLGetXMLValue(psTree, "MaxError", nullptr);
4274 150 : if (pszMaxError != nullptr)
4275 : {
4276 148 : dfMaxErrorForward = CPLAtof(pszMaxError);
4277 148 : dfMaxErrorReverse = dfMaxErrorForward;
4278 : }
4279 : const char *pszMaxErrorForward =
4280 150 : CPLGetXMLValue(psTree, "MaxErrorForward", nullptr);
4281 150 : if (pszMaxErrorForward != nullptr)
4282 : {
4283 2 : dfMaxErrorForward = CPLAtof(pszMaxErrorForward);
4284 : }
4285 : const char *pszMaxErrorReverse =
4286 150 : CPLGetXMLValue(psTree, "MaxErrorReverse", nullptr);
4287 150 : if (pszMaxErrorReverse != nullptr)
4288 : {
4289 2 : dfMaxErrorReverse = CPLAtof(pszMaxErrorReverse);
4290 : }
4291 :
4292 150 : GDALTransformerFunc pfnBaseTransform = nullptr;
4293 150 : void *pBaseCBData = nullptr;
4294 :
4295 150 : CPLXMLNode *psContainer = CPLGetXMLNode(psTree, "BaseTransformer");
4296 :
4297 150 : if (psContainer != nullptr && psContainer->psChild != nullptr)
4298 : {
4299 150 : GDALDeserializeTransformer(psContainer->psChild, &pfnBaseTransform,
4300 : &pBaseCBData);
4301 : }
4302 :
4303 150 : if (pfnBaseTransform == nullptr)
4304 : {
4305 0 : CPLError(CE_Failure, CPLE_AppDefined,
4306 : "Cannot get base transform for approx transformer.");
4307 0 : return nullptr;
4308 : }
4309 :
4310 150 : void *pApproxCBData = GDALCreateApproxTransformer2(
4311 : pfnBaseTransform, pBaseCBData, dfMaxErrorForward, dfMaxErrorReverse);
4312 150 : GDALApproxTransformerOwnsSubtransformer(pApproxCBData, TRUE);
4313 :
4314 150 : return pApproxCBData;
4315 : }
4316 :
4317 : /************************************************************************/
4318 : /* GDALTransformLonLatToDestApproxTransformer() */
4319 : /************************************************************************/
4320 :
4321 2130 : int GDALTransformLonLatToDestApproxTransformer(void *hTransformArg,
4322 : double *pdfX, double *pdfY)
4323 : {
4324 2130 : ApproxTransformInfo *psInfo =
4325 : static_cast<ApproxTransformInfo *>(hTransformArg);
4326 :
4327 2130 : if (GDALIsTransformer(psInfo->pBaseCBData,
4328 : GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
4329 : {
4330 2130 : return GDALTransformLonLatToDestGenImgProjTransformer(
4331 2130 : psInfo->pBaseCBData, pdfX, pdfY);
4332 : }
4333 0 : return false;
4334 : }
4335 :
4336 : /************************************************************************/
4337 : /* GDALApplyGeoTransform() */
4338 : /************************************************************************/
4339 :
4340 : /**
4341 : * Apply GeoTransform to x/y coordinate.
4342 : *
4343 : * Applies the following computation, converting a (pixel, line) coordinate
4344 : * into a georeferenced (geo_x, geo_y) location.
4345 : * \code{.c}
4346 : * *pdfGeoX = padfGeoTransform[0] + dfPixel * padfGeoTransform[1]
4347 : * + dfLine * padfGeoTransform[2];
4348 : * *pdfGeoY = padfGeoTransform[3] + dfPixel * padfGeoTransform[4]
4349 : * + dfLine * padfGeoTransform[5];
4350 : * \endcode
4351 : *
4352 : * @param padfGeoTransform Six coefficient GeoTransform to apply.
4353 : * @param dfPixel Input pixel position.
4354 : * @param dfLine Input line position.
4355 : * @param pdfGeoX output location where geo_x (easting/longitude)
4356 : * location is placed.
4357 : * @param pdfGeoY output location where geo_y (northing/latitude)
4358 : * location is placed.
4359 : */
4360 :
4361 304748 : void CPL_STDCALL GDALApplyGeoTransform(const double *padfGeoTransform,
4362 : double dfPixel, double dfLine,
4363 : double *pdfGeoX, double *pdfGeoY)
4364 : {
4365 304748 : *pdfGeoX = padfGeoTransform[0] + dfPixel * padfGeoTransform[1] +
4366 304748 : dfLine * padfGeoTransform[2];
4367 304748 : *pdfGeoY = padfGeoTransform[3] + dfPixel * padfGeoTransform[4] +
4368 304748 : dfLine * padfGeoTransform[5];
4369 304748 : }
4370 :
4371 : /************************************************************************/
4372 : /* GDALInvGeoTransform() */
4373 : /************************************************************************/
4374 :
4375 : /**
4376 : * Invert Geotransform.
4377 : *
4378 : * This function will invert a standard 3x2 set of GeoTransform coefficients.
4379 : * This converts the equation from being pixel to geo to being geo to pixel.
4380 : *
4381 : * @param gt_in Input geotransform (six doubles - unaltered).
4382 : * @param gt_out Output geotransform (six doubles - updated).
4383 : *
4384 : * @return TRUE on success or FALSE if the equation is uninvertable.
4385 : */
4386 :
4387 2986 : int CPL_STDCALL GDALInvGeoTransform(const double *gt_in, double *gt_out)
4388 :
4389 : {
4390 : // Special case - no rotation - to avoid computing determinate
4391 : // and potential precision issues.
4392 2986 : if (gt_in[2] == 0.0 && gt_in[4] == 0.0 && gt_in[1] != 0.0 &&
4393 2932 : gt_in[5] != 0.0)
4394 : {
4395 : /*X = gt_in[0] + x * gt_in[1]
4396 : Y = gt_in[3] + y * gt_in[5]
4397 : -->
4398 : x = -gt_in[0] / gt_in[1] + (1 / gt_in[1]) * X
4399 : y = -gt_in[3] / gt_in[5] + (1 / gt_in[5]) * Y
4400 : */
4401 2932 : gt_out[0] = -gt_in[0] / gt_in[1];
4402 2932 : gt_out[1] = 1.0 / gt_in[1];
4403 2932 : gt_out[2] = 0.0;
4404 2932 : gt_out[3] = -gt_in[3] / gt_in[5];
4405 2932 : gt_out[4] = 0.0;
4406 2932 : gt_out[5] = 1.0 / gt_in[5];
4407 2932 : return 1;
4408 : }
4409 :
4410 : // Assume a 3rd row that is [1 0 0].
4411 :
4412 : // Compute determinate.
4413 :
4414 54 : const double det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4];
4415 108 : const double magnitude = std::max(std::max(fabs(gt_in[1]), fabs(gt_in[2])),
4416 54 : std::max(fabs(gt_in[4]), fabs(gt_in[5])));
4417 :
4418 54 : if (fabs(det) <= 1e-10 * magnitude * magnitude)
4419 5 : return 0;
4420 :
4421 49 : const double inv_det = 1.0 / det;
4422 :
4423 : // Compute adjoint, and divide by determinate.
4424 :
4425 49 : gt_out[1] = gt_in[5] * inv_det;
4426 49 : gt_out[4] = -gt_in[4] * inv_det;
4427 :
4428 49 : gt_out[2] = -gt_in[2] * inv_det;
4429 49 : gt_out[5] = gt_in[1] * inv_det;
4430 :
4431 49 : gt_out[0] = (gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det;
4432 49 : gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det;
4433 :
4434 49 : return 1;
4435 : }
4436 :
4437 : /************************************************************************/
4438 : /* GDALSerializeTransformer() */
4439 : /************************************************************************/
4440 :
4441 263 : CPLXMLNode *GDALSerializeTransformer(GDALTransformerFunc /* pfnFunc */,
4442 : void *pTransformArg)
4443 : {
4444 263 : VALIDATE_POINTER1(pTransformArg, "GDALSerializeTransformer", nullptr);
4445 :
4446 263 : GDALTransformerInfo *psInfo =
4447 : static_cast<GDALTransformerInfo *>(pTransformArg);
4448 :
4449 263 : if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4450 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4451 : {
4452 0 : CPLError(CE_Failure, CPLE_AppDefined,
4453 : "Attempt to serialize non-GTI2 transformer.");
4454 0 : return nullptr;
4455 : }
4456 263 : else if (psInfo->pfnSerialize == nullptr)
4457 : {
4458 0 : CPLError(CE_Failure, CPLE_AppDefined,
4459 : "No serialization function available for this transformer.");
4460 0 : return nullptr;
4461 : }
4462 :
4463 263 : return psInfo->pfnSerialize(pTransformArg);
4464 : }
4465 :
4466 : /************************************************************************/
4467 : /* GDALRegisterTransformDeserializer() */
4468 : /************************************************************************/
4469 :
4470 : static CPLList *psListDeserializer = nullptr;
4471 : static CPLMutex *hDeserializerMutex = nullptr;
4472 :
4473 : typedef struct
4474 : {
4475 : char *pszTransformName;
4476 : GDALTransformerFunc pfnTransformerFunc;
4477 : GDALTransformDeserializeFunc pfnDeserializeFunc;
4478 : } TransformDeserializerInfo;
4479 :
4480 0 : void *GDALRegisterTransformDeserializer(
4481 : const char *pszTransformName, GDALTransformerFunc pfnTransformerFunc,
4482 : GDALTransformDeserializeFunc pfnDeserializeFunc)
4483 : {
4484 : TransformDeserializerInfo *psInfo =
4485 : static_cast<TransformDeserializerInfo *>(
4486 0 : CPLMalloc(sizeof(TransformDeserializerInfo)));
4487 0 : psInfo->pszTransformName = CPLStrdup(pszTransformName);
4488 0 : psInfo->pfnTransformerFunc = pfnTransformerFunc;
4489 0 : psInfo->pfnDeserializeFunc = pfnDeserializeFunc;
4490 :
4491 0 : CPLMutexHolderD(&hDeserializerMutex);
4492 0 : psListDeserializer = CPLListInsert(psListDeserializer, psInfo, 0);
4493 :
4494 0 : return psInfo;
4495 : }
4496 :
4497 : /************************************************************************/
4498 : /* GDALUnregisterTransformDeserializer() */
4499 : /************************************************************************/
4500 :
4501 0 : void GDALUnregisterTransformDeserializer(void *pData)
4502 : {
4503 0 : CPLMutexHolderD(&hDeserializerMutex);
4504 0 : CPLList *psList = psListDeserializer;
4505 0 : CPLList *psLast = nullptr;
4506 0 : while (psList)
4507 : {
4508 0 : if (psList->pData == pData)
4509 : {
4510 0 : TransformDeserializerInfo *psInfo =
4511 : static_cast<TransformDeserializerInfo *>(pData);
4512 0 : CPLFree(psInfo->pszTransformName);
4513 0 : CPLFree(pData);
4514 0 : if (psLast)
4515 0 : psLast->psNext = psList->psNext;
4516 : else
4517 0 : psListDeserializer = nullptr;
4518 0 : CPLFree(psList);
4519 0 : break;
4520 : }
4521 0 : psLast = psList;
4522 0 : psList = psList->psNext;
4523 : }
4524 0 : }
4525 :
4526 : /************************************************************************/
4527 : /* GDALUnregisterTransformDeserializer() */
4528 : /************************************************************************/
4529 :
4530 943 : void GDALCleanupTransformDeserializerMutex()
4531 : {
4532 943 : if (hDeserializerMutex != nullptr)
4533 : {
4534 0 : CPLDestroyMutex(hDeserializerMutex);
4535 0 : hDeserializerMutex = nullptr;
4536 : }
4537 943 : }
4538 :
4539 : /************************************************************************/
4540 : /* GDALDeserializeTransformer() */
4541 : /************************************************************************/
4542 :
4543 538 : CPLErr GDALDeserializeTransformer(CPLXMLNode *psTree,
4544 : GDALTransformerFunc *ppfnFunc,
4545 : void **ppTransformArg)
4546 :
4547 : {
4548 538 : *ppfnFunc = nullptr;
4549 538 : *ppTransformArg = nullptr;
4550 :
4551 538 : CPLErrorReset();
4552 :
4553 538 : if (psTree == nullptr || psTree->eType != CXT_Element)
4554 0 : CPLError(CE_Failure, CPLE_AppDefined,
4555 : "Malformed element in GDALDeserializeTransformer");
4556 538 : else if (EQUAL(psTree->pszValue, "GenImgProjTransformer"))
4557 : {
4558 198 : *ppfnFunc = GDALGenImgProjTransform;
4559 198 : *ppTransformArg = GDALDeserializeGenImgProjTransformer(psTree);
4560 : }
4561 340 : else if (EQUAL(psTree->pszValue, "ReprojectionTransformer"))
4562 : {
4563 174 : *ppfnFunc = GDALReprojectionTransform;
4564 174 : *ppTransformArg = GDALDeserializeReprojectionTransformer(psTree);
4565 : }
4566 166 : else if (EQUAL(psTree->pszValue, "GCPTransformer"))
4567 : {
4568 12 : *ppfnFunc = GDALGCPTransform;
4569 12 : *ppTransformArg = GDALDeserializeGCPTransformer(psTree);
4570 : }
4571 154 : else if (EQUAL(psTree->pszValue, "TPSTransformer"))
4572 : {
4573 3 : *ppfnFunc = GDALTPSTransform;
4574 3 : *ppTransformArg = GDALDeserializeTPSTransformer(psTree);
4575 : }
4576 151 : else if (EQUAL(psTree->pszValue, "GeoLocTransformer"))
4577 : {
4578 1 : *ppfnFunc = GDALGeoLocTransform;
4579 1 : *ppTransformArg = GDALDeserializeGeoLocTransformer(psTree);
4580 : }
4581 150 : else if (EQUAL(psTree->pszValue, "RPCTransformer"))
4582 : {
4583 0 : *ppfnFunc = GDALRPCTransform;
4584 0 : *ppTransformArg = GDALDeserializeRPCTransformer(psTree);
4585 : }
4586 150 : else if (EQUAL(psTree->pszValue, "ApproxTransformer"))
4587 : {
4588 150 : *ppfnFunc = GDALApproxTransform;
4589 150 : *ppTransformArg = GDALDeserializeApproxTransformer(psTree);
4590 : }
4591 : else
4592 : {
4593 0 : GDALTransformDeserializeFunc pfnDeserializeFunc = nullptr;
4594 : {
4595 0 : CPLMutexHolderD(&hDeserializerMutex);
4596 0 : CPLList *psList = psListDeserializer;
4597 0 : while (psList)
4598 : {
4599 0 : TransformDeserializerInfo *psInfo =
4600 : static_cast<TransformDeserializerInfo *>(psList->pData);
4601 0 : if (strcmp(psInfo->pszTransformName, psTree->pszValue) == 0)
4602 : {
4603 0 : *ppfnFunc = psInfo->pfnTransformerFunc;
4604 0 : pfnDeserializeFunc = psInfo->pfnDeserializeFunc;
4605 0 : break;
4606 : }
4607 0 : psList = psList->psNext;
4608 : }
4609 : }
4610 :
4611 0 : if (pfnDeserializeFunc != nullptr)
4612 : {
4613 0 : *ppTransformArg = pfnDeserializeFunc(psTree);
4614 : }
4615 : else
4616 : {
4617 0 : CPLError(CE_Failure, CPLE_AppDefined,
4618 : "Unrecognized element '%s' GDALDeserializeTransformer",
4619 : psTree->pszValue);
4620 : }
4621 : }
4622 :
4623 538 : return CPLGetLastErrorType();
4624 : }
4625 :
4626 : /************************************************************************/
4627 : /* GDALDestroyTransformer() */
4628 : /************************************************************************/
4629 :
4630 3535 : void GDALDestroyTransformer(void *pTransformArg)
4631 :
4632 : {
4633 3535 : if (pTransformArg == nullptr)
4634 10 : return;
4635 :
4636 3525 : GDALTransformerInfo *psInfo =
4637 : static_cast<GDALTransformerInfo *>(pTransformArg);
4638 :
4639 3525 : if (memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4640 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4641 : {
4642 0 : CPLError(CE_Failure, CPLE_AppDefined,
4643 : "Attempt to destroy non-GTI2 transformer.");
4644 0 : return;
4645 : }
4646 :
4647 3525 : psInfo->pfnCleanup(pTransformArg);
4648 : }
4649 :
4650 : /************************************************************************/
4651 : /* GDALUseTransformer() */
4652 : /************************************************************************/
4653 :
4654 8680 : int GDALUseTransformer(void *pTransformArg, int bDstToSrc, int nPointCount,
4655 : double *x, double *y, double *z, int *panSuccess)
4656 : {
4657 8680 : GDALTransformerInfo *psInfo =
4658 : static_cast<GDALTransformerInfo *>(pTransformArg);
4659 :
4660 8680 : if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4661 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4662 : {
4663 0 : CPLError(CE_Failure, CPLE_AppDefined,
4664 : "Attempt to use non-GTI2 transformer.");
4665 0 : return FALSE;
4666 : }
4667 :
4668 8680 : return psInfo->pfnTransform(pTransformArg, bDstToSrc, nPointCount, x, y, z,
4669 8680 : panSuccess);
4670 : }
4671 :
4672 : /************************************************************************/
4673 : /* GDALCloneTransformer() */
4674 : /************************************************************************/
4675 :
4676 23 : void *GDALCloneTransformer(void *pTransformArg)
4677 : {
4678 23 : GDALTransformerInfo *psInfo =
4679 : static_cast<GDALTransformerInfo *>(pTransformArg);
4680 :
4681 23 : if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4682 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4683 : {
4684 0 : CPLError(CE_Failure, CPLE_AppDefined,
4685 : "Attempt to clone non-GTI2 transformer.");
4686 0 : return nullptr;
4687 : }
4688 :
4689 23 : if (psInfo->pfnCreateSimilar != nullptr)
4690 : {
4691 12 : return psInfo->pfnCreateSimilar(psInfo, 1.0, 1.0);
4692 : }
4693 :
4694 11 : if (psInfo->pfnSerialize == nullptr)
4695 : {
4696 0 : CPLError(CE_Failure, CPLE_AppDefined,
4697 : "No serialization function available for this transformer.");
4698 0 : return nullptr;
4699 : }
4700 :
4701 11 : CPLXMLNode *pSerialized = psInfo->pfnSerialize(pTransformArg);
4702 11 : if (pSerialized == nullptr)
4703 0 : return nullptr;
4704 11 : GDALTransformerFunc pfnTransformer = nullptr;
4705 11 : void *pClonedTransformArg = nullptr;
4706 11 : if (GDALDeserializeTransformer(pSerialized, &pfnTransformer,
4707 11 : &pClonedTransformArg) != CE_None)
4708 : {
4709 0 : CPLDestroyXMLNode(pSerialized);
4710 0 : CPLFree(pClonedTransformArg);
4711 0 : return nullptr;
4712 : }
4713 :
4714 11 : CPLDestroyXMLNode(pSerialized);
4715 11 : return pClonedTransformArg;
4716 : }
4717 :
4718 : /************************************************************************/
4719 : /* GDALCreateSimilarTransformer() */
4720 : /************************************************************************/
4721 :
4722 44 : void *GDALCreateSimilarTransformer(void *pTransformArg, double dfRatioX,
4723 : double dfRatioY)
4724 : {
4725 44 : GDALTransformerInfo *psInfo =
4726 : static_cast<GDALTransformerInfo *>(pTransformArg);
4727 :
4728 44 : if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4729 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4730 : {
4731 0 : CPLError(CE_Failure, CPLE_AppDefined,
4732 : "Attempt to call CreateSimilar on a non-GTI2 transformer.");
4733 0 : return nullptr;
4734 : }
4735 :
4736 44 : if (psInfo->pfnCreateSimilar == nullptr)
4737 : {
4738 0 : CPLError(CE_Failure, CPLE_AppDefined,
4739 : "No CreateSimilar function available for this transformer.");
4740 0 : return nullptr;
4741 : }
4742 :
4743 44 : return psInfo->pfnCreateSimilar(psInfo, dfRatioX, dfRatioY);
4744 : }
4745 :
4746 : /************************************************************************/
4747 : /* GetGenImgProjTransformInfo() */
4748 : /************************************************************************/
4749 :
4750 44 : static GDALTransformerInfo *GetGenImgProjTransformInfo(const char *pszFunc,
4751 : void *pTransformArg)
4752 : {
4753 44 : GDALTransformerInfo *psInfo =
4754 : static_cast<GDALTransformerInfo *>(pTransformArg);
4755 :
4756 44 : if (psInfo == nullptr || memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4757 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4758 : {
4759 0 : CPLError(CE_Failure, CPLE_AppDefined,
4760 : "Attempt to call %s on "
4761 : "a non-GTI2 transformer.",
4762 : pszFunc);
4763 0 : return nullptr;
4764 : }
4765 :
4766 44 : if (EQUAL(psInfo->pszClassName, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
4767 : {
4768 14 : ApproxTransformInfo *psATInfo =
4769 : static_cast<ApproxTransformInfo *>(pTransformArg);
4770 14 : psInfo = static_cast<GDALTransformerInfo *>(psATInfo->pBaseCBData);
4771 :
4772 14 : if (psInfo == nullptr ||
4773 14 : memcmp(psInfo->abySignature, GDAL_GTI2_SIGNATURE,
4774 : strlen(GDAL_GTI2_SIGNATURE)) != 0)
4775 : {
4776 0 : CPLError(CE_Failure, CPLE_AppDefined,
4777 : "Attempt to call %s on "
4778 : "a non-GTI2 transformer.",
4779 : pszFunc);
4780 0 : return nullptr;
4781 : }
4782 : }
4783 :
4784 44 : if (EQUAL(psInfo->pszClassName, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
4785 : {
4786 44 : return psInfo;
4787 : }
4788 :
4789 0 : return nullptr;
4790 : }
4791 :
4792 : /************************************************************************/
4793 : /* GDALSetTransformerDstGeoTransform() */
4794 : /************************************************************************/
4795 :
4796 : /**
4797 : * Set ApproxTransformer or GenImgProj output geotransform.
4798 : *
4799 : * This is a layer above GDALSetGenImgProjTransformerDstGeoTransform() that
4800 : * checks that the passed hTransformArg is compatible.
4801 : *
4802 : * Normally the "destination geotransform", or transformation between
4803 : * georeferenced output coordinates and pixel/line coordinates on the
4804 : * destination file is extracted from the destination file by
4805 : * GDALCreateGenImgProjTransformer() and stored in the GenImgProj private
4806 : * info. However, sometimes it is inconvenient to have an output file
4807 : * handle with appropriate geotransform information when creating the
4808 : * transformation. For these cases, this function can be used to apply
4809 : * the destination geotransform.
4810 : *
4811 : * @param pTransformArg the handle to update.
4812 : * @param padfGeoTransform the destination geotransform to apply (six doubles).
4813 : */
4814 :
4815 22 : void GDALSetTransformerDstGeoTransform(void *pTransformArg,
4816 : const double *padfGeoTransform)
4817 : {
4818 22 : VALIDATE_POINTER0(pTransformArg, "GDALSetTransformerDstGeoTransform");
4819 :
4820 22 : GDALTransformerInfo *psInfo = GetGenImgProjTransformInfo(
4821 : "GDALSetTransformerDstGeoTransform", pTransformArg);
4822 22 : if (psInfo)
4823 : {
4824 22 : GDALSetGenImgProjTransformerDstGeoTransform(psInfo, padfGeoTransform);
4825 : }
4826 : }
4827 :
4828 : /************************************************************************/
4829 : /* GDALGetTransformerDstGeoTransform() */
4830 : /************************************************************************/
4831 :
4832 : /**
4833 : * Get ApproxTransformer or GenImgProj output geotransform.
4834 : *
4835 : * @param pTransformArg transformer handle.
4836 : * @param padfGeoTransform (output) the destination geotransform to return (six
4837 : * doubles).
4838 : */
4839 :
4840 22 : void GDALGetTransformerDstGeoTransform(void *pTransformArg,
4841 : double *padfGeoTransform)
4842 : {
4843 22 : VALIDATE_POINTER0(pTransformArg, "GDALGetTransformerDstGeoTransform");
4844 :
4845 22 : GDALTransformerInfo *psInfo = GetGenImgProjTransformInfo(
4846 : "GDALGetTransformerDstGeoTransform", pTransformArg);
4847 22 : if (psInfo)
4848 : {
4849 22 : GDALGenImgProjTransformInfo *psGenImgProjInfo =
4850 : reinterpret_cast<GDALGenImgProjTransformInfo *>(psInfo);
4851 :
4852 22 : memcpy(padfGeoTransform, psGenImgProjInfo->adfDstGeoTransform,
4853 : sizeof(double) * 6);
4854 : }
4855 : }
4856 :
4857 : /************************************************************************/
4858 : /* GDALTransformIsTranslationOnPixelBoundaries() */
4859 : /************************************************************************/
4860 :
4861 1421 : bool GDALTransformIsTranslationOnPixelBoundaries(GDALTransformerFunc,
4862 : void *pTransformerArg)
4863 : {
4864 1421 : if (GDALIsTransformer(pTransformerArg, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
4865 : {
4866 1065 : const auto *pApproxInfo =
4867 : static_cast<const ApproxTransformInfo *>(pTransformerArg);
4868 1065 : pTransformerArg = pApproxInfo->pBaseCBData;
4869 : }
4870 1421 : if (GDALIsTransformer(pTransformerArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
4871 : {
4872 1268 : const auto *pGenImgpProjInfo =
4873 : static_cast<GDALGenImgProjTransformInfo *>(pTransformerArg);
4874 394 : const auto IsCloseToInteger = [](double dfVal)
4875 394 : { return std::fabs(dfVal - std::round(dfVal)) <= 1e-6; };
4876 2457 : return pGenImgpProjInfo->pSrcTransformArg == nullptr &&
4877 1189 : pGenImgpProjInfo->pDstTransformArg == nullptr &&
4878 1187 : pGenImgpProjInfo->pReproject == nullptr &&
4879 492 : pGenImgpProjInfo->adfSrcGeoTransform[1] ==
4880 492 : pGenImgpProjInfo->adfDstGeoTransform[1] &&
4881 249 : pGenImgpProjInfo->adfSrcGeoTransform[5] ==
4882 249 : pGenImgpProjInfo->adfDstGeoTransform[5] &&
4883 215 : pGenImgpProjInfo->adfSrcGeoTransform[2] ==
4884 215 : pGenImgpProjInfo->adfDstGeoTransform[2] &&
4885 215 : pGenImgpProjInfo->adfSrcGeoTransform[4] ==
4886 430 : pGenImgpProjInfo->adfDstGeoTransform[4] &&
4887 : // Check that the georeferenced origin of the destination
4888 : // geotransform is close to be an integer value when transformed
4889 : // to source image coordinates
4890 215 : IsCloseToInteger(
4891 215 : pGenImgpProjInfo->adfSrcInvGeoTransform[0] +
4892 215 : pGenImgpProjInfo->adfDstGeoTransform[0] *
4893 215 : pGenImgpProjInfo->adfSrcInvGeoTransform[1] +
4894 215 : pGenImgpProjInfo->adfDstGeoTransform[3] *
4895 2672 : pGenImgpProjInfo->adfSrcInvGeoTransform[2]) &&
4896 179 : IsCloseToInteger(pGenImgpProjInfo->adfSrcInvGeoTransform[3] +
4897 179 : pGenImgpProjInfo->adfDstGeoTransform[0] *
4898 179 : pGenImgpProjInfo->adfSrcInvGeoTransform[4] +
4899 179 : pGenImgpProjInfo->adfDstGeoTransform[3] *
4900 1447 : pGenImgpProjInfo->adfSrcInvGeoTransform[5]);
4901 : }
4902 153 : return false;
4903 : }
4904 :
4905 : /************************************************************************/
4906 : /* GDALTransformIsAffineNoRotation() */
4907 : /************************************************************************/
4908 :
4909 18 : bool GDALTransformIsAffineNoRotation(GDALTransformerFunc, void *pTransformerArg)
4910 : {
4911 18 : if (GDALIsTransformer(pTransformerArg, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
4912 : {
4913 18 : const auto *pApproxInfo =
4914 : static_cast<const ApproxTransformInfo *>(pTransformerArg);
4915 18 : pTransformerArg = pApproxInfo->pBaseCBData;
4916 : }
4917 18 : if (GDALIsTransformer(pTransformerArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
4918 : {
4919 18 : const auto *pGenImgpProjInfo =
4920 : static_cast<GDALGenImgProjTransformInfo *>(pTransformerArg);
4921 36 : return pGenImgpProjInfo->pSrcTransformArg == nullptr &&
4922 18 : pGenImgpProjInfo->pDstTransformArg == nullptr &&
4923 18 : pGenImgpProjInfo->pReproject == nullptr &&
4924 8 : pGenImgpProjInfo->adfSrcGeoTransform[2] == 0 &&
4925 8 : pGenImgpProjInfo->adfSrcGeoTransform[4] == 0 &&
4926 44 : pGenImgpProjInfo->adfDstGeoTransform[2] == 0 &&
4927 26 : pGenImgpProjInfo->adfDstGeoTransform[4] == 0;
4928 : }
4929 0 : return false;
4930 : }
4931 :
4932 : /************************************************************************/
4933 : /* GDALTransformHasFastClone() */
4934 : /************************************************************************/
4935 :
4936 : /** Returns whether GDALCloneTransformer() on this transformer is
4937 : * "fast"
4938 : * Counter-examples are GCPs or TPSs transformers.
4939 : */
4940 2 : bool GDALTransformHasFastClone(void *pTransformerArg)
4941 : {
4942 2 : if (GDALIsTransformer(pTransformerArg, GDAL_APPROX_TRANSFORMER_CLASS_NAME))
4943 : {
4944 1 : const auto *pApproxInfo =
4945 : static_cast<const ApproxTransformInfo *>(pTransformerArg);
4946 1 : pTransformerArg = pApproxInfo->pBaseCBData;
4947 : // Fallback to next lines
4948 : }
4949 :
4950 2 : if (GDALIsTransformer(pTransformerArg, GDAL_GEN_IMG_TRANSFORMER_CLASS_NAME))
4951 : {
4952 2 : const auto *pGenImgpProjInfo =
4953 : static_cast<GDALGenImgProjTransformInfo *>(pTransformerArg);
4954 2 : return (pGenImgpProjInfo->pSrcTransformArg == nullptr ||
4955 0 : GDALTransformHasFastClone(
4956 4 : pGenImgpProjInfo->pSrcTransformArg)) &&
4957 2 : (pGenImgpProjInfo->pDstTransformArg == nullptr ||
4958 2 : GDALTransformHasFastClone(pGenImgpProjInfo->pDstTransformArg));
4959 : }
4960 0 : else if (GDALIsTransformer(pTransformerArg,
4961 : GDAL_RPC_TRANSFORMER_CLASS_NAME))
4962 : {
4963 0 : return true;
4964 : }
4965 : else
4966 : {
4967 0 : return false;
4968 : }
4969 : }
|
