← Back to branch summary

~ubuntu-branches/debian/sid/gdal/sid

« back to all changes in this revision

Viewing changes to alg/gdalgrid.cpp

  • Committer: Package Import Robot
  • Author(s): Francesco Paolo Lovergine
  • Date: 2012-05-07 15:04:42 UTC
  • mfrom: (5.5.16 experimental)
  • Revision ID: package-import@ubuntu.com-20120507150442-2eks97loeh6rq005
Tags: 1.9.0-1
* Ready for sid, starting transition.
* All symfiles updated to latest builds.
* Added dh_numpy call in debian/rules to depend on numpy ABI.
* Policy bumped to 3.9.3, no changes required.

Show diffs side-by-side

added added

removed removed

Lines of Context:
alg/gdalgrid.cpp
1
1
/******************************************************************************
2
 
 * $Id: gdalgrid.cpp 18041 2009-11-17 14:14:43Z dron $
 
2
 * $Id: gdalgrid.cpp 22889 2011-08-07 13:07:14Z rouault $
3
3
 *
4
4
 * Project:  GDAL Gridding API.
5
5
 * Purpose:  Implementation of GDAL scattered data gridder.
28
28
 ****************************************************************************/
29
29
 
30
30
#include "cpl_vsi.h"
 
31
#include "cpl_string.h"
31
32
#include "gdalgrid.h"
 
33
#include <float.h>
 
34
#include <limits.h>
32
35
 
33
 
CPL_CVSID("$Id: gdalgrid.cpp 18041 2009-11-17 14:14:43Z dron $");
 
36
CPL_CVSID("$Id: gdalgrid.cpp 22889 2011-08-07 13:07:14Z rouault $");
34
37
 
35
38
#define TO_RADIANS (3.14159265358979323846 / 180.0)
36
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
 
37
48
/************************************************************************/
38
49
/*                   GDALGridInverseDistanceToAPower()                  */
39
50
/************************************************************************/
145
156
        // Is this point located inside the search ellipse?
146
157
        if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
147
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.
148
161
            if ( CPLIsEqual(dfR2, 0.0) )
149
162
            {
150
163
                (*pdfValue) = padfZ[i];
210
223
        const double dfR2 =
211
224
            dfRX * dfRX + dfRY * dfRY + dfSmoothing * dfSmoothing;
212
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.
213
228
        if ( CPLIsEqual(dfR2, 0.0) )
214
229
        {
215
230
            (*pdfValue) = padfZ[i];
402
417
    // If the nearest point will not be found, its value remains as NODATA.
403
418
    double      dfNearestValue =
404
419
        ((GDALGridNearestNeighborOptions *)poOptions)->dfNoDataValue;
405
 
    // Nearest distance will be initialized with a largest ellipse semi-axis.
406
 
    // All nearest points should be located in this range.
407
 
    double      dfNearestR = MAX(dfRadius1, dfRadius2);
 
420
    // Nearest distance will be initialized with the distance to the first
 
421
    // point in array.
 
422
    double      dfNearestR = DBL_MAX;
408
423
    GUInt32 i = 0;
409
424
 
410
425
    while ( i < nPoints )
425
440
        if ( dfRadius2 * dfRX * dfRX + dfRadius1 * dfRY * dfRY <= dfR12 )
426
441
        {
427
442
            const double    dfR2 = dfRX * dfRX + dfRY * dfRY;
428
 
            if ( dfNearestR == 0.0 || dfR2 < dfNearestR )
 
443
            if ( dfR2 <= dfNearestR )
429
444
            {
430
445
                dfNearestR = dfR2;
431
446
                dfNearestValue = padfZ[i];
1125
1140
/**
1126
1141
 * Create regular grid from the scattered data.
1127
1142
 *
1128
 
 * This fucntion takes the arrays of X and Y coordinates and corresponding Z
 
1143
 * This function takes the arrays of X and Y coordinates and corresponding Z
1129
1144
 * values as input and computes regular grid (or call it a raster) from these
1130
1145
 * scattered data. You should supply geometry and extent of the output grid
1131
1146
 * and allocate array sufficient to hold such a grid.
1278
1293
    return CE_None;
1279
1294
}
1280
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
CPLErr ParseAlgorithmAndOptions( const char *pszAlgoritm,
 
1305
                                 GDALGridAlgorithm *peAlgorithm,
 
1306
                                 void **ppOptions )
 
1307
{
 
1308
    CPLAssert( pszAlgoritm );
 
1309
    CPLAssert( peAlgorithm );
 
1310
    CPLAssert( ppOptions );
 
1311
 
 
1312
    *ppOptions = NULL;
 
1313
 
 
1314
    char **papszParms = CSLTokenizeString2( pszAlgoritm, ":", FALSE );
 
1315
 
 
1316
    if ( CSLCount(papszParms) < 1 )
 
1317
        return CE_Failure;
 
1318
 
 
1319
    if ( EQUAL(papszParms[0], szAlgNameInvDist) )
 
1320
        *peAlgorithm = GGA_InverseDistanceToAPower;
 
1321
    else if ( EQUAL(papszParms[0], szAlgNameAverage) )
 
1322
        *peAlgorithm = GGA_MovingAverage;
 
1323
    else if ( EQUAL(papszParms[0], szAlgNameNearest) )
 
1324
        *peAlgorithm = GGA_NearestNeighbor;
 
1325
    else if ( EQUAL(papszParms[0], szAlgNameMinimum) )
 
1326
        *peAlgorithm = GGA_MetricMinimum;
 
1327
    else if ( EQUAL(papszParms[0], szAlgNameMaximum) )
 
1328
        *peAlgorithm = GGA_MetricMaximum;
 
1329
    else if ( EQUAL(papszParms[0], szAlgNameRange) )
 
1330
        *peAlgorithm = GGA_MetricRange;
 
1331
    else if ( EQUAL(papszParms[0], szAlgNameCount) )
 
1332
        *peAlgorithm = GGA_MetricCount;
 
1333
    else if ( EQUAL(papszParms[0], szAlgNameAverageDistance) )
 
1334
        *peAlgorithm = GGA_MetricAverageDistance;
 
1335
    else if ( EQUAL(papszParms[0], szAlgNameAverageDistancePts) )
 
1336
        *peAlgorithm = GGA_MetricAverageDistancePts;
 
1337
    else
 
1338
    {
 
1339
        fprintf( stderr, "Unsupported gridding method \"%s\".\n",
 
1340
                 papszParms[0] );
 
1341
        CSLDestroy( papszParms );
 
1342
        return CE_Failure;
 
1343
    }
 
1344
 
 
1345
/* -------------------------------------------------------------------- */
 
1346
/*      Parse algorithm parameters and assign defaults.                 */
 
1347
/* -------------------------------------------------------------------- */
 
1348
    const char  *pszValue;
 
1349
 
 
1350
    switch ( *peAlgorithm )
 
1351
    {
 
1352
        case GGA_InverseDistanceToAPower:
 
1353
        default:
 
1354
            *ppOptions =
 
1355
                CPLMalloc( sizeof(GDALGridInverseDistanceToAPowerOptions) );
 
1356
 
 
1357
            pszValue = CSLFetchNameValue( papszParms, "power" );
 
1358
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
 
1359
                dfPower = (pszValue) ? CPLAtofM(pszValue) : 2.0;
 
1360
 
 
1361
            pszValue = CSLFetchNameValue( papszParms, "smoothing" );
 
1362
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
 
1363
                dfSmoothing = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1364
 
 
1365
            pszValue = CSLFetchNameValue( papszParms, "radius1" );
 
1366
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
 
1367
                dfRadius1 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1368
 
 
1369
            pszValue = CSLFetchNameValue( papszParms, "radius2" );
 
1370
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
 
1371
                dfRadius2 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1372
 
 
1373
            pszValue = CSLFetchNameValue( papszParms, "angle" );
 
1374
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
 
1375
                dfAngle = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1376
 
 
1377
            pszValue = CSLFetchNameValue( papszParms, "max_points" );
 
1378
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
 
1379
                nMaxPoints = (GUInt32) ((pszValue) ? CPLAtofM(pszValue) : 0);
 
1380
 
 
1381
            pszValue = CSLFetchNameValue( papszParms, "min_points" );
 
1382
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
 
1383
                nMinPoints = (GUInt32) ((pszValue) ? CPLAtofM(pszValue) : 0);
 
1384
 
 
1385
            pszValue = CSLFetchNameValue( papszParms, "nodata" );
 
1386
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
 
1387
                dfNoDataValue = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1388
            break;
 
1389
 
 
1390
        case GGA_MovingAverage:
 
1391
            *ppOptions =
 
1392
                CPLMalloc( sizeof(GDALGridMovingAverageOptions) );
 
1393
 
 
1394
            pszValue = CSLFetchNameValue( papszParms, "radius1" );
 
1395
            ((GDALGridMovingAverageOptions *)*ppOptions)->
 
1396
                dfRadius1 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1397
 
 
1398
            pszValue = CSLFetchNameValue( papszParms, "radius2" );
 
1399
            ((GDALGridMovingAverageOptions *)*ppOptions)->
 
1400
                dfRadius2 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1401
 
 
1402
            pszValue = CSLFetchNameValue( papszParms, "angle" );
 
1403
            ((GDALGridMovingAverageOptions *)*ppOptions)->
 
1404
                dfAngle = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1405
 
 
1406
            pszValue = CSLFetchNameValue( papszParms, "min_points" );
 
1407
            ((GDALGridMovingAverageOptions *)*ppOptions)->
 
1408
                nMinPoints = (GUInt32) ((pszValue) ? CPLAtofM(pszValue) : 0);
 
1409
 
 
1410
            pszValue = CSLFetchNameValue( papszParms, "nodata" );
 
1411
            ((GDALGridMovingAverageOptions *)*ppOptions)->
 
1412
                dfNoDataValue = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1413
            break;
 
1414
 
 
1415
        case GGA_NearestNeighbor:
 
1416
            *ppOptions =
 
1417
                CPLMalloc( sizeof(GDALGridNearestNeighborOptions) );
 
1418
 
 
1419
            pszValue = CSLFetchNameValue( papszParms, "radius1" );
 
1420
            ((GDALGridNearestNeighborOptions *)*ppOptions)->
 
1421
                dfRadius1 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1422
 
 
1423
            pszValue = CSLFetchNameValue( papszParms, "radius2" );
 
1424
            ((GDALGridNearestNeighborOptions *)*ppOptions)->
 
1425
                dfRadius2 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1426
 
 
1427
            pszValue = CSLFetchNameValue( papszParms, "angle" );
 
1428
            ((GDALGridNearestNeighborOptions *)*ppOptions)->
 
1429
                dfAngle = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1430
 
 
1431
            pszValue = CSLFetchNameValue( papszParms, "nodata" );
 
1432
            ((GDALGridNearestNeighborOptions *)*ppOptions)->
 
1433
                dfNoDataValue = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1434
            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
                CPLMalloc( sizeof(GDALGridDataMetricsOptions) );
 
1444
 
 
1445
            pszValue = CSLFetchNameValue( papszParms, "radius1" );
 
1446
            ((GDALGridDataMetricsOptions *)*ppOptions)->
 
1447
                dfRadius1 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1448
 
 
1449
            pszValue = CSLFetchNameValue( papszParms, "radius2" );
 
1450
            ((GDALGridDataMetricsOptions *)*ppOptions)->
 
1451
                dfRadius2 = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1452
 
 
1453
            pszValue = CSLFetchNameValue( papszParms, "angle" );
 
1454
            ((GDALGridDataMetricsOptions *)*ppOptions)->
 
1455
                dfAngle = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1456
 
 
1457
            pszValue = CSLFetchNameValue( papszParms, "min_points" );
 
1458
            ((GDALGridDataMetricsOptions *)*ppOptions)->
 
1459
                nMinPoints = (pszValue) ? atol(pszValue) : 0;
 
1460
 
 
1461
            pszValue = CSLFetchNameValue( papszParms, "nodata" );
 
1462
            ((GDALGridDataMetricsOptions *)*ppOptions)->
 
1463
                dfNoDataValue = (pszValue) ? CPLAtofM(pszValue) : 0.0;
 
1464
            break;
 
1465
 
 
1466
   }
 
1467
 
 
1468
    CSLDestroy( papszParms );
 
1469
    return CE_None;
 
1470
}
 
1471