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