1 : /******************************************************************************
2 : * $Id: gdalgeoloc.cpp 25403 2012-12-30 18:14:28Z rouault $
3 : *
4 : * Project: GDAL
5 : * Purpose: Implements Geolocation array based transformer.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2006, Frank Warmerdam <warmerdam@pobox.com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "gdal_priv.h"
31 : #include "gdal_alg.h"
32 :
33 : #ifdef SHAPE_DEBUG
34 : #include "/u/pkg/shapelib/shapefil.h"
35 :
36 : SHPHandle hSHP = NULL;
37 : DBFHandle hDBF = NULL;
38 : #endif
39 :
40 : CPL_CVSID("$Id: gdalgeoloc.cpp 25403 2012-12-30 18:14:28Z rouault $");
41 :
42 : CPL_C_START
43 : CPLXMLNode *GDALSerializeGeoLocTransformer( void *pTransformArg );
44 : void *GDALDeserializeGeoLocTransformer( CPLXMLNode *psTree );
45 : CPL_C_END
46 :
47 : /************************************************************************/
48 : /* ==================================================================== */
49 : /* GDALGeoLocTransformer */
50 : /* ==================================================================== */
51 : /************************************************************************/
52 :
53 : typedef struct {
54 :
55 : GDALTransformerInfo sTI;
56 :
57 : int bReversed;
58 :
59 : // Map from target georef coordinates back to geolocation array
60 : // pixel line coordinates. Built only if needed.
61 :
62 : int nBackMapWidth;
63 : int nBackMapHeight;
64 : double adfBackMapGeoTransform[6]; // maps georef to pixel/line.
65 : float *pafBackMapX;
66 : float *pafBackMapY;
67 :
68 : // geolocation bands.
69 :
70 : GDALDatasetH hDS_X;
71 : GDALRasterBandH hBand_X;
72 : GDALDatasetH hDS_Y;
73 : GDALRasterBandH hBand_Y;
74 :
75 : // Located geolocation data.
76 : int nGeoLocXSize;
77 : int nGeoLocYSize;
78 : double *padfGeoLocX;
79 : double *padfGeoLocY;
80 :
81 : double dfNoDataX;
82 : double dfNoDataY;
83 :
84 : // geolocation <-> base image mapping.
85 : double dfPIXEL_OFFSET;
86 : double dfPIXEL_STEP;
87 : double dfLINE_OFFSET;
88 : double dfLINE_STEP;
89 :
90 : char ** papszGeolocationInfo;
91 :
92 : } GDALGeoLocTransformInfo;
93 :
94 : /************************************************************************/
95 : /* GeoLocLoadFullData() */
96 : /************************************************************************/
97 :
98 2 : static int GeoLocLoadFullData( GDALGeoLocTransformInfo *psTransform )
99 :
100 : {
101 : int nXSize, nYSize;
102 :
103 2 : int nXSize_XBand = GDALGetRasterXSize( psTransform->hDS_X );
104 2 : int nYSize_XBand = GDALGetRasterYSize( psTransform->hDS_X );
105 2 : int nXSize_YBand = GDALGetRasterXSize( psTransform->hDS_Y );
106 2 : int nYSize_YBand = GDALGetRasterYSize( psTransform->hDS_Y );
107 2 : if (nYSize_XBand == 1 && nYSize_YBand == 1)
108 : {
109 0 : nXSize = nXSize_XBand;
110 0 : nYSize = nXSize_YBand;
111 : }
112 : else
113 : {
114 2 : nXSize = nXSize_XBand;
115 2 : nYSize = nYSize_XBand;
116 : }
117 :
118 2 : psTransform->nGeoLocXSize = nXSize;
119 2 : psTransform->nGeoLocYSize = nYSize;
120 :
121 : psTransform->padfGeoLocY = (double *)
122 2 : VSIMalloc3(sizeof(double), nXSize, nYSize);
123 : psTransform->padfGeoLocX = (double *)
124 2 : VSIMalloc3(sizeof(double), nXSize, nYSize);
125 :
126 2 : if( psTransform->padfGeoLocX == NULL ||
127 : psTransform->padfGeoLocY == NULL )
128 : {
129 : CPLError(CE_Failure, CPLE_OutOfMemory,
130 0 : "GeoLocLoadFullData : Out of memory");
131 0 : return FALSE;
132 : }
133 :
134 2 : if (nYSize_XBand == 1 && nYSize_YBand == 1)
135 : {
136 : /* Case of regular grid */
137 : /* The XBAND contains the x coordinates for all lines */
138 : /* The YBAND contains the y coordinates for all columns */
139 :
140 0 : double* padfTempX = (double*)VSIMalloc2(nXSize, sizeof(double));
141 0 : double* padfTempY = (double*)VSIMalloc2(nYSize, sizeof(double));
142 0 : if (padfTempX == NULL || padfTempY == NULL)
143 : {
144 0 : CPLFree(padfTempX);
145 0 : CPLFree(padfTempY);
146 : CPLError(CE_Failure, CPLE_OutOfMemory,
147 0 : "GeoLocLoadFullData : Out of memory");
148 0 : return FALSE;
149 : }
150 :
151 0 : CPLErr eErr = CE_None;
152 :
153 : eErr = GDALRasterIO( psTransform->hBand_X, GF_Read,
154 : 0, 0, nXSize, 1,
155 : padfTempX, nXSize, 1,
156 0 : GDT_Float64, 0, 0 );
157 :
158 : int i,j;
159 0 : for(j=0;j<nYSize;j++)
160 : {
161 : memcpy( psTransform->padfGeoLocX + j * nXSize,
162 : padfTempX,
163 0 : nXSize * sizeof(double) );
164 : }
165 :
166 0 : if (eErr == CE_None)
167 : {
168 : eErr = GDALRasterIO( psTransform->hBand_Y, GF_Read,
169 : 0, 0, nYSize, 1,
170 : padfTempY, nYSize, 1,
171 0 : GDT_Float64, 0, 0 );
172 :
173 0 : for(j=0;j<nYSize;j++)
174 : {
175 0 : for(i=0;i<nXSize;i++)
176 : {
177 0 : psTransform->padfGeoLocY[j * nXSize + i] = padfTempY[j];
178 : }
179 : }
180 : }
181 :
182 0 : CPLFree(padfTempX);
183 0 : CPLFree(padfTempY);
184 :
185 0 : if (eErr != CE_None)
186 0 : return FALSE;
187 : }
188 : else
189 : {
190 2 : if( GDALRasterIO( psTransform->hBand_X, GF_Read,
191 : 0, 0, nXSize, nYSize,
192 : psTransform->padfGeoLocX, nXSize, nYSize,
193 : GDT_Float64, 0, 0 ) != CE_None
194 : || GDALRasterIO( psTransform->hBand_Y, GF_Read,
195 : 0, 0, nXSize, nYSize,
196 : psTransform->padfGeoLocY, nXSize, nYSize,
197 : GDT_Float64, 0, 0 ) != CE_None )
198 0 : return FALSE;
199 : }
200 :
201 : psTransform->dfNoDataX = GDALGetRasterNoDataValue( psTransform->hBand_X,
202 2 : NULL );
203 : psTransform->dfNoDataY = GDALGetRasterNoDataValue( psTransform->hBand_Y,
204 2 : NULL );
205 :
206 2 : return TRUE;
207 : }
208 :
209 : /************************************************************************/
210 : /* GeoLocGenerateBackMap() */
211 : /************************************************************************/
212 :
213 2 : static int GeoLocGenerateBackMap( GDALGeoLocTransformInfo *psTransform )
214 :
215 : {
216 2 : int nXSize = psTransform->nGeoLocXSize;
217 2 : int nYSize = psTransform->nGeoLocYSize;
218 2 : int nMaxIter = 3;
219 :
220 : /* -------------------------------------------------------------------- */
221 : /* Scan forward map for lat/long extents. */
222 : /* -------------------------------------------------------------------- */
223 2 : double dfMinX=0, dfMaxX=0, dfMinY=0, dfMaxY=0;
224 2 : int i, bInit = FALSE;
225 :
226 4682 : for( i = nXSize * nYSize - 1; i >= 0; i-- )
227 : {
228 4680 : if( psTransform->padfGeoLocX[i] != psTransform->dfNoDataX )
229 : {
230 4680 : if( bInit )
231 : {
232 4678 : dfMinX = MIN(dfMinX,psTransform->padfGeoLocX[i]);
233 4678 : dfMaxX = MAX(dfMaxX,psTransform->padfGeoLocX[i]);
234 4678 : dfMinY = MIN(dfMinY,psTransform->padfGeoLocY[i]);
235 4678 : dfMaxY = MAX(dfMaxY,psTransform->padfGeoLocY[i]);
236 : }
237 : else
238 : {
239 2 : bInit = TRUE;
240 2 : dfMinX = dfMaxX = psTransform->padfGeoLocX[i];
241 2 : dfMinY = dfMaxY = psTransform->padfGeoLocY[i];
242 : }
243 : }
244 : }
245 :
246 : /* -------------------------------------------------------------------- */
247 : /* Decide on resolution for backmap. We aim for slightly */
248 : /* higher resolution than the source but we can't easily */
249 : /* establish how much dead space there is in the backmap, so it */
250 : /* is approximate. */
251 : /* -------------------------------------------------------------------- */
252 2 : double dfTargetPixels = (nXSize * nYSize * 1.3);
253 : double dfPixelSize = sqrt((dfMaxX - dfMinX) * (dfMaxY - dfMinY)
254 2 : / dfTargetPixels);
255 : int nBMXSize, nBMYSize;
256 :
257 : nBMYSize = psTransform->nBackMapHeight =
258 2 : (int) ((dfMaxY - dfMinY) / dfPixelSize + 1);
259 : nBMXSize= psTransform->nBackMapWidth =
260 2 : (int) ((dfMaxX - dfMinX) / dfPixelSize + 1);
261 :
262 2 : if (nBMXSize > INT_MAX / nBMYSize)
263 : {
264 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %d x %d",
265 0 : nBMXSize, nBMYSize);
266 0 : return FALSE;
267 : }
268 :
269 2 : dfMinX -= dfPixelSize/2.0;
270 2 : dfMaxY += dfPixelSize/2.0;
271 :
272 2 : psTransform->adfBackMapGeoTransform[0] = dfMinX;
273 2 : psTransform->adfBackMapGeoTransform[1] = dfPixelSize;
274 2 : psTransform->adfBackMapGeoTransform[2] = 0.0;
275 2 : psTransform->adfBackMapGeoTransform[3] = dfMaxY;
276 2 : psTransform->adfBackMapGeoTransform[4] = 0.0;
277 2 : psTransform->adfBackMapGeoTransform[5] = -dfPixelSize;
278 :
279 : /* -------------------------------------------------------------------- */
280 : /* Allocate backmap, and initialize to nodata value (-1.0). */
281 : /* -------------------------------------------------------------------- */
282 : GByte *pabyValidFlag;
283 :
284 : pabyValidFlag = (GByte *)
285 2 : VSICalloc(nBMXSize, nBMYSize);
286 :
287 : psTransform->pafBackMapX = (float *)
288 2 : VSIMalloc3(nBMXSize, nBMYSize, sizeof(float));
289 : psTransform->pafBackMapY = (float *)
290 2 : VSIMalloc3(nBMXSize, nBMYSize, sizeof(float));
291 :
292 2 : if( pabyValidFlag == NULL ||
293 : psTransform->pafBackMapX == NULL ||
294 : psTransform->pafBackMapY == NULL )
295 : {
296 : CPLError( CE_Failure, CPLE_OutOfMemory,
297 : "Unable to allocate %dx%d back-map for geolocation array transformer.",
298 0 : nBMXSize, nBMYSize );
299 0 : CPLFree( pabyValidFlag );
300 0 : return FALSE;
301 : }
302 :
303 6258 : for( i = nBMXSize * nBMYSize - 1; i >= 0; i-- )
304 : {
305 6256 : psTransform->pafBackMapX[i] = -1.0;
306 6256 : psTransform->pafBackMapY[i] = -1.0;
307 : }
308 :
309 : /* -------------------------------------------------------------------- */
310 : /* Run through the whole geoloc array forward projecting and */
311 : /* pushing into the backmap. */
312 : /* Initialise to the nMaxIter+1 value so we can spot genuinely */
313 : /* valid pixels in the hole-filling loop. */
314 : /* -------------------------------------------------------------------- */
315 : int iBMX, iBMY;
316 : int iX, iY;
317 :
318 80 : for( iY = 0; iY < nYSize; iY++ )
319 : {
320 4758 : for( iX = 0; iX < nXSize; iX++ )
321 : {
322 4680 : if( psTransform->padfGeoLocX[iX + iY * nXSize]
323 : == psTransform->dfNoDataX )
324 0 : continue;
325 :
326 4680 : i = iX + iY * nXSize;
327 :
328 4680 : iBMX = (int) ((psTransform->padfGeoLocX[i] - dfMinX) / dfPixelSize);
329 4680 : iBMY = (int) ((dfMaxY - psTransform->padfGeoLocY[i]) / dfPixelSize);
330 :
331 4680 : if( iBMX < 0 || iBMY < 0 || iBMX >= nBMXSize || iBMY >= nBMYSize )
332 2 : continue;
333 :
334 4678 : psTransform->pafBackMapX[iBMX + iBMY * nBMXSize] =
335 4678 : (float)(iX * psTransform->dfPIXEL_STEP + psTransform->dfPIXEL_OFFSET);
336 4678 : psTransform->pafBackMapY[iBMX + iBMY * nBMXSize] =
337 4678 : (float)(iY * psTransform->dfLINE_STEP + psTransform->dfLINE_OFFSET);
338 :
339 4678 : pabyValidFlag[iBMX + iBMY * nBMXSize] = (GByte) (nMaxIter+1);
340 :
341 : }
342 : }
343 :
344 : /* -------------------------------------------------------------------- */
345 : /* Now, loop over the backmap trying to fill in holes with */
346 : /* nearby values. */
347 : /* -------------------------------------------------------------------- */
348 : int iIter;
349 : int nNumValid;
350 :
351 8 : for( iIter = 0; iIter < nMaxIter; iIter++ )
352 : {
353 6 : nNumValid = 0;
354 282 : for( iBMY = 0; iBMY < nBMYSize; iBMY++ )
355 : {
356 19044 : for( iBMX = 0; iBMX < nBMXSize; iBMX++ )
357 : {
358 : // if this point is already set, ignore it.
359 18768 : if( pabyValidFlag[iBMX + iBMY*nBMXSize] )
360 : {
361 13078 : nNumValid++;
362 13078 : continue;
363 : }
364 :
365 5690 : int nCount = 0;
366 5690 : double dfXSum = 0.0, dfYSum = 0.0;
367 5690 : int nMarkedAsGood = nMaxIter - iIter;
368 :
369 : // left?
370 11144 : if( iBMX > 0 &&
371 5454 : pabyValidFlag[iBMX-1+iBMY*nBMXSize] > nMarkedAsGood )
372 : {
373 824 : dfXSum += psTransform->pafBackMapX[iBMX-1+iBMY*nBMXSize];
374 824 : dfYSum += psTransform->pafBackMapY[iBMX-1+iBMY*nBMXSize];
375 824 : nCount++;
376 : }
377 : // right?
378 11162 : if( iBMX + 1 < nBMXSize &&
379 5472 : pabyValidFlag[iBMX+1+iBMY*nBMXSize] > nMarkedAsGood )
380 : {
381 842 : dfXSum += psTransform->pafBackMapX[iBMX+1+iBMY*nBMXSize];
382 842 : dfYSum += psTransform->pafBackMapY[iBMX+1+iBMY*nBMXSize];
383 842 : nCount++;
384 : }
385 : // top?
386 11040 : if( iBMY > 0 &&
387 5350 : pabyValidFlag[iBMX+(iBMY-1)*nBMXSize] > nMarkedAsGood )
388 : {
389 796 : dfXSum += psTransform->pafBackMapX[iBMX+(iBMY-1)*nBMXSize];
390 796 : dfYSum += psTransform->pafBackMapY[iBMX+(iBMY-1)*nBMXSize];
391 796 : nCount++;
392 : }
393 : // bottom?
394 11024 : if( iBMY + 1 < nBMYSize &&
395 5334 : pabyValidFlag[iBMX+(iBMY+1)*nBMXSize] > nMarkedAsGood )
396 : {
397 780 : dfXSum += psTransform->pafBackMapX[iBMX+(iBMY+1)*nBMXSize];
398 780 : dfYSum += psTransform->pafBackMapY[iBMX+(iBMY+1)*nBMXSize];
399 780 : nCount++;
400 : }
401 : // top-left?
402 10810 : if( iBMX > 0 && iBMY > 0 &&
403 5120 : pabyValidFlag[iBMX-1+(iBMY-1)*nBMXSize] > nMarkedAsGood )
404 : {
405 966 : dfXSum += psTransform->pafBackMapX[iBMX-1+(iBMY-1)*nBMXSize];
406 966 : dfYSum += psTransform->pafBackMapY[iBMX-1+(iBMY-1)*nBMXSize];
407 966 : nCount++;
408 : }
409 : // top-right?
410 10828 : if( iBMX + 1 < nBMXSize && iBMY > 0 &&
411 5138 : pabyValidFlag[iBMX+1+(iBMY-1)*nBMXSize] > nMarkedAsGood )
412 : {
413 830 : dfXSum += psTransform->pafBackMapX[iBMX+1+(iBMY-1)*nBMXSize];
414 830 : dfYSum += psTransform->pafBackMapY[iBMX+1+(iBMY-1)*nBMXSize];
415 830 : nCount++;
416 : }
417 : // bottom-left?
418 10794 : if( iBMX > 0 && iBMY + 1 < nBMYSize &&
419 5104 : pabyValidFlag[iBMX-1+(iBMY+1)*nBMXSize] > nMarkedAsGood )
420 : {
421 796 : dfXSum += psTransform->pafBackMapX[iBMX-1+(iBMY+1)*nBMXSize];
422 796 : dfYSum += psTransform->pafBackMapY[iBMX-1+(iBMY+1)*nBMXSize];
423 796 : nCount++;
424 : }
425 : // bottom-right?
426 10812 : if( iBMX + 1 < nBMXSize && iBMY + 1 < nBMYSize &&
427 5122 : pabyValidFlag[iBMX+1+(iBMY+1)*nBMXSize] > nMarkedAsGood )
428 : {
429 968 : dfXSum += psTransform->pafBackMapX[iBMX+1+(iBMY+1)*nBMXSize];
430 968 : dfYSum += psTransform->pafBackMapY[iBMX+1+(iBMY+1)*nBMXSize];
431 968 : nCount++;
432 : }
433 :
434 5690 : if( nCount > 0 )
435 : {
436 1850 : psTransform->pafBackMapX[iBMX + iBMY * nBMXSize] = (float)(dfXSum/nCount);
437 1850 : psTransform->pafBackMapY[iBMX + iBMY * nBMXSize] = (float)(dfYSum/nCount);
438 : // genuinely valid points will have value iMaxIter+1
439 : // On each iteration mark newly valid points with a
440 : // descending value so that it will not be used on the
441 : // current iteration only on subsequent ones.
442 1850 : pabyValidFlag[iBMX+iBMY*nBMXSize] = (GByte) (nMaxIter - iIter);
443 : }
444 : }
445 : }
446 6 : if (nNumValid == nBMXSize * nBMYSize)
447 0 : break;
448 : }
449 :
450 : #ifdef notdef
451 : GDALDatasetH hBMDS = GDALCreate( GDALGetDriverByName( "GTiff" ),
452 : "backmap.tif", nBMXSize, nBMYSize, 2,
453 : GDT_Float32, NULL );
454 : GDALSetGeoTransform( hBMDS, psTransform->adfBackMapGeoTransform );
455 : GDALRasterIO( GDALGetRasterBand(hBMDS,1), GF_Write,
456 : 0, 0, nBMXSize, nBMYSize,
457 : psTransform->pafBackMapX, nBMXSize, nBMYSize,
458 : GDT_Float32, 0, 0 );
459 : GDALRasterIO( GDALGetRasterBand(hBMDS,2), GF_Write,
460 : 0, 0, nBMXSize, nBMYSize,
461 : psTransform->pafBackMapY, nBMXSize, nBMYSize,
462 : GDT_Float32, 0, 0 );
463 : GDALClose( hBMDS );
464 : #endif
465 :
466 2 : CPLFree( pabyValidFlag );
467 :
468 2 : return TRUE;
469 : }
470 :
471 : /************************************************************************/
472 : /* FindGeoLocPosition() */
473 : /************************************************************************/
474 :
475 : #ifdef notdef
476 :
477 : /*
478 : This searching approach has been abandoned because it is too sensitive
479 : to discontinuities in the data. Left in case it might be revived in
480 : the future.
481 : */
482 :
483 : static int FindGeoLocPosition( GDALGeoLocTransformInfo *psTransform,
484 : double dfGeoX, double dfGeoY,
485 : int nStartX, int nStartY,
486 : double *pdfFoundX, double *pdfFoundY )
487 :
488 : {
489 : double adfPathX[5000], adfPathY[5000];
490 :
491 : if( psTransform->padfGeoLocX == NULL )
492 : return FALSE;
493 :
494 : int nXSize = psTransform->nGeoLocXSize;
495 : int nYSize = psTransform->nGeoLocYSize;
496 : int nStepCount = 0;
497 :
498 : // Start in center if we don't have any provided info.
499 : if( nStartX < 0 || nStartY < 0
500 : || nStartX >= nXSize || nStartY >= nYSize )
501 : {
502 : nStartX = nXSize / 2;
503 : nStartY = nYSize / 2;
504 : }
505 :
506 : nStartX = MIN(nStartX,nXSize-2);
507 : nStartY = MIN(nStartY,nYSize-2);
508 :
509 : int iX = nStartX, iY = nStartY;
510 : int iLastX = -1, iLastY = -1;
511 : int iSecondLastX = -1, iSecondLastY = -1;
512 :
513 : while( nStepCount < MAX(nXSize,nYSize) )
514 : {
515 : int iXNext = -1, iYNext = -1;
516 : double dfDeltaXRight, dfDeltaYRight, dfDeltaXDown, dfDeltaYDown;
517 :
518 : double *padfThisX = psTransform->padfGeoLocX + iX + iY * nXSize;
519 : double *padfThisY = psTransform->padfGeoLocY + iX + iY * nXSize;
520 :
521 : double dfDeltaX = dfGeoX - *padfThisX;
522 : double dfDeltaY = dfGeoY - *padfThisY;
523 :
524 : if( iX == nXSize-1 )
525 : {
526 : dfDeltaXRight = *(padfThisX) - *(padfThisX-1);
527 : dfDeltaYRight = *(padfThisY) - *(padfThisY-1);
528 : }
529 : else
530 : {
531 : dfDeltaXRight = *(padfThisX+1) - *padfThisX;
532 : dfDeltaYRight = *(padfThisY+1) - *padfThisY;
533 : }
534 :
535 : if( iY == nYSize - 1 )
536 : {
537 : dfDeltaXDown = *(padfThisX) - *(padfThisX-nXSize);
538 : dfDeltaYDown = *(padfThisY) - *(padfThisY-nXSize);
539 : }
540 : else
541 : {
542 : dfDeltaXDown = *(padfThisX+nXSize) - *padfThisX;
543 : dfDeltaYDown = *(padfThisY+nXSize) - *padfThisY;
544 : }
545 :
546 : double dfRightProjection =
547 : (dfDeltaXRight * dfDeltaX + dfDeltaYRight * dfDeltaY)
548 : / (dfDeltaXRight*dfDeltaXRight + dfDeltaYRight*dfDeltaYRight);
549 :
550 : double dfDownProjection =
551 : (dfDeltaXDown * dfDeltaX + dfDeltaYDown * dfDeltaY)
552 : / (dfDeltaXDown*dfDeltaXDown + dfDeltaYDown*dfDeltaYDown);
553 :
554 : // Are we in our target cell?
555 : if( dfRightProjection >= 0.0 && dfRightProjection < 1.0
556 : && dfDownProjection >= 0.0 && dfDownProjection < 1.0 )
557 : {
558 : *pdfFoundX = iX + dfRightProjection;
559 : *pdfFoundY = iY + dfDownProjection;
560 :
561 : return TRUE;
562 : }
563 :
564 : if( ABS(dfRightProjection) > ABS(dfDownProjection) )
565 : {
566 : // Do we want to move right?
567 : if( dfRightProjection > 1.0 && iX < nXSize-1 )
568 : {
569 : iXNext = iX + MAX(1,(int)(dfRightProjection - nStepCount)/2);
570 : iYNext = iY;
571 : }
572 :
573 : // Do we want to move left?
574 : else if( dfRightProjection < 0.0 && iX > 0 )
575 : {
576 : iXNext = iX - MAX(1,(int)(ABS(dfRightProjection) - nStepCount)/2);
577 : iYNext = iY;
578 : }
579 :
580 : // Do we want to move down.
581 : else if( dfDownProjection > 1.0 && iY < nYSize-1 )
582 : {
583 : iXNext = iX;
584 : iYNext = iY + MAX(1,(int)(dfDownProjection - nStepCount)/2);
585 : }
586 :
587 : // Do we want to move up?
588 : else if( dfDownProjection < 0.0 && iY > 0 )
589 : {
590 : iXNext = iX;
591 : iYNext = iY - MAX(1,(int)(ABS(dfDownProjection) - nStepCount)/2);
592 : }
593 :
594 : // We aren't there, and we have no where to go
595 : else
596 : {
597 : return FALSE;
598 : }
599 : }
600 : else
601 : {
602 : // Do we want to move down.
603 : if( dfDownProjection > 1.0 && iY < nYSize-1 )
604 : {
605 : iXNext = iX;
606 : iYNext = iY + MAX(1,(int)(dfDownProjection - nStepCount)/2);
607 : }
608 :
609 : // Do we want to move up?
610 : else if( dfDownProjection < 0.0 && iY > 0 )
611 : {
612 : iXNext = iX;
613 : iYNext = iY - MAX(1,(int)(ABS(dfDownProjection) - nStepCount)/2);
614 : }
615 :
616 : // Do we want to move right?
617 : else if( dfRightProjection > 1.0 && iX < nXSize-1 )
618 : {
619 : iXNext = iX + MAX(1,(int)(dfRightProjection - nStepCount)/2);
620 : iYNext = iY;
621 : }
622 :
623 : // Do we want to move left?
624 : else if( dfRightProjection < 0.0 && iX > 0 )
625 : {
626 : iXNext = iX - MAX(1,(int)(ABS(dfRightProjection) - nStepCount)/2);
627 : iYNext = iY;
628 : }
629 :
630 : // We aren't there, and we have no where to go
631 : else
632 : {
633 : return FALSE;
634 : }
635 : }
636 : adfPathX[nStepCount] = iX;
637 : adfPathY[nStepCount] = iY;
638 :
639 : nStepCount++;
640 : iX = MAX(0,MIN(iXNext,nXSize-1));
641 : iY = MAX(0,MIN(iYNext,nYSize-1));
642 :
643 : if( iX == iSecondLastX && iY == iSecondLastY )
644 : {
645 : // Are we *near* our target cell?
646 : if( dfRightProjection >= -1.0 && dfRightProjection < 2.0
647 : && dfDownProjection >= -1.0 && dfDownProjection < 2.0 )
648 : {
649 : *pdfFoundX = iX + dfRightProjection;
650 : *pdfFoundY = iY + dfDownProjection;
651 :
652 : return TRUE;
653 : }
654 :
655 : #ifdef SHAPE_DEBUG
656 : if( hSHP != NULL )
657 : {
658 : SHPObject *hObj;
659 :
660 : hObj = SHPCreateSimpleObject( SHPT_ARC, nStepCount,
661 : adfPathX, adfPathY, NULL );
662 : SHPWriteObject( hSHP, -1, hObj );
663 : SHPDestroyObject( hObj );
664 :
665 : int iShape = DBFGetRecordCount( hDBF );
666 : DBFWriteDoubleAttribute( hDBF, iShape, 0, dfGeoX );
667 : DBFWriteDoubleAttribute( hDBF, iShape, 1, dfGeoY );
668 : }
669 : #endif
670 : //CPLDebug( "GeoL", "Looping at step (%d) on search for %g,%g.",
671 : // nStepCount, dfGeoX, dfGeoY );
672 : return FALSE;
673 : }
674 :
675 : iSecondLastX = iLastX;
676 : iSecondLastY = iLastY;
677 :
678 : iLastX = iX;
679 : iLastY = iY;
680 :
681 : }
682 :
683 : //CPLDebug( "GeoL", "Exceeded step count max (%d) on search for %g,%g.",
684 : // MAX(nXSize,nYSize),
685 : // dfGeoX, dfGeoY );
686 :
687 : #ifdef SHAPE_DEBUG
688 : if( hSHP != NULL )
689 : {
690 : SHPObject *hObj;
691 :
692 : hObj = SHPCreateSimpleObject( SHPT_ARC, nStepCount,
693 : adfPathX, adfPathY, NULL );
694 : SHPWriteObject( hSHP, -1, hObj );
695 : SHPDestroyObject( hObj );
696 :
697 : int iShape = DBFGetRecordCount( hDBF );
698 : DBFWriteDoubleAttribute( hDBF, iShape, 0, dfGeoX );
699 : DBFWriteDoubleAttribute( hDBF, iShape, 1, dfGeoY );
700 : }
701 : #endif
702 :
703 : return FALSE;
704 : }
705 : #endif /* def notdef */
706 :
707 :
708 :
709 : /************************************************************************/
710 : /* GDALCreateGeoLocTransformer() */
711 : /************************************************************************/
712 :
713 2 : void *GDALCreateGeoLocTransformer( GDALDatasetH hBaseDS,
714 : char **papszGeolocationInfo,
715 : int bReversed )
716 :
717 : {
718 : GDALGeoLocTransformInfo *psTransform;
719 :
720 2 : if( CSLFetchNameValue(papszGeolocationInfo,"PIXEL_OFFSET") == NULL
721 : || CSLFetchNameValue(papszGeolocationInfo,"LINE_OFFSET") == NULL
722 : || CSLFetchNameValue(papszGeolocationInfo,"PIXEL_STEP") == NULL
723 : || CSLFetchNameValue(papszGeolocationInfo,"LINE_STEP") == NULL
724 : || CSLFetchNameValue(papszGeolocationInfo,"X_BAND") == NULL
725 : || CSLFetchNameValue(papszGeolocationInfo,"Y_BAND") == NULL )
726 : {
727 : CPLError( CE_Failure, CPLE_AppDefined,
728 0 : "Missing some geolocation fields in GDALCreateGeoLocTransformer()" );
729 0 : return NULL;
730 : }
731 :
732 : /* -------------------------------------------------------------------- */
733 : /* Initialize core info. */
734 : /* -------------------------------------------------------------------- */
735 : psTransform = (GDALGeoLocTransformInfo *)
736 2 : CPLCalloc(sizeof(GDALGeoLocTransformInfo),1);
737 :
738 2 : psTransform->bReversed = bReversed;
739 :
740 2 : strcpy( psTransform->sTI.szSignature, "GTI" );
741 2 : psTransform->sTI.pszClassName = "GDALGeoLocTransformer";
742 2 : psTransform->sTI.pfnTransform = GDALGeoLocTransform;
743 2 : psTransform->sTI.pfnCleanup = GDALDestroyGeoLocTransformer;
744 2 : psTransform->sTI.pfnSerialize = GDALSerializeGeoLocTransformer;
745 2 : psTransform->papszGeolocationInfo = CSLDuplicate( papszGeolocationInfo );
746 :
747 : /* -------------------------------------------------------------------- */
748 : /* Pull geolocation info from the options/metadata. */
749 : /* -------------------------------------------------------------------- */
750 : psTransform->dfPIXEL_OFFSET = atof(CSLFetchNameValue( papszGeolocationInfo,
751 2 : "PIXEL_OFFSET" ));
752 : psTransform->dfLINE_OFFSET = atof(CSLFetchNameValue( papszGeolocationInfo,
753 2 : "LINE_OFFSET" ));
754 : psTransform->dfPIXEL_STEP = atof(CSLFetchNameValue( papszGeolocationInfo,
755 2 : "PIXEL_STEP" ));
756 : psTransform->dfLINE_STEP = atof(CSLFetchNameValue( papszGeolocationInfo,
757 2 : "LINE_STEP" ));
758 :
759 : /* -------------------------------------------------------------------- */
760 : /* Establish access to geolocation dataset(s). */
761 : /* -------------------------------------------------------------------- */
762 : const char *pszDSName = CSLFetchNameValue( papszGeolocationInfo,
763 2 : "X_DATASET" );
764 2 : if( pszDSName != NULL )
765 : {
766 2 : psTransform->hDS_X = GDALOpenShared( pszDSName, GA_ReadOnly );
767 : }
768 : else
769 : {
770 0 : psTransform->hDS_X = hBaseDS;
771 0 : GDALReferenceDataset( psTransform->hDS_X );
772 : psTransform->papszGeolocationInfo =
773 : CSLSetNameValue( psTransform->papszGeolocationInfo,
774 : "X_DATASET",
775 0 : GDALGetDescription( hBaseDS ) );
776 : }
777 :
778 2 : pszDSName = CSLFetchNameValue( papszGeolocationInfo, "Y_DATASET" );
779 2 : if( pszDSName != NULL )
780 : {
781 2 : psTransform->hDS_Y = GDALOpenShared( pszDSName, GA_ReadOnly );
782 : }
783 : else
784 : {
785 0 : psTransform->hDS_Y = hBaseDS;
786 0 : GDALReferenceDataset( psTransform->hDS_Y );
787 : psTransform->papszGeolocationInfo =
788 : CSLSetNameValue( psTransform->papszGeolocationInfo,
789 : "Y_DATASET",
790 0 : GDALGetDescription( hBaseDS ) );
791 : }
792 :
793 2 : if (psTransform->hDS_X == NULL ||
794 : psTransform->hDS_Y == NULL)
795 : {
796 0 : GDALDestroyGeoLocTransformer( psTransform );
797 0 : return NULL;
798 : }
799 :
800 : /* -------------------------------------------------------------------- */
801 : /* Get the band handles. */
802 : /* -------------------------------------------------------------------- */
803 : int nBand;
804 :
805 2 : nBand = MAX(1,atoi(CSLFetchNameValue( papszGeolocationInfo, "X_BAND" )));
806 2 : psTransform->hBand_X = GDALGetRasterBand( psTransform->hDS_X, nBand );
807 :
808 2 : nBand = MAX(1,atoi(CSLFetchNameValue( papszGeolocationInfo, "Y_BAND" )));
809 2 : psTransform->hBand_Y = GDALGetRasterBand( psTransform->hDS_Y, nBand );
810 :
811 2 : if (psTransform->hBand_X == NULL ||
812 : psTransform->hBand_Y == NULL)
813 : {
814 0 : GDALDestroyGeoLocTransformer( psTransform );
815 0 : return NULL;
816 : }
817 :
818 : /* -------------------------------------------------------------------- */
819 : /* Check that X and Y bands have the same dimensions */
820 : /* -------------------------------------------------------------------- */
821 2 : int nXSize_XBand = GDALGetRasterXSize( psTransform->hDS_X );
822 2 : int nYSize_XBand = GDALGetRasterYSize( psTransform->hDS_X );
823 2 : int nXSize_YBand = GDALGetRasterXSize( psTransform->hDS_Y );
824 2 : int nYSize_YBand = GDALGetRasterYSize( psTransform->hDS_Y );
825 2 : if (nYSize_XBand == 1 || nYSize_YBand == 1)
826 : {
827 0 : if (nYSize_XBand != 1 || nYSize_YBand != 1)
828 : {
829 : CPLError(CE_Failure, CPLE_AppDefined,
830 0 : "X_BAND and Y_BAND should have both nYSize == 1");
831 0 : GDALDestroyGeoLocTransformer( psTransform );
832 0 : return NULL;
833 : }
834 : }
835 2 : else if (nXSize_XBand != nXSize_YBand ||
836 : nYSize_XBand != nYSize_YBand )
837 : {
838 : CPLError(CE_Failure, CPLE_AppDefined,
839 0 : "X_BAND and Y_BAND do not have the same dimensions");
840 0 : GDALDestroyGeoLocTransformer( psTransform );
841 0 : return NULL;
842 : }
843 :
844 2 : if (nXSize_XBand > INT_MAX / nYSize_XBand)
845 : {
846 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow : %d x %d",
847 0 : nXSize_XBand, nYSize_XBand);
848 0 : GDALDestroyGeoLocTransformer( psTransform );
849 0 : return NULL;
850 : }
851 :
852 : /* -------------------------------------------------------------------- */
853 : /* Load the geolocation array. */
854 : /* -------------------------------------------------------------------- */
855 2 : if( !GeoLocLoadFullData( psTransform )
856 : || !GeoLocGenerateBackMap( psTransform ) )
857 : {
858 0 : GDALDestroyGeoLocTransformer( psTransform );
859 0 : return NULL;
860 : }
861 :
862 2 : return psTransform;
863 : }
864 :
865 : /************************************************************************/
866 : /* GDALDestroyGeoLocTransformer() */
867 : /************************************************************************/
868 :
869 2 : void GDALDestroyGeoLocTransformer( void *pTransformAlg )
870 :
871 : {
872 : GDALGeoLocTransformInfo *psTransform =
873 2 : (GDALGeoLocTransformInfo *) pTransformAlg;
874 :
875 2 : CPLFree( psTransform->pafBackMapX );
876 2 : CPLFree( psTransform->pafBackMapY );
877 2 : CSLDestroy( psTransform->papszGeolocationInfo );
878 2 : CPLFree( psTransform->padfGeoLocX );
879 2 : CPLFree( psTransform->padfGeoLocY );
880 :
881 2 : if( psTransform->hDS_X != NULL
882 : && GDALDereferenceDataset( psTransform->hDS_X ) == 0 )
883 0 : GDALClose( psTransform->hDS_X );
884 :
885 2 : if( psTransform->hDS_Y != NULL
886 : && GDALDereferenceDataset( psTransform->hDS_Y ) == 0 )
887 2 : GDALClose( psTransform->hDS_Y );
888 :
889 2 : CPLFree( pTransformAlg );
890 2 : }
891 :
892 : /************************************************************************/
893 : /* GDALGeoLocTransform() */
894 : /************************************************************************/
895 :
896 260 : int GDALGeoLocTransform( void *pTransformArg, int bDstToSrc,
897 : int nPointCount,
898 : double *padfX, double *padfY, double *padfZ,
899 : int *panSuccess )
900 :
901 : {
902 : GDALGeoLocTransformInfo *psTransform =
903 260 : (GDALGeoLocTransformInfo *) pTransformArg;
904 :
905 260 : if( psTransform->bReversed )
906 0 : bDstToSrc = !bDstToSrc;
907 :
908 : /* -------------------------------------------------------------------- */
909 : /* Do original pixel line to target geox/geoy. */
910 : /* -------------------------------------------------------------------- */
911 260 : if( !bDstToSrc )
912 : {
913 1 : int i, nXSize = psTransform->nGeoLocXSize;
914 :
915 2 : for( i = 0; i < nPointCount; i++ )
916 : {
917 1 : if( padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL )
918 : {
919 0 : panSuccess[i] = FALSE;
920 0 : continue;
921 : }
922 :
923 1 : double dfGeoLocPixel = (padfX[i] - psTransform->dfPIXEL_OFFSET)
924 1 : / psTransform->dfPIXEL_STEP;
925 1 : double dfGeoLocLine = (padfY[i] - psTransform->dfLINE_OFFSET)
926 1 : / psTransform->dfLINE_STEP;
927 :
928 : int iX, iY;
929 :
930 1 : iX = MAX(0,(int) dfGeoLocPixel);
931 1 : iX = MIN(iX,psTransform->nGeoLocXSize-1);
932 1 : iY = MAX(0,(int) dfGeoLocLine);
933 1 : iY = MIN(iY,psTransform->nGeoLocYSize-1);
934 :
935 1 : double *padfGLX = psTransform->padfGeoLocX + iX + iY * nXSize;
936 1 : double *padfGLY = psTransform->padfGeoLocY + iX + iY * nXSize;
937 :
938 : // This assumes infinite extension beyond borders of available
939 : // data based on closest grid square.
940 :
941 2 : if( iX + 1 < psTransform->nGeoLocXSize &&
942 : iY + 1 < psTransform->nGeoLocYSize )
943 : {
944 1 : padfX[i] = padfGLX[0]
945 1 : + (dfGeoLocPixel-iX) * (padfGLX[1] - padfGLX[0])
946 2 : + (dfGeoLocLine -iY) * (padfGLX[nXSize] - padfGLX[0]);
947 1 : padfY[i] = padfGLY[0]
948 1 : + (dfGeoLocPixel-iX) * (padfGLY[1] - padfGLY[0])
949 2 : + (dfGeoLocLine -iY) * (padfGLY[nXSize] - padfGLY[0]);
950 : }
951 0 : else if( iX + 1 < psTransform->nGeoLocXSize )
952 : {
953 0 : padfX[i] = padfGLX[0]
954 0 : + (dfGeoLocPixel-iX) * (padfGLX[1] - padfGLX[0]);
955 0 : padfY[i] = padfGLY[0]
956 0 : + (dfGeoLocPixel-iX) * (padfGLY[1] - padfGLY[0]);
957 : }
958 0 : else if( iY + 1 < psTransform->nGeoLocYSize )
959 : {
960 0 : padfX[i] = padfGLX[0]
961 0 : + (dfGeoLocLine -iY) * (padfGLX[nXSize] - padfGLX[0]);
962 0 : padfY[i] = padfGLY[0]
963 0 : + (dfGeoLocLine -iY) * (padfGLY[nXSize] - padfGLY[0]);
964 : }
965 : else
966 : {
967 0 : padfX[i] = padfGLX[0];
968 0 : padfY[i] = padfGLY[0];
969 : }
970 :
971 1 : panSuccess[i] = TRUE;
972 : }
973 : }
974 :
975 : /* -------------------------------------------------------------------- */
976 : /* geox/geoy to pixel/line using backmap. */
977 : /* -------------------------------------------------------------------- */
978 : else
979 : {
980 : int i;
981 :
982 66705 : for( i = 0; i < nPointCount; i++ )
983 : {
984 66446 : if( padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL )
985 : {
986 0 : panSuccess[i] = FALSE;
987 0 : continue;
988 : }
989 :
990 : double dfBMX, dfBMY;
991 : int iBMX, iBMY;
992 :
993 66446 : dfBMX = ((padfX[i] - psTransform->adfBackMapGeoTransform[0])
994 66446 : / psTransform->adfBackMapGeoTransform[1]);
995 66446 : dfBMY = ((padfY[i] - psTransform->adfBackMapGeoTransform[3])
996 66446 : / psTransform->adfBackMapGeoTransform[5]);
997 :
998 66446 : iBMX = (int) dfBMX;
999 66446 : iBMY = (int) dfBMY;
1000 :
1001 66446 : int iBM = iBMX + iBMY * psTransform->nBackMapWidth;
1002 :
1003 69480 : if( iBMX < 0 || iBMY < 0
1004 : || iBMX >= psTransform->nBackMapWidth
1005 : || iBMY >= psTransform->nBackMapHeight
1006 3034 : || psTransform->pafBackMapX[iBM] < 0 )
1007 : {
1008 63879 : panSuccess[i] = FALSE;
1009 63879 : padfX[i] = HUGE_VAL;
1010 63879 : padfY[i] = HUGE_VAL;
1011 63879 : continue;
1012 : }
1013 :
1014 2567 : float* pafBMX = psTransform->pafBackMapX + iBM;
1015 2567 : float* pafBMY = psTransform->pafBackMapY + iBM;
1016 10055 : if( iBMX + 1 < psTransform->nBackMapWidth &&
1017 : iBMY + 1 < psTransform->nBackMapHeight &&
1018 5021 : pafBMX[1] >=0 && pafBMX[psTransform->nBackMapWidth] >= 0 )
1019 : {
1020 2467 : padfX[i] = pafBMX[0] +
1021 2467 : (dfBMX - iBMX) * (pafBMX[1] - pafBMX[0]) +
1022 4934 : (dfBMY - iBMY) * (pafBMX[psTransform->nBackMapWidth] - pafBMX[0]);
1023 2467 : padfY[i] = pafBMY[0] +
1024 2467 : (dfBMX - iBMX) * (pafBMY[1] - pafBMY[0]) +
1025 4934 : (dfBMY - iBMY) * (pafBMY[psTransform->nBackMapWidth] - pafBMY[0]);
1026 : }
1027 234 : else if( iBMX + 1 < psTransform->nBackMapWidth &&
1028 79 : pafBMX[1] >=0)
1029 : {
1030 55 : padfX[i] = pafBMX[0] +
1031 55 : (dfBMX - iBMX) * (pafBMX[1] - pafBMX[0]);
1032 55 : padfY[i] = pafBMY[0] +
1033 55 : (dfBMX - iBMX) * (pafBMY[1] - pafBMY[0]);
1034 : }
1035 127 : else if( iBMY + 1 < psTransform->nBackMapHeight &&
1036 44 : pafBMX[psTransform->nBackMapWidth] >= 0 )
1037 : {
1038 38 : padfX[i] = pafBMX[0] +
1039 38 : (dfBMY - iBMY) * (pafBMX[psTransform->nBackMapWidth] - pafBMX[0]);
1040 38 : padfY[i] = pafBMY[0] +
1041 38 : (dfBMY - iBMY) * (pafBMY[psTransform->nBackMapWidth] - pafBMY[0]);
1042 : }
1043 : else
1044 : {
1045 7 : padfX[i] = pafBMX[0];
1046 7 : padfY[i] = pafBMY[0];
1047 : }
1048 2567 : panSuccess[i] = TRUE;
1049 : }
1050 : }
1051 :
1052 : /* -------------------------------------------------------------------- */
1053 : /* geox/geoy to pixel/line using search algorithm. */
1054 : /* -------------------------------------------------------------------- */
1055 : #ifdef notdef
1056 : else
1057 : {
1058 : int i;
1059 : int nStartX = -1, nStartY = -1;
1060 :
1061 : #ifdef SHAPE_DEBUG
1062 : hSHP = SHPCreate( "tracks.shp", SHPT_ARC );
1063 : hDBF = DBFCreate( "tracks.dbf" );
1064 : DBFAddField( hDBF, "GEOX", FTDouble, 10, 4 );
1065 : DBFAddField( hDBF, "GEOY", FTDouble, 10, 4 );
1066 : #endif
1067 : for( i = 0; i < nPointCount; i++ )
1068 : {
1069 : double dfGeoLocX, dfGeoLocY;
1070 :
1071 : if( padfX[i] == HUGE_VAL || padfY[i] == HUGE_VAL )
1072 : {
1073 : panSuccess[i] = FALSE;
1074 : continue;
1075 : }
1076 :
1077 : if( !FindGeoLocPosition( psTransform, padfX[i], padfY[i],
1078 : -1, -1, &dfGeoLocX, &dfGeoLocY ) )
1079 : {
1080 : padfX[i] = HUGE_VAL;
1081 : padfY[i] = HUGE_VAL;
1082 :
1083 : panSuccess[i] = FALSE;
1084 : continue;
1085 : }
1086 : nStartX = (int) dfGeoLocX;
1087 : nStartY = (int) dfGeoLocY;
1088 :
1089 : padfX[i] = dfGeoLocX * psTransform->dfPIXEL_STEP
1090 : + psTransform->dfPIXEL_OFFSET;
1091 : padfY[i] = dfGeoLocY * psTransform->dfLINE_STEP
1092 : + psTransform->dfLINE_OFFSET;
1093 :
1094 : panSuccess[i] = TRUE;
1095 : }
1096 :
1097 : #ifdef SHAPE_DEBUG
1098 : if( hSHP != NULL )
1099 : {
1100 : DBFClose( hDBF );
1101 : hDBF = NULL;
1102 :
1103 : SHPClose( hSHP );
1104 : hSHP = NULL;
1105 : }
1106 : #endif
1107 : }
1108 : #endif
1109 :
1110 260 : return TRUE;
1111 : }
1112 :
1113 : /************************************************************************/
1114 : /* GDALSerializeGeoLocTransformer() */
1115 : /************************************************************************/
1116 :
1117 0 : CPLXMLNode *GDALSerializeGeoLocTransformer( void *pTransformArg )
1118 :
1119 : {
1120 0 : VALIDATE_POINTER1( pTransformArg, "GDALSerializeGeoLocTransformer", NULL );
1121 :
1122 : CPLXMLNode *psTree;
1123 : GDALGeoLocTransformInfo *psInfo =
1124 0 : (GDALGeoLocTransformInfo *)(pTransformArg);
1125 :
1126 0 : psTree = CPLCreateXMLNode( NULL, CXT_Element, "GeoLocTransformer" );
1127 :
1128 : /* -------------------------------------------------------------------- */
1129 : /* Serialize bReversed. */
1130 : /* -------------------------------------------------------------------- */
1131 : CPLCreateXMLElementAndValue(
1132 : psTree, "Reversed",
1133 0 : CPLString().Printf( "%d", psInfo->bReversed ) );
1134 :
1135 : /* -------------------------------------------------------------------- */
1136 : /* geoloc metadata. */
1137 : /* -------------------------------------------------------------------- */
1138 0 : char **papszMD = psInfo->papszGeolocationInfo;
1139 : CPLXMLNode *psMD= CPLCreateXMLNode( psTree, CXT_Element,
1140 0 : "Metadata" );
1141 :
1142 0 : for( int i = 0; papszMD != NULL && papszMD[i] != NULL; i++ )
1143 : {
1144 : const char *pszRawValue;
1145 : char *pszKey;
1146 : CPLXMLNode *psMDI;
1147 :
1148 0 : pszRawValue = CPLParseNameValue( papszMD[i], &pszKey );
1149 :
1150 0 : psMDI = CPLCreateXMLNode( psMD, CXT_Element, "MDI" );
1151 0 : CPLSetXMLValue( psMDI, "#key", pszKey );
1152 0 : CPLCreateXMLNode( psMDI, CXT_Text, pszRawValue );
1153 :
1154 0 : CPLFree( pszKey );
1155 : }
1156 :
1157 0 : return psTree;
1158 : }
1159 :
1160 : /************************************************************************/
1161 : /* GDALDeserializeGeoLocTransformer() */
1162 : /************************************************************************/
1163 :
1164 1 : void *GDALDeserializeGeoLocTransformer( CPLXMLNode *psTree )
1165 :
1166 : {
1167 : void *pResult;
1168 : int bReversed;
1169 1 : char **papszMD = NULL;
1170 : CPLXMLNode *psMDI, *psMetadata;
1171 :
1172 : /* -------------------------------------------------------------------- */
1173 : /* Collect metadata. */
1174 : /* -------------------------------------------------------------------- */
1175 1 : psMetadata = CPLGetXMLNode( psTree, "Metadata" );
1176 :
1177 1 : if( psMetadata == NULL ||
1178 : psMetadata->eType != CXT_Element
1179 : || !EQUAL(psMetadata->pszValue,"Metadata") )
1180 0 : return NULL;
1181 :
1182 10 : for( psMDI = psMetadata->psChild; psMDI != NULL;
1183 : psMDI = psMDI->psNext )
1184 : {
1185 9 : if( !EQUAL(psMDI->pszValue,"MDI")
1186 : || psMDI->eType != CXT_Element
1187 : || psMDI->psChild == NULL
1188 : || psMDI->psChild->psNext == NULL
1189 : || psMDI->psChild->eType != CXT_Attribute
1190 : || psMDI->psChild->psChild == NULL )
1191 0 : continue;
1192 :
1193 : papszMD =
1194 : CSLSetNameValue( papszMD,
1195 : psMDI->psChild->psChild->pszValue,
1196 9 : psMDI->psChild->psNext->pszValue );
1197 : }
1198 :
1199 : /* -------------------------------------------------------------------- */
1200 : /* Get other flags. */
1201 : /* -------------------------------------------------------------------- */
1202 1 : bReversed = atoi(CPLGetXMLValue(psTree,"Reversed","0"));
1203 :
1204 : /* -------------------------------------------------------------------- */
1205 : /* Generate transformation. */
1206 : /* -------------------------------------------------------------------- */
1207 1 : pResult = GDALCreateGeoLocTransformer( NULL, papszMD, bReversed );
1208 :
1209 : /* -------------------------------------------------------------------- */
1210 : /* Cleanup GCP copy. */
1211 : /* -------------------------------------------------------------------- */
1212 1 : CSLDestroy( papszMD );
1213 :
1214 1 : return pResult;
1215 : }
|
