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