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