1 : /******************************************************************************
2 : * $Id: gdalgrid.cpp 18041 2009-11-17 14:14:43Z dron $
3 : *
4 : * Project: GDAL Gridding API.
5 : * Purpose: Implementation of GDAL scattered data gridder.
6 : * Author: Andrey Kiselev, dron@ak4719.spb.edu
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2007, Andrey Kiselev <dron@ak4719.spb.edu>
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 "cpl_vsi.h"
31 : #include "gdalgrid.h"
32 :
33 : CPL_CVSID("$Id: gdalgrid.cpp 18041 2009-11-17 14:14:43Z dron $");
34 :
35 : #define TO_RADIANS (3.14159265358979323846 / 180.0)
36 :
37 : /************************************************************************/
38 : /* GDALGridInverseDistanceToAPower() */
39 : /************************************************************************/
40 :
41 : /**
42 : * Inverse distance to a power.
43 : *
44 : * The Inverse Distance to a Power gridding method is a weighted average
45 : * interpolator. You should supply the input arrays with the scattered data
46 : * values including coordinates of every data point and output grid geometry.
47 : * The function will compute interpolated value for the given position in
48 : * output grid.
49 : *
50 : * For every grid node the resulting value \f$Z\f$ will be calculated using
51 : * formula:
52 : *
53 : * \f[
54 : * Z=\frac{\sum_{i=1}^n{\frac{Z_i}{r_i^p}}}{\sum_{i=1}^n{\frac{1}{r_i^p}}}
55 : * \f]
56 : *
57 : * where
58 : * <ul>
59 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
60 : * <li> \f$r_i\f$ is an Euclidean distance from the grid node
61 : * to point \f$i\f$,
62 : * <li> \f$p\f$ is a weighting power,
63 : * <li> \f$n\f$ is a total number of points in search ellipse.
64 : * </ul>
65 : *
66 : * In this method the weighting factor \f$w\f$ is
67 : *
68 : * \f[
69 : * w=\frac{1}{r^p}
70 : * \f]
71 : *
72 : * @param poOptions Algorithm parameters. This should point to
73 : * GDALGridInverseDistanceToAPowerOptions object.
74 : * @param nPoints Number of elements in input arrays.
75 : * @param padfX Input array of X coordinates.
76 : * @param padfY Input array of Y coordinates.
77 : * @param padfZ Input array of Z values.
78 : * @param dfXPoint X coordinate of the point to compute.
79 : * @param dfYPoint Y coordinate of the point to compute.
80 : * @param nXPoint X position of the point to compute.
81 : * @param nYPoint Y position of the point to compute.
82 : * @param pdfValue Pointer to variable where the computed grid node value
83 : * will be returned.
84 : *
85 : * @return CE_None on success or CE_Failure if something goes wrong.
86 : */
87 :
88 : CPLErr
89 0 : GDALGridInverseDistanceToAPower( const void *poOptions, GUInt32 nPoints,
90 : const double *padfX, const double *padfY,
91 : const double *padfZ,
92 : double dfXPoint, double dfYPoint,
93 : double *pdfValue )
94 : {
95 : // TODO: For optimization purposes pre-computed parameters should be moved
96 : // out of this routine to the calling function.
97 :
98 : // Pre-compute search ellipse parameters
99 : double dfRadius1 =
100 0 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfRadius1;
101 : double dfRadius2 =
102 0 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfRadius2;
103 : double dfR12;
104 :
105 0 : dfRadius1 *= dfRadius1;
106 0 : dfRadius2 *= dfRadius2;
107 0 : dfR12 = dfRadius1 * dfRadius2;
108 :
109 : // Compute coefficients for coordinate system rotation.
110 0 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
111 : const double dfAngle = TO_RADIANS
112 0 : * ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfAngle;
113 0 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
114 0 : if ( bRotated )
115 : {
116 0 : dfCoeff1 = cos(dfAngle);
117 0 : dfCoeff2 = sin(dfAngle);
118 : }
119 :
120 : const double dfPower =
121 0 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfPower;
122 : const double dfSmoothing =
123 0 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfSmoothing;
124 : const GUInt32 nMaxPoints =
125 0 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->nMaxPoints;
126 0 : double dfNominator = 0.0, dfDenominator = 0.0;
127 0 : GUInt32 i, n = 0;
128 :
129 0 : for ( i = 0; i < nPoints; i++ )
130 : {
131 0 : double dfRX = padfX[i] - dfXPoint;
132 0 : double dfRY = padfY[i] - dfYPoint;
133 : const double dfR2 =
134 0 : dfRX * dfRX + dfRY * dfRY + dfSmoothing * dfSmoothing;
135 :
136 0 : if ( bRotated )
137 : {
138 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
139 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
140 :
141 0 : dfRX = dfRXRotated;
142 0 : dfRY = dfRYRotated;
143 : }
144 :
145 : // Is this point located inside the search ellipse?
146 0 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
147 : {
148 0 : if ( CPLIsEqual(dfR2, 0.0) )
149 : {
150 0 : (*pdfValue) = padfZ[i];
151 0 : return CE_None;
152 : }
153 : else
154 : {
155 0 : const double dfW = pow( sqrt(dfR2), dfPower );
156 0 : dfNominator += padfZ[i] / dfW;
157 0 : dfDenominator += 1.0 / dfW;
158 0 : n++;
159 0 : if ( nMaxPoints > 0 && n > nMaxPoints )
160 0 : break;
161 : }
162 : }
163 : }
164 :
165 0 : if ( n < ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->nMinPoints
166 : || dfDenominator == 0.0 )
167 : {
168 : (*pdfValue) =
169 0 : ((GDALGridInverseDistanceToAPowerOptions*)poOptions)->dfNoDataValue;
170 : }
171 : else
172 0 : (*pdfValue) = dfNominator / dfDenominator;
173 :
174 0 : return CE_None;
175 : }
176 :
177 : /************************************************************************/
178 : /* GDALGridInverseDistanceToAPowerNoSearch() */
179 : /************************************************************************/
180 :
181 : /**
182 : * Inverse distance to a power for whole data set.
183 : *
184 : * This is somewhat optimized version of the Inverse Distance to a Power
185 : * method. It is used when the search ellips is not set. The algorithm and
186 : * parameters are the same as in GDALGridInverseDistanceToAPower(), but this
187 : * implementation works faster, because of no search.
188 : *
189 : * @see GDALGridInverseDistanceToAPower()
190 : */
191 :
192 : CPLErr
193 0 : GDALGridInverseDistanceToAPowerNoSearch( const void *poOptions, GUInt32 nPoints,
194 : const double *padfX, const double *padfY,
195 : const double *padfZ,
196 : double dfXPoint, double dfYPoint,
197 : double *pdfValue )
198 : {
199 : const double dfPower =
200 0 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfPower;
201 : const double dfSmoothing =
202 0 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->dfSmoothing;
203 0 : double dfNominator = 0.0, dfDenominator = 0.0;
204 : GUInt32 i;
205 :
206 0 : for ( i = 0; i < nPoints; i++ )
207 : {
208 0 : const double dfRX = padfX[i] - dfXPoint;
209 0 : const double dfRY = padfY[i] - dfYPoint;
210 : const double dfR2 =
211 0 : dfRX * dfRX + dfRY * dfRY + dfSmoothing * dfSmoothing;
212 :
213 0 : if ( CPLIsEqual(dfR2, 0.0) )
214 : {
215 0 : (*pdfValue) = padfZ[i];
216 0 : return CE_None;
217 : }
218 : else
219 : {
220 0 : const double dfW = pow( sqrt(dfR2), dfPower );
221 0 : dfNominator += padfZ[i] / dfW;
222 0 : dfDenominator += 1.0 / dfW;
223 : }
224 : }
225 :
226 0 : if ( dfDenominator == 0.0 )
227 : {
228 : (*pdfValue) =
229 0 : ((GDALGridInverseDistanceToAPowerOptions*)poOptions)->dfNoDataValue;
230 : }
231 : else
232 0 : (*pdfValue) = dfNominator / dfDenominator;
233 :
234 0 : return CE_None;
235 : }
236 : /************************************************************************/
237 : /* GDALGridMovingAverage() */
238 : /************************************************************************/
239 :
240 : /**
241 : * Moving average.
242 : *
243 : * The Moving Average is a simple data averaging algorithm. It uses a moving
244 : * window of elliptic form to search values and averages all data points
245 : * within the window. Search ellipse can be rotated by specified angle, the
246 : * center of ellipse located at the grid node. Also the minimum number of data
247 : * points to average can be set, if there are not enough points in window, the
248 : * grid node considered empty and will be filled with specified NODATA value.
249 : *
250 : * Mathematically it can be expressed with the formula:
251 : *
252 : * \f[
253 : * Z=\frac{\sum_{i=1}^n{Z_i}}{n}
254 : * \f]
255 : *
256 : * where
257 : * <ul>
258 : * <li> \f$Z\f$ is a resulting value at the grid node,
259 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
260 : * <li> \f$n\f$ is a total number of points in search ellipse.
261 : * </ul>
262 : *
263 : * @param poOptions Algorithm parameters. This should point to
264 : * GDALGridMovingAverageOptions object.
265 : * @param nPoints Number of elements in input arrays.
266 : * @param padfX Input array of X coordinates.
267 : * @param padfY Input array of Y coordinates.
268 : * @param padfZ Input array of Z values.
269 : * @param dfXPoint X coordinate of the point to compute.
270 : * @param dfYPoint Y coordinate of the point to compute.
271 : * @param pdfValue Pointer to variable where the computed grid node value
272 : * will be returned.
273 : *
274 : * @return CE_None on success or CE_Failure if something goes wrong.
275 : */
276 :
277 : CPLErr
278 0 : GDALGridMovingAverage( const void *poOptions, GUInt32 nPoints,
279 : const double *padfX, const double *padfY,
280 : const double *padfZ,
281 : double dfXPoint, double dfYPoint, double *pdfValue )
282 : {
283 : // TODO: For optimization purposes pre-computed parameters should be moved
284 : // out of this routine to the calling function.
285 :
286 : // Pre-compute search ellipse parameters
287 0 : double dfRadius1 = ((GDALGridMovingAverageOptions *)poOptions)->dfRadius1;
288 0 : double dfRadius2 = ((GDALGridMovingAverageOptions *)poOptions)->dfRadius2;
289 : double dfR12;
290 :
291 0 : dfRadius1 *= dfRadius1;
292 0 : dfRadius2 *= dfRadius2;
293 0 : dfR12 = dfRadius1 * dfRadius2;
294 :
295 : // Compute coefficients for coordinate system rotation.
296 0 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
297 : const double dfAngle =
298 0 : TO_RADIANS * ((GDALGridMovingAverageOptions *)poOptions)->dfAngle;
299 0 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
300 0 : if ( bRotated )
301 : {
302 0 : dfCoeff1 = cos(dfAngle);
303 0 : dfCoeff2 = sin(dfAngle);
304 : }
305 :
306 0 : double dfAccumulator = 0.0;
307 0 : GUInt32 i = 0, n = 0;
308 :
309 0 : while ( i < nPoints )
310 : {
311 0 : double dfRX = padfX[i] - dfXPoint;
312 0 : double dfRY = padfY[i] - dfYPoint;
313 :
314 0 : if ( bRotated )
315 : {
316 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
317 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
318 :
319 0 : dfRX = dfRXRotated;
320 0 : dfRY = dfRYRotated;
321 : }
322 :
323 : // Is this point located inside the search ellipse?
324 0 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
325 : {
326 0 : dfAccumulator += padfZ[i];
327 0 : n++;
328 : }
329 :
330 0 : i++;
331 : }
332 :
333 0 : if ( n < ((GDALGridMovingAverageOptions *)poOptions)->nMinPoints
334 : || n == 0 )
335 : {
336 : (*pdfValue) =
337 0 : ((GDALGridMovingAverageOptions *)poOptions)->dfNoDataValue;
338 : }
339 : else
340 0 : (*pdfValue) = dfAccumulator / n;
341 :
342 0 : return CE_None;
343 : }
344 :
345 : /************************************************************************/
346 : /* GDALGridNearestNeighbor() */
347 : /************************************************************************/
348 :
349 : /**
350 : * Nearest neighbor.
351 : *
352 : * The Nearest Neighbor method doesn't perform any interpolation or smoothing,
353 : * it just takes the value of nearest point found in grid node search ellipse
354 : * and returns it as a result. If there are no points found, the specified
355 : * NODATA value will be returned.
356 : *
357 : * @param poOptions Algorithm parameters. This should point to
358 : * GDALGridNearestNeighborOptions object.
359 : * @param nPoints Number of elements in input arrays.
360 : * @param padfX Input array of X coordinates.
361 : * @param padfY Input array of Y coordinates.
362 : * @param padfZ Input array of Z values.
363 : * @param dfXPoint X coordinate of the point to compute.
364 : * @param dfYPoint Y coordinate of the point to compute.
365 : * @param pdfValue Pointer to variable where the computed grid node value
366 : * will be returned.
367 : *
368 : * @return CE_None on success or CE_Failure if something goes wrong.
369 : */
370 :
371 : CPLErr
372 14641 : GDALGridNearestNeighbor( const void *poOptions, GUInt32 nPoints,
373 : const double *padfX, const double *padfY,
374 : const double *padfZ,
375 : double dfXPoint, double dfYPoint, double *pdfValue )
376 : {
377 : // TODO: For optimization purposes pre-computed parameters should be moved
378 : // out of this routine to the calling function.
379 :
380 : // Pre-compute search ellipse parameters
381 : double dfRadius1 =
382 14641 : ((GDALGridNearestNeighborOptions *)poOptions)->dfRadius1;
383 : double dfRadius2 =
384 14641 : ((GDALGridNearestNeighborOptions *)poOptions)->dfRadius2;
385 : double dfR12;
386 :
387 14641 : dfRadius1 *= dfRadius1;
388 14641 : dfRadius2 *= dfRadius2;
389 14641 : dfR12 = dfRadius1 * dfRadius2;
390 :
391 : // Compute coefficients for coordinate system rotation.
392 14641 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
393 : const double dfAngle =
394 14641 : TO_RADIANS * ((GDALGridNearestNeighborOptions *)poOptions)->dfAngle;
395 14641 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
396 14641 : if ( bRotated )
397 : {
398 0 : dfCoeff1 = cos(dfAngle);
399 0 : dfCoeff2 = sin(dfAngle);
400 : }
401 :
402 : // If the nearest point will not be found, its value remains as NODATA.
403 : double dfNearestValue =
404 14641 : ((GDALGridNearestNeighborOptions *)poOptions)->dfNoDataValue;
405 : // Nearest distance will be initialized with a largest ellipse semi-axis.
406 : // All nearest points should be located in this range.
407 14641 : double dfNearestR = MAX(dfRadius1, dfRadius2);
408 14641 : GUInt32 i = 0;
409 :
410 214388163 : while ( i < nPoints )
411 : {
412 214358881 : double dfRX = padfX[i] - dfXPoint;
413 214358881 : double dfRY = padfY[i] - dfYPoint;
414 :
415 214358881 : if ( bRotated )
416 : {
417 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
418 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
419 :
420 0 : dfRX = dfRXRotated;
421 0 : dfRY = dfRYRotated;
422 : }
423 :
424 : // Is this point located inside the search ellipse?
425 214358881 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
426 : {
427 214358881 : const double dfR2 = dfRX * dfRX + dfRY * dfRY;
428 214358881 : if ( dfNearestR == 0.0 || dfR2 < dfNearestR )
429 : {
430 8305158 : dfNearestR = dfR2;
431 8305158 : dfNearestValue = padfZ[i];
432 : }
433 : }
434 :
435 214358881 : i++;
436 : }
437 :
438 14641 : (*pdfValue) = dfNearestValue;
439 :
440 14641 : return CE_None;
441 : }
442 :
443 : /************************************************************************/
444 : /* GDALGridDataMetricMinimum() */
445 : /************************************************************************/
446 :
447 : /**
448 : * Minimum data value (data metric).
449 : *
450 : * Minimum value found in grid node search ellipse. If there are no points
451 : * found, the specified NODATA value will be returned.
452 : *
453 : * \f[
454 : * Z=\min{(Z_1,Z_2,\ldots,Z_n)}
455 : * \f]
456 : *
457 : * where
458 : * <ul>
459 : * <li> \f$Z\f$ is a resulting value at the grid node,
460 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
461 : * <li> \f$n\f$ is a total number of points in search ellipse.
462 : * </ul>
463 : *
464 : * @param poOptions Algorithm parameters. This should point to
465 : * GDALGridDataMetricsOptions object.
466 : * @param nPoints Number of elements in input arrays.
467 : * @param padfX Input array of X coordinates.
468 : * @param padfY Input array of Y coordinates.
469 : * @param padfZ Input array of Z values.
470 : * @param dfXPoint X coordinate of the point to compute.
471 : * @param dfYPoint Y coordinate of the point to compute.
472 : * @param pdfValue Pointer to variable where the computed grid node value
473 : * will be returned.
474 : *
475 : * @return CE_None on success or CE_Failure if something goes wrong.
476 : */
477 :
478 : CPLErr
479 0 : GDALGridDataMetricMinimum( const void *poOptions, GUInt32 nPoints,
480 : const double *padfX, const double *padfY,
481 : const double *padfZ,
482 : double dfXPoint, double dfYPoint, double *pdfValue )
483 : {
484 : // TODO: For optimization purposes pre-computed parameters should be moved
485 : // out of this routine to the calling function.
486 :
487 : // Pre-compute search ellipse parameters
488 : double dfRadius1 =
489 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
490 : double dfRadius2 =
491 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
492 : double dfR12;
493 :
494 0 : dfRadius1 *= dfRadius1;
495 0 : dfRadius2 *= dfRadius2;
496 0 : dfR12 = dfRadius1 * dfRadius2;
497 :
498 : // Compute coefficients for coordinate system rotation.
499 0 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
500 : const double dfAngle =
501 0 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
502 0 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
503 0 : if ( bRotated )
504 : {
505 0 : dfCoeff1 = cos(dfAngle);
506 0 : dfCoeff2 = sin(dfAngle);
507 : }
508 :
509 0 : double dfMinimumValue=0.0;
510 0 : GUInt32 i = 0, n = 0;
511 :
512 0 : while ( i < nPoints )
513 : {
514 0 : double dfRX = padfX[i] - dfXPoint;
515 0 : double dfRY = padfY[i] - dfYPoint;
516 :
517 0 : if ( bRotated )
518 : {
519 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
520 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
521 :
522 0 : dfRX = dfRXRotated;
523 0 : dfRY = dfRYRotated;
524 : }
525 :
526 : // Is this point located inside the search ellipse?
527 0 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
528 : {
529 0 : if ( n )
530 : {
531 0 : if ( dfMinimumValue > padfZ[i] )
532 0 : dfMinimumValue = padfZ[i];
533 : }
534 : else
535 0 : dfMinimumValue = padfZ[i];
536 0 : n++;
537 : }
538 :
539 0 : i++;
540 : }
541 :
542 0 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints
543 : || n == 0 )
544 : {
545 : (*pdfValue) =
546 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
547 : }
548 : else
549 0 : (*pdfValue) = dfMinimumValue;
550 :
551 0 : return CE_None;
552 : }
553 :
554 : /************************************************************************/
555 : /* GDALGridDataMetricMaximum() */
556 : /************************************************************************/
557 :
558 : /**
559 : * Maximum data value (data metric).
560 : *
561 : * Maximum value found in grid node search ellipse. If there are no points
562 : * found, the specified NODATA value will be returned.
563 : *
564 : * \f[
565 : * Z=\max{(Z_1,Z_2,\ldots,Z_n)}
566 : * \f]
567 : *
568 : * where
569 : * <ul>
570 : * <li> \f$Z\f$ is a resulting value at the grid node,
571 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
572 : * <li> \f$n\f$ is a total number of points in search ellipse.
573 : * </ul>
574 : *
575 : * @param poOptions Algorithm parameters. This should point to
576 : * GDALGridDataMetricsOptions object.
577 : * @param nPoints Number of elements in input arrays.
578 : * @param padfX Input array of X coordinates.
579 : * @param padfY Input array of Y coordinates.
580 : * @param padfZ Input array of Z values.
581 : * @param dfXPoint X coordinate of the point to compute.
582 : * @param dfYPoint Y coordinate of the point to compute.
583 : * @param pdfValue Pointer to variable where the computed grid node value
584 : * will be returned.
585 : *
586 : * @return CE_None on success or CE_Failure if something goes wrong.
587 : */
588 :
589 : CPLErr
590 0 : GDALGridDataMetricMaximum( const void *poOptions, GUInt32 nPoints,
591 : const double *padfX, const double *padfY,
592 : const double *padfZ,
593 : double dfXPoint, double dfYPoint, double *pdfValue )
594 : {
595 : // TODO: For optimization purposes pre-computed parameters should be moved
596 : // out of this routine to the calling function.
597 :
598 : // Pre-compute search ellipse parameters
599 : double dfRadius1 =
600 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
601 : double dfRadius2 =
602 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
603 : double dfR12;
604 :
605 0 : dfRadius1 *= dfRadius1;
606 0 : dfRadius2 *= dfRadius2;
607 0 : dfR12 = dfRadius1 * dfRadius2;
608 :
609 : // Compute coefficients for coordinate system rotation.
610 0 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
611 : const double dfAngle =
612 0 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
613 0 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
614 0 : if ( bRotated )
615 : {
616 0 : dfCoeff1 = cos(dfAngle);
617 0 : dfCoeff2 = sin(dfAngle);
618 : }
619 :
620 0 : double dfMaximumValue=0.0;
621 0 : GUInt32 i = 0, n = 0;
622 :
623 0 : while ( i < nPoints )
624 : {
625 0 : double dfRX = padfX[i] - dfXPoint;
626 0 : double dfRY = padfY[i] - dfYPoint;
627 :
628 0 : if ( bRotated )
629 : {
630 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
631 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
632 :
633 0 : dfRX = dfRXRotated;
634 0 : dfRY = dfRYRotated;
635 : }
636 :
637 : // Is this point located inside the search ellipse?
638 0 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
639 : {
640 0 : if ( n )
641 : {
642 0 : if ( dfMaximumValue < padfZ[i] )
643 0 : dfMaximumValue = padfZ[i];
644 : }
645 : else
646 0 : dfMaximumValue = padfZ[i];
647 0 : n++;
648 : }
649 :
650 0 : i++;
651 : }
652 :
653 0 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints
654 : || n == 0 )
655 : {
656 : (*pdfValue) =
657 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
658 : }
659 : else
660 0 : (*pdfValue) = dfMaximumValue;
661 :
662 0 : return CE_None;
663 : }
664 :
665 : /************************************************************************/
666 : /* GDALGridDataMetricRange() */
667 : /************************************************************************/
668 :
669 : /**
670 : * Data range (data metric).
671 : *
672 : * A difference between the minimum and maximum values found in grid node
673 : * search ellipse. If there are no points found, the specified NODATA
674 : * value will be returned.
675 : *
676 : * \f[
677 : * Z=\max{(Z_1,Z_2,\ldots,Z_n)}-\min{(Z_1,Z_2,\ldots,Z_n)}
678 : * \f]
679 : *
680 : * where
681 : * <ul>
682 : * <li> \f$Z\f$ is a resulting value at the grid node,
683 : * <li> \f$Z_i\f$ is a known value at point \f$i\f$,
684 : * <li> \f$n\f$ is a total number of points in search ellipse.
685 : * </ul>
686 : *
687 : * @param poOptions Algorithm parameters. This should point to
688 : * GDALGridDataMetricsOptions object.
689 : * @param nPoints Number of elements in input arrays.
690 : * @param padfX Input array of X coordinates.
691 : * @param padfY Input array of Y coordinates.
692 : * @param padfZ Input array of Z values.
693 : * @param dfXPoint X coordinate of the point to compute.
694 : * @param dfYPoint Y coordinate of the point to compute.
695 : * @param pdfValue Pointer to variable where the computed grid node value
696 : * will be returned.
697 : *
698 : * @return CE_None on success or CE_Failure if something goes wrong.
699 : */
700 :
701 : CPLErr
702 0 : GDALGridDataMetricRange( const void *poOptions, GUInt32 nPoints,
703 : const double *padfX, const double *padfY,
704 : const double *padfZ,
705 : double dfXPoint, double dfYPoint, double *pdfValue )
706 : {
707 : // TODO: For optimization purposes pre-computed parameters should be moved
708 : // out of this routine to the calling function.
709 :
710 : // Pre-compute search ellipse parameters
711 : double dfRadius1 =
712 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
713 : double dfRadius2 =
714 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
715 : double dfR12;
716 :
717 0 : dfRadius1 *= dfRadius1;
718 0 : dfRadius2 *= dfRadius2;
719 0 : dfR12 = dfRadius1 * dfRadius2;
720 :
721 : // Compute coefficients for coordinate system rotation.
722 0 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
723 : const double dfAngle =
724 0 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
725 0 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
726 0 : if ( bRotated )
727 : {
728 0 : dfCoeff1 = cos(dfAngle);
729 0 : dfCoeff2 = sin(dfAngle);
730 : }
731 :
732 0 : double dfMaximumValue=0.0, dfMinimumValue=0.0;
733 0 : GUInt32 i = 0, n = 0;
734 :
735 0 : while ( i < nPoints )
736 : {
737 0 : double dfRX = padfX[i] - dfXPoint;
738 0 : double dfRY = padfY[i] - dfYPoint;
739 :
740 0 : if ( bRotated )
741 : {
742 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
743 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
744 :
745 0 : dfRX = dfRXRotated;
746 0 : dfRY = dfRYRotated;
747 : }
748 :
749 : // Is this point located inside the search ellipse?
750 0 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
751 : {
752 0 : if ( n )
753 : {
754 0 : if ( dfMinimumValue > padfZ[i] )
755 0 : dfMinimumValue = padfZ[i];
756 0 : if ( dfMaximumValue < padfZ[i] )
757 0 : dfMaximumValue = padfZ[i];
758 : }
759 : else
760 0 : dfMinimumValue = dfMaximumValue = padfZ[i];
761 0 : n++;
762 : }
763 :
764 0 : i++;
765 : }
766 :
767 0 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints
768 : || n == 0 )
769 : {
770 : (*pdfValue) =
771 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
772 : }
773 : else
774 0 : (*pdfValue) = dfMaximumValue - dfMinimumValue;
775 :
776 0 : return CE_None;
777 : }
778 :
779 : /************************************************************************/
780 : /* GDALGridDataMetricCount() */
781 : /************************************************************************/
782 :
783 : /**
784 : * Number of data points (data metric).
785 : *
786 : * A number of data points found in grid node search ellipse.
787 : *
788 : * \f[
789 : * Z=n
790 : * \f]
791 : *
792 : * where
793 : * <ul>
794 : * <li> \f$Z\f$ is a resulting value at the grid node,
795 : * <li> \f$n\f$ is a total number of points in search ellipse.
796 : * </ul>
797 : *
798 : * @param poOptions Algorithm parameters. This should point to
799 : * GDALGridDataMetricsOptions object.
800 : * @param nPoints Number of elements in input arrays.
801 : * @param padfX Input array of X coordinates.
802 : * @param padfY Input array of Y coordinates.
803 : * @param padfZ Input array of Z values.
804 : * @param dfXPoint X coordinate of the point to compute.
805 : * @param dfYPoint Y coordinate of the point to compute.
806 : * @param pdfValue Pointer to variable where the computed grid node value
807 : * will be returned.
808 : *
809 : * @return CE_None on success or CE_Failure if something goes wrong.
810 : */
811 :
812 : CPLErr
813 0 : GDALGridDataMetricCount( const void *poOptions, GUInt32 nPoints,
814 : const double *padfX, const double *padfY,
815 : const double *padfZ,
816 : double dfXPoint, double dfYPoint, double *pdfValue )
817 : {
818 : // TODO: For optimization purposes pre-computed parameters should be moved
819 : // out of this routine to the calling function.
820 :
821 : // Pre-compute search ellipse parameters
822 : double dfRadius1 =
823 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
824 : double dfRadius2 =
825 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
826 : double dfR12;
827 :
828 0 : dfRadius1 *= dfRadius1;
829 0 : dfRadius2 *= dfRadius2;
830 0 : dfR12 = dfRadius1 * dfRadius2;
831 :
832 : // Compute coefficients for coordinate system rotation.
833 0 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
834 : const double dfAngle =
835 0 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
836 0 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
837 0 : if ( bRotated )
838 : {
839 0 : dfCoeff1 = cos(dfAngle);
840 0 : dfCoeff2 = sin(dfAngle);
841 : }
842 :
843 0 : GUInt32 i = 0, n = 0;
844 :
845 0 : while ( i < nPoints )
846 : {
847 0 : double dfRX = padfX[i] - dfXPoint;
848 0 : double dfRY = padfY[i] - dfYPoint;
849 :
850 0 : if ( bRotated )
851 : {
852 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
853 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
854 :
855 0 : dfRX = dfRXRotated;
856 0 : dfRY = dfRYRotated;
857 : }
858 :
859 : // Is this point located inside the search ellipse?
860 0 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
861 0 : n++;
862 :
863 0 : i++;
864 : }
865 :
866 0 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints )
867 : {
868 : (*pdfValue) =
869 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
870 : }
871 : else
872 0 : (*pdfValue) = (double)n;
873 :
874 0 : return CE_None;
875 : }
876 :
877 : /************************************************************************/
878 : /* GDALGridDataMetricAverageDistance() */
879 : /************************************************************************/
880 :
881 : /**
882 : * Average distance (data metric).
883 : *
884 : * An average distance between the grid node (center of the search ellipse)
885 : * and all of the data points found in grid node search ellipse. If there are
886 : * no points found, the specified NODATA value will be returned.
887 : *
888 : * \f[
889 : * Z=\frac{\sum_{i = 1}^n r_i}{n}
890 : * \f]
891 : *
892 : * where
893 : * <ul>
894 : * <li> \f$Z\f$ is a resulting value at the grid node,
895 : * <li> \f$r_i\f$ is an Euclidean distance from the grid node
896 : * to point \f$i\f$,
897 : * <li> \f$n\f$ is a total number of points in search ellipse.
898 : * </ul>
899 : *
900 : * @param poOptions Algorithm parameters. This should point to
901 : * GDALGridDataMetricsOptions object.
902 : * @param nPoints Number of elements in input arrays.
903 : * @param padfX Input array of X coordinates.
904 : * @param padfY Input array of Y coordinates.
905 : * @param padfZ Input array of Z values.
906 : * @param dfXPoint X coordinate of the point to compute.
907 : * @param dfYPoint Y coordinate of the point to compute.
908 : * @param pdfValue Pointer to variable where the computed grid node value
909 : * will be returned.
910 : *
911 : * @return CE_None on success or CE_Failure if something goes wrong.
912 : */
913 :
914 : CPLErr
915 0 : GDALGridDataMetricAverageDistance( const void *poOptions, GUInt32 nPoints,
916 : const double *padfX, const double *padfY,
917 : const double *padfZ,
918 : double dfXPoint, double dfYPoint,
919 : double *pdfValue )
920 : {
921 : // TODO: For optimization purposes pre-computed parameters should be moved
922 : // out of this routine to the calling function.
923 :
924 : // Pre-compute search ellipse parameters
925 : double dfRadius1 =
926 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
927 : double dfRadius2 =
928 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
929 : double dfR12;
930 :
931 0 : dfRadius1 *= dfRadius1;
932 0 : dfRadius2 *= dfRadius2;
933 0 : dfR12 = dfRadius1 * dfRadius2;
934 :
935 : // Compute coefficients for coordinate system rotation.
936 0 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
937 : const double dfAngle =
938 0 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
939 0 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
940 0 : if ( bRotated )
941 : {
942 0 : dfCoeff1 = cos(dfAngle);
943 0 : dfCoeff2 = sin(dfAngle);
944 : }
945 :
946 0 : double dfAccumulator = 0.0;
947 0 : GUInt32 i = 0, n = 0;
948 :
949 0 : while ( i < nPoints )
950 : {
951 0 : double dfRX = padfX[i] - dfXPoint;
952 0 : double dfRY = padfY[i] - dfYPoint;
953 :
954 0 : if ( bRotated )
955 : {
956 0 : double dfRXRotated = dfRX * dfCoeff1 + dfRY * dfCoeff2;
957 0 : double dfRYRotated = dfRY * dfCoeff1 - dfRX * dfCoeff2;
958 :
959 0 : dfRX = dfRXRotated;
960 0 : dfRY = dfRYRotated;
961 : }
962 :
963 : // Is this point located inside the search ellipse?
964 0 : if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
965 : {
966 0 : dfAccumulator += sqrt( dfRX * dfRX + dfRY * dfRY );
967 0 : n++;
968 : }
969 :
970 0 : i++;
971 : }
972 :
973 0 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints
974 : || n == 0 )
975 : {
976 : (*pdfValue) =
977 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
978 : }
979 : else
980 0 : (*pdfValue) = dfAccumulator / n;
981 :
982 0 : return CE_None;
983 : }
984 :
985 : /************************************************************************/
986 : /* GDALGridDataMetricAverageDistance() */
987 : /************************************************************************/
988 :
989 : /**
990 : * Average distance between points (data metric).
991 : *
992 : * An average distance between the data points found in grid node search
993 : * ellipse. The distance between each pair of points within ellipse is
994 : * calculated and average of all distances is set as a grid node value. If
995 : * there are no points found, the specified NODATA value will be returned.
996 :
997 : *
998 : * \f[
999 : * Z=\frac{\sum_{i = 1}^{n-1}\sum_{j=i+1}^{n} r_{ij}}{\left(n-1\right)\,n-\frac{n+{\left(n-1\right)}^{2}-1}{2}}
1000 : * \f]
1001 : *
1002 : * where
1003 : * <ul>
1004 : * <li> \f$Z\f$ is a resulting value at the grid node,
1005 : * <li> \f$r_{ij}\f$ is an Euclidean distance between points
1006 : * \f$i\f$ and \f$j\f$,
1007 : * <li> \f$n\f$ is a total number of points in search ellipse.
1008 : * </ul>
1009 : *
1010 : * @param poOptions Algorithm parameters. This should point to
1011 : * GDALGridDataMetricsOptions object.
1012 : * @param nPoints Number of elements in input arrays.
1013 : * @param padfX Input array of X coordinates.
1014 : * @param padfY Input array of Y coordinates.
1015 : * @param padfZ Input array of Z values.
1016 : * @param dfXPoint X coordinate of the point to compute.
1017 : * @param dfYPoint Y coordinate of the point to compute.
1018 : * @param pdfValue Pointer to variable where the computed grid node value
1019 : * will be returned.
1020 : *
1021 : * @return CE_None on success or CE_Failure if something goes wrong.
1022 : */
1023 :
1024 : CPLErr
1025 0 : GDALGridDataMetricAverageDistancePts( const void *poOptions, GUInt32 nPoints,
1026 : const double *padfX, const double *padfY,
1027 : const double *padfZ,
1028 : double dfXPoint, double dfYPoint,
1029 : double *pdfValue )
1030 : {
1031 : // TODO: For optimization purposes pre-computed parameters should be moved
1032 : // out of this routine to the calling function.
1033 :
1034 : // Pre-compute search ellipse parameters
1035 : double dfRadius1 =
1036 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius1;
1037 : double dfRadius2 =
1038 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfRadius2;
1039 : double dfR12;
1040 :
1041 0 : dfRadius1 *= dfRadius1;
1042 0 : dfRadius2 *= dfRadius2;
1043 0 : dfR12 = dfRadius1 * dfRadius2;
1044 :
1045 : // Compute coefficients for coordinate system rotation.
1046 0 : double dfCoeff1 = 0.0, dfCoeff2 = 0.0;
1047 : const double dfAngle =
1048 0 : TO_RADIANS * ((GDALGridDataMetricsOptions *)poOptions)->dfAngle;
1049 0 : const bool bRotated = ( dfAngle == 0.0 ) ? false : true;
1050 0 : if ( bRotated )
1051 : {
1052 0 : dfCoeff1 = cos(dfAngle);
1053 0 : dfCoeff2 = sin(dfAngle);
1054 : }
1055 :
1056 0 : double dfAccumulator = 0.0;
1057 0 : GUInt32 i = 0, n = 0;
1058 :
1059 : // Search for the first point within the search ellipse
1060 0 : while ( i < nPoints - 1 )
1061 : {
1062 0 : double dfRX1 = padfX[i] - dfXPoint;
1063 0 : double dfRY1 = padfY[i] - dfYPoint;
1064 :
1065 0 : if ( bRotated )
1066 : {
1067 0 : double dfRXRotated = dfRX1 * dfCoeff1 + dfRY1 * dfCoeff2;
1068 0 : double dfRYRotated = dfRY1 * dfCoeff1 - dfRX1 * dfCoeff2;
1069 :
1070 0 : dfRX1 = dfRXRotated;
1071 0 : dfRY1 = dfRYRotated;
1072 : }
1073 :
1074 : // Is this point located inside the search ellipse?
1075 0 : if ( dfRadius2 * dfRX1 * dfRX1 + dfRadius1 * dfRY1 * dfRY1 <= dfR12 )
1076 : {
1077 : GUInt32 j;
1078 :
1079 : // Search all the remaining points within the ellipse and compute
1080 : // distances between them and the first point
1081 0 : for ( j = i + 1; j < nPoints; j++ )
1082 : {
1083 0 : double dfRX2 = padfX[j] - dfXPoint;
1084 0 : double dfRY2 = padfY[j] - dfYPoint;
1085 :
1086 0 : if ( bRotated )
1087 : {
1088 0 : double dfRXRotated = dfRX2 * dfCoeff1 + dfRY2 * dfCoeff2;
1089 0 : double dfRYRotated = dfRY2 * dfCoeff1 - dfRX2 * dfCoeff2;
1090 :
1091 0 : dfRX2 = dfRXRotated;
1092 0 : dfRY2 = dfRYRotated;
1093 : }
1094 :
1095 0 : if ( dfRadius2 * dfRX2 * dfRX2 + dfRadius1 * dfRY2 * dfRY2 <= dfR12 )
1096 : {
1097 0 : const double dfRX = padfX[j] - padfX[i];
1098 0 : const double dfRY = padfY[j] - padfY[i];
1099 :
1100 0 : dfAccumulator += sqrt( dfRX * dfRX + dfRY * dfRY );
1101 0 : n++;
1102 : }
1103 : }
1104 : }
1105 :
1106 0 : i++;
1107 : }
1108 :
1109 0 : if ( n < ((GDALGridDataMetricsOptions *)poOptions)->nMinPoints
1110 : || n == 0 )
1111 : {
1112 : (*pdfValue) =
1113 0 : ((GDALGridDataMetricsOptions *)poOptions)->dfNoDataValue;
1114 : }
1115 : else
1116 0 : (*pdfValue) = dfAccumulator / n;
1117 :
1118 0 : return CE_None;
1119 : }
1120 :
1121 : /************************************************************************/
1122 : /* GDALGridCreate() */
1123 : /************************************************************************/
1124 :
1125 : /**
1126 : * Create regular grid from the scattered data.
1127 : *
1128 : * This fucntion takes the arrays of X and Y coordinates and corresponding Z
1129 : * values as input and computes regular grid (or call it a raster) from these
1130 : * scattered data. You should supply geometry and extent of the output grid
1131 : * and allocate array sufficient to hold such a grid.
1132 : *
1133 : * @param eAlgorithm Gridding method.
1134 : * @param poOptions Options to control choosen gridding method.
1135 : * @param nPoints Number of elements in input arrays.
1136 : * @param padfX Input array of X coordinates.
1137 : * @param padfY Input array of Y coordinates.
1138 : * @param padfZ Input array of Z values.
1139 : * @param dfXMin Lowest X border of output grid.
1140 : * @param dfXMax Highest X border of output grid.
1141 : * @param dfYMin Lowest Y border of output grid.
1142 : * @param dfYMax Highest Y border of output grid.
1143 : * @param nXSize Number of columns in output grid.
1144 : * @param nYSize Number of rows in output grid.
1145 : * @param eType Data type of output array.
1146 : * @param pData Pointer to array where the computed grid will be stored.
1147 : * @param pfnProgress a GDALProgressFunc() compatible callback function for
1148 : * reporting progress or NULL.
1149 : * @param pProgressArg argument to be passed to pfnProgress. May be NULL.
1150 : *
1151 : * @return CE_None on success or CE_Failure if something goes wrong.
1152 : */
1153 :
1154 : CPLErr
1155 1 : GDALGridCreate( GDALGridAlgorithm eAlgorithm, const void *poOptions,
1156 : GUInt32 nPoints,
1157 : const double *padfX, const double *padfY, const double *padfZ,
1158 : double dfXMin, double dfXMax, double dfYMin, double dfYMax,
1159 : GUInt32 nXSize, GUInt32 nYSize, GDALDataType eType, void *pData,
1160 : GDALProgressFunc pfnProgress, void *pProgressArg )
1161 : {
1162 : CPLAssert( poOptions );
1163 : CPLAssert( padfX );
1164 : CPLAssert( padfY );
1165 : CPLAssert( padfZ );
1166 : CPLAssert( pData );
1167 :
1168 1 : if ( pfnProgress == NULL )
1169 0 : pfnProgress = GDALDummyProgress;
1170 :
1171 1 : if ( nXSize == 0 || nYSize == 0 )
1172 : {
1173 : CPLError( CE_Failure, CPLE_IllegalArg,
1174 0 : "Output raster dimesions should have non-zero size." );
1175 0 : return CE_Failure;
1176 : }
1177 :
1178 : GDALGridFunction pfnGDALGridMethod;
1179 :
1180 1 : switch ( eAlgorithm )
1181 : {
1182 : case GGA_InverseDistanceToAPower:
1183 0 : if ( ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->
1184 : dfRadius1 == 0.0 &&
1185 : ((GDALGridInverseDistanceToAPowerOptions *)poOptions)->
1186 : dfRadius2 == 0.0 )
1187 0 : pfnGDALGridMethod = GDALGridInverseDistanceToAPowerNoSearch;
1188 : else
1189 0 : pfnGDALGridMethod = GDALGridInverseDistanceToAPower;
1190 0 : break;
1191 :
1192 : case GGA_MovingAverage:
1193 0 : pfnGDALGridMethod = GDALGridMovingAverage;
1194 0 : break;
1195 :
1196 : case GGA_NearestNeighbor:
1197 1 : pfnGDALGridMethod = GDALGridNearestNeighbor;
1198 1 : break;
1199 :
1200 : case GGA_MetricMinimum:
1201 0 : pfnGDALGridMethod = GDALGridDataMetricMinimum;
1202 0 : break;
1203 :
1204 : case GGA_MetricMaximum:
1205 0 : pfnGDALGridMethod = GDALGridDataMetricMaximum;
1206 0 : break;
1207 :
1208 : case GGA_MetricRange:
1209 0 : pfnGDALGridMethod = GDALGridDataMetricRange;
1210 0 : break;
1211 :
1212 : case GGA_MetricCount:
1213 0 : pfnGDALGridMethod = GDALGridDataMetricCount;
1214 0 : break;
1215 :
1216 : case GGA_MetricAverageDistance:
1217 0 : pfnGDALGridMethod = GDALGridDataMetricAverageDistance;
1218 0 : break;
1219 :
1220 : case GGA_MetricAverageDistancePts:
1221 0 : pfnGDALGridMethod = GDALGridDataMetricAverageDistancePts;
1222 0 : break;
1223 :
1224 : default:
1225 : CPLError( CE_Failure, CPLE_IllegalArg,
1226 0 : "GDAL does not support gridding method %d", eAlgorithm );
1227 0 : return CE_Failure;
1228 : }
1229 :
1230 : GUInt32 nXPoint, nYPoint;
1231 1 : const double dfDeltaX = ( dfXMax - dfXMin ) / nXSize;
1232 1 : const double dfDeltaY = ( dfYMax - dfYMin ) / nYSize;
1233 :
1234 : /* -------------------------------------------------------------------- */
1235 : /* Allocate a buffer of scanline size, fill it with gridded values */
1236 : /* and use GDALCopyWords() to copy values into output data array with */
1237 : /* appropriate data type conversion. */
1238 : /* -------------------------------------------------------------------- */
1239 1 : double *padfValues = (double *)VSIMalloc( sizeof(double) * nXSize );
1240 1 : GByte *pabyData = (GByte *)pData;
1241 1 : int nDataTypeSize = GDALGetDataTypeSize(eType) / 8;
1242 1 : int nLineSpace = nXSize * nDataTypeSize;
1243 :
1244 122 : for ( nYPoint = 0; nYPoint < nYSize; nYPoint++ )
1245 : {
1246 121 : const double dfYPoint = dfYMin + ( nYPoint + 0.5 ) * dfDeltaY;
1247 :
1248 14762 : for ( nXPoint = 0; nXPoint < nXSize; nXPoint++ )
1249 : {
1250 14641 : const double dfXPoint = dfXMin + ( nXPoint + 0.5 ) * dfDeltaX;
1251 :
1252 14641 : if ( (*pfnGDALGridMethod)( poOptions, nPoints, padfX, padfY, padfZ,
1253 : dfXPoint, dfYPoint,
1254 : padfValues + nXPoint ) != CE_None )
1255 : {
1256 : CPLError( CE_Failure, CPLE_AppDefined,
1257 : "Gridding failed at X position %lu, Y position %lu",
1258 : (long unsigned int)nXPoint,
1259 0 : (long unsigned int)nYPoint );
1260 0 : return CE_Failure;
1261 : }
1262 : }
1263 :
1264 : GDALCopyWords( padfValues, GDT_Float64, sizeof(double),
1265 : pabyData, eType, nDataTypeSize,
1266 121 : nXSize );
1267 121 : pabyData += nLineSpace;
1268 :
1269 121 : if( !pfnProgress( (double)(nYPoint + 1) / nYSize, NULL, pProgressArg ) )
1270 : {
1271 0 : CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
1272 0 : return CE_Failure;
1273 : }
1274 : }
1275 :
1276 1 : VSIFree( padfValues );
1277 :
1278 1 : return CE_None;
1279 : }
1280 :
|
