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