← Back to branch summary

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

« back to all changes in this revision

Viewing changes to apps/gdal_grid.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:
apps/gdal_grid.cpp
1
1
/* ****************************************************************************
2
 
 * $Id: gdal_grid.cpp 18482 2010-01-08 23:22:58Z mloskot $
 
2
 * $Id: gdal_grid.cpp 22783 2011-07-23 19:28:16Z rouault $
3
3
 *
4
4
 * Project:  GDAL Utilities
5
5
 * Purpose:  GDAL scattered data gridding (interpolation) tool
37
37
#include "ogr_spatialref.h"
38
38
#include "ogr_api.h"
39
39
#include "ogrsf_frmts.h"
40
 
 
41
 
CPL_CVSID("$Id: gdal_grid.cpp 18482 2010-01-08 23:22:58Z mloskot $");
42
 
 
43
 
static const char szAlgNameInvDist[] = "invdist";
44
 
static const char szAlgNameAverage[] = "average";
45
 
static const char szAlgNameNearest[] = "nearest";
46
 
static const char szAlgNameMinimum[] = "minimum";
47
 
static const char szAlgNameMaximum[] = "maximum";
48
 
static const char szAlgNameRange[] = "range";
49
 
static const char szAlgNameCount[] = "count";
50
 
static const char szAlgNameAverageDistance[] = "average_distance";
51
 
static const char szAlgNameAverageDistancePts[] = "average_distance_pts";
 
40
#include "gdalgrid.h"
 
41
#include "commonutils.h"
 
42
 
 
43
CPL_CVSID("$Id: gdal_grid.cpp 22783 2011-07-23 19:28:16Z rouault $");
52
44
 
53
45
/************************************************************************/
54
46
/*                               Usage()                                */
213
205
}
214
206
 
215
207
/************************************************************************/
216
 
/*                      ParseAlgorithmAndOptions()                      */
217
 
/*                                                                      */
218
 
/*      Translates mnemonic gridding algorithm names into               */
219
 
/*      GDALGridAlgorithm code, parse control parameters and assign     */
220
 
/*      defaults.                                                       */
221
 
/************************************************************************/
222
 
 
223
 
static CPLErr ParseAlgorithmAndOptions( const char *pszAlgoritm,
224
 
                                        GDALGridAlgorithm *peAlgorithm,
225
 
                                        void **ppOptions )
226
 
{
227
 
    char **papszParms = CSLTokenizeString2( pszAlgoritm, ":", FALSE );
228
 
 
229
 
    if ( CSLCount(papszParms) < 1 )
230
 
        return CE_Failure;
231
 
 
232
 
    if ( EQUAL(papszParms[0], szAlgNameInvDist) )
233
 
        *peAlgorithm = GGA_InverseDistanceToAPower;
234
 
    else if ( EQUAL(papszParms[0], szAlgNameAverage) )
235
 
        *peAlgorithm = GGA_MovingAverage;
236
 
    else if ( EQUAL(papszParms[0], szAlgNameNearest) )
237
 
        *peAlgorithm = GGA_NearestNeighbor;
238
 
    else if ( EQUAL(papszParms[0], szAlgNameMinimum) )
239
 
        *peAlgorithm = GGA_MetricMinimum;
240
 
    else if ( EQUAL(papszParms[0], szAlgNameMaximum) )
241
 
        *peAlgorithm = GGA_MetricMaximum;
242
 
    else if ( EQUAL(papszParms[0], szAlgNameRange) )
243
 
        *peAlgorithm = GGA_MetricRange;
244
 
    else if ( EQUAL(papszParms[0], szAlgNameCount) )
245
 
        *peAlgorithm = GGA_MetricCount;
246
 
    else if ( EQUAL(papszParms[0], szAlgNameAverageDistance) )
247
 
        *peAlgorithm = GGA_MetricAverageDistance;
248
 
    else if ( EQUAL(papszParms[0], szAlgNameAverageDistancePts) )
249
 
        *peAlgorithm = GGA_MetricAverageDistancePts;
250
 
    else
251
 
    {
252
 
        fprintf( stderr, "Unsupported gridding method \"%s\".\n",
253
 
                 papszParms[0] );
254
 
        CSLDestroy( papszParms );
255
 
        return CE_Failure;
256
 
    }
257
 
 
258
 
/* -------------------------------------------------------------------- */
259
 
/*      Parse algorithm parameters and assign defaults.                 */
260
 
/* -------------------------------------------------------------------- */
261
 
    const char  *pszValue;
262
 
 
263
 
    switch ( *peAlgorithm )
264
 
    {
265
 
        case GGA_InverseDistanceToAPower:
266
 
        default:
267
 
            *ppOptions =
268
 
                CPLMalloc( sizeof(GDALGridInverseDistanceToAPowerOptions) );
269
 
 
270
 
            pszValue = CSLFetchNameValue( papszParms, "power" );
271
 
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
272
 
                dfPower = (pszValue) ? atof(pszValue) : 2.0;
273
 
 
274
 
            pszValue = CSLFetchNameValue( papszParms, "smoothing" );
275
 
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
276
 
                dfSmoothing = (pszValue) ? atof(pszValue) : 0.0;
277
 
 
278
 
            pszValue = CSLFetchNameValue( papszParms, "radius1" );
279
 
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
280
 
                dfRadius1 = (pszValue) ? atof(pszValue) : 0.0;
281
 
 
282
 
            pszValue = CSLFetchNameValue( papszParms, "radius2" );
283
 
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
284
 
                dfRadius2 = (pszValue) ? atof(pszValue) : 0.0;
285
 
 
286
 
            pszValue = CSLFetchNameValue( papszParms, "angle" );
287
 
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
288
 
                dfAngle = (pszValue) ? atof(pszValue) : 0.0;
289
 
 
290
 
            pszValue = CSLFetchNameValue( papszParms, "max_points" );
291
 
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
292
 
                nMaxPoints = (pszValue) ? atol(pszValue) : 0;
293
 
 
294
 
            pszValue = CSLFetchNameValue( papszParms, "min_points" );
295
 
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
296
 
                nMinPoints = (pszValue) ? atol(pszValue) : 0;
297
 
 
298
 
            pszValue = CSLFetchNameValue( papszParms, "nodata" );
299
 
            ((GDALGridInverseDistanceToAPowerOptions *)*ppOptions)->
300
 
                dfNoDataValue = (pszValue) ? atof(pszValue) : 0.0;
301
 
            break;
302
 
 
303
 
        case GGA_MovingAverage:
304
 
            *ppOptions =
305
 
                CPLMalloc( sizeof(GDALGridMovingAverageOptions) );
306
 
 
307
 
            pszValue = CSLFetchNameValue( papszParms, "radius1" );
308
 
            ((GDALGridMovingAverageOptions *)*ppOptions)->
309
 
                dfRadius1 = (pszValue) ? atof(pszValue) : 0.0;
310
 
 
311
 
            pszValue = CSLFetchNameValue( papszParms, "radius2" );
312
 
            ((GDALGridMovingAverageOptions *)*ppOptions)->
313
 
                dfRadius2 = (pszValue) ? atof(pszValue) : 0.0;
314
 
 
315
 
            pszValue = CSLFetchNameValue( papszParms, "angle" );
316
 
            ((GDALGridMovingAverageOptions *)*ppOptions)->
317
 
                dfAngle = (pszValue) ? atof(pszValue) : 0.0;
318
 
 
319
 
            pszValue = CSLFetchNameValue( papszParms, "min_points" );
320
 
            ((GDALGridMovingAverageOptions *)*ppOptions)->
321
 
                nMinPoints = (pszValue) ? atol(pszValue) : 0;
322
 
 
323
 
            pszValue = CSLFetchNameValue( papszParms, "nodata" );
324
 
            ((GDALGridMovingAverageOptions *)*ppOptions)->
325
 
                dfNoDataValue = (pszValue) ? atof(pszValue) : 0.0;
326
 
            break;
327
 
 
328
 
        case GGA_NearestNeighbor:
329
 
            *ppOptions =
330
 
                CPLMalloc( sizeof(GDALGridNearestNeighborOptions) );
331
 
 
332
 
            pszValue = CSLFetchNameValue( papszParms, "radius1" );
333
 
            ((GDALGridNearestNeighborOptions *)*ppOptions)->
334
 
                dfRadius1 = (pszValue) ? atof(pszValue) : 0.0;
335
 
 
336
 
            pszValue = CSLFetchNameValue( papszParms, "radius2" );
337
 
            ((GDALGridNearestNeighborOptions *)*ppOptions)->
338
 
                dfRadius2 = (pszValue) ? atof(pszValue) : 0.0;
339
 
 
340
 
            pszValue = CSLFetchNameValue( papszParms, "angle" );
341
 
            ((GDALGridNearestNeighborOptions *)*ppOptions)->
342
 
                dfAngle = (pszValue) ? atof(pszValue) : 0.0;
343
 
 
344
 
            pszValue = CSLFetchNameValue( papszParms, "nodata" );
345
 
            ((GDALGridNearestNeighborOptions *)*ppOptions)->
346
 
                dfNoDataValue = (pszValue) ? atof(pszValue) : 0.0;
347
 
            break;
348
 
 
349
 
        case GGA_MetricMinimum:
350
 
        case GGA_MetricMaximum:
351
 
        case GGA_MetricRange:
352
 
        case GGA_MetricCount:
353
 
        case GGA_MetricAverageDistance:
354
 
        case GGA_MetricAverageDistancePts:
355
 
            *ppOptions =
356
 
                CPLMalloc( sizeof(GDALGridDataMetricsOptions) );
357
 
 
358
 
            pszValue = CSLFetchNameValue( papszParms, "radius1" );
359
 
            ((GDALGridDataMetricsOptions *)*ppOptions)->
360
 
                dfRadius1 = (pszValue) ? atof(pszValue) : 0.0;
361
 
 
362
 
            pszValue = CSLFetchNameValue( papszParms, "radius2" );
363
 
            ((GDALGridDataMetricsOptions *)*ppOptions)->
364
 
                dfRadius2 = (pszValue) ? atof(pszValue) : 0.0;
365
 
 
366
 
            pszValue = CSLFetchNameValue( papszParms, "angle" );
367
 
            ((GDALGridDataMetricsOptions *)*ppOptions)->
368
 
                dfAngle = (pszValue) ? atof(pszValue) : 0.0;
369
 
 
370
 
            pszValue = CSLFetchNameValue( papszParms, "min_points" );
371
 
            ((GDALGridDataMetricsOptions *)*ppOptions)->
372
 
                nMinPoints = (pszValue) ? atol(pszValue) : 0;
373
 
 
374
 
            pszValue = CSLFetchNameValue( papszParms, "nodata" );
375
 
            ((GDALGridDataMetricsOptions *)*ppOptions)->
376
 
                dfNoDataValue = (pszValue) ? atof(pszValue) : 0.0;
377
 
            break;
378
 
 
379
 
   }
380
 
 
381
 
    CSLDestroy( papszParms );
382
 
    return CE_None;
383
 
}
384
 
 
385
 
/************************************************************************/
386
208
/*                          ProcessGeometry()                           */
387
209
/*                                                                      */
388
210
/*  Extract point coordinates from the geometry reference and set the   */
394
216
                             int iBurnField, double dfBurnValue,
395
217
                             std::vector<double> &adfX,
396
218
                             std::vector<double> &adfY,
397
 
                             std::vector<double> &adfZ)
 
219
                             std::vector<double> &adfZ )
398
220
 
399
221
{
400
222
    if ( poClipSrc && !poGeom->Within(poClipSrc) )
488
310
        OGRFeature::DestroyFeature( poFeat );
489
311
    }
490
312
 
491
 
    if (adfX.size() == 0)
 
313
    if ( adfX.size() == 0 )
492
314
    {
493
315
        printf( "No point geometry found on layer %s, skipping.\n",
494
316
                OGR_FD_GetName( OGR_L_GetLayerDefn( hSrcLayer ) ) );
498
320
/* -------------------------------------------------------------------- */
499
321
/*      Compute grid geometry.                                          */
500
322
/* -------------------------------------------------------------------- */
501
 
 
502
 
    if ( !bIsXExtentSet )
503
 
    {
504
 
        dfXMin = *std::min_element(adfX.begin(), adfX.end());
505
 
        dfXMax = *std::max_element(adfX.begin(), adfX.end());
506
 
        bIsXExtentSet = TRUE;
507
 
    }
508
 
 
509
 
    if ( !bIsYExtentSet )
510
 
    {
511
 
        dfYMin = *std::min_element(adfY.begin(), adfY.end());
512
 
        dfYMax = *std::max_element(adfY.begin(), adfY.end());
513
 
        bIsYExtentSet = TRUE;
 
323
    if ( !bIsXExtentSet || !bIsYExtentSet )
 
324
    {
 
325
        OGREnvelope sEnvelope;
 
326
        OGR_L_GetExtent( hSrcLayer, &sEnvelope, TRUE );
 
327
 
 
328
        if ( !bIsXExtentSet )
 
329
        {
 
330
            dfXMin = sEnvelope.MinX;
 
331
            dfXMax = sEnvelope.MaxX;
 
332
            bIsXExtentSet = TRUE;
 
333
        }
 
334
 
 
335
        if ( !bIsYExtentSet )
 
336
        {
 
337
            dfYMin = sEnvelope.MinY;
 
338
            dfYMax = sEnvelope.MaxY;
 
339
            bIsYExtentSet = TRUE;
 
340
        }
514
341
    }
515
342
 
516
343
/* -------------------------------------------------------------------- */
686
513
{
687
514
    GDALDriverH     hDriver;
688
515
    const char      *pszSource=NULL, *pszDest=NULL, *pszFormat = "GTiff";
 
516
    int             bFormatExplicitelySet = FALSE;
689
517
    char            **papszLayers = NULL;
690
518
    const char      *pszBurnAttribute = NULL;
691
519
    const char      *pszWHERE = NULL, *pszSQL = NULL;
733
561
        else if( EQUAL(argv[i],"-of") && i < argc-1 )
734
562
        {
735
563
            pszFormat = argv[++i];
 
564
            bFormatExplicitelySet = TRUE;
736
565
        }
737
566
 
738
567
        else if( EQUAL(argv[i],"-q") || EQUAL(argv[i],"-quiet") )
961
790
            Usage();
962
791
        }
963
792
    }
964
 
    else if ( bClipSrc && poClipSrc == NULL )
965
 
    {
966
 
        if ( poSpatialFilter )
967
 
            poClipSrc = poSpatialFilter->clone();
968
 
        if ( poClipSrc == NULL )
969
 
        {
970
 
            fprintf( stderr,
971
 
                     "FAILURE: -clipsrc must be used with -spat option or \n"
972
 
                     "a bounding box, WKT string or datasource must be "
973
 
                     "specified\n\n" );
974
 
            Usage();
 
793
    else if ( bClipSrc && poClipSrc == NULL && !poSpatialFilter )
 
794
    {
 
795
        fprintf( stderr,
 
796
                 "FAILURE: -clipsrc must be used with -spat option or \n"
 
797
                 "a bounding box, WKT string or datasource must be "
 
798
                 "specified\n\n" );
 
799
        Usage();
 
800
    }
 
801
 
 
802
    if ( poSpatialFilter )
 
803
    {
 
804
        if ( poClipSrc )
 
805
        {
 
806
            OGRGeometry *poTemp = poSpatialFilter->Intersection( poClipSrc );
 
807
 
 
808
            if ( poTemp )
 
809
            {
 
810
                OGRGeometryFactory::destroyGeometry( poSpatialFilter );
 
811
                poSpatialFilter = poTemp;
 
812
            }
 
813
 
 
814
            OGRGeometryFactory::destroyGeometry( poClipSrc );
 
815
            poClipSrc = NULL;
 
816
        }
 
817
    }
 
818
    else
 
819
    {
 
820
        if ( poClipSrc )
 
821
        {
 
822
            poSpatialFilter = poClipSrc;
 
823
            poClipSrc = NULL;
975
824
        }
976
825
    }
977
826
 
1034
883
    if ( nYSize == 0 )
1035
884
        nYSize = 256;
1036
885
 
 
886
    if (!bQuiet && !bFormatExplicitelySet)
 
887
        CheckExtensionConsistency(pszDest, pszFormat);
 
888
 
1037
889
    hDstDS = GDALCreate( hDriver, pszDest, nXSize, nYSize, nBands,
1038
890
                         eOutputType, papszCreateOptions );
1039
891
    if ( hDstDS == NULL )
1062
914
        if( hLayer != NULL )
1063
915
        {
1064
916
            // Custom layer will be rasterized in the first band.
1065
 
            ProcessLayer( hLayer, hDstDS, poClipSrc, nXSize, nYSize, 1,
 
917
            ProcessLayer( hLayer, hDstDS, poSpatialFilter, nXSize, nYSize, 1,
1066
918
                          bIsXExtentSet, bIsYExtentSet,
1067
919
                          dfXMin, dfXMax, dfYMin, dfYMax, pszBurnAttribute,
1068
920
                          eOutputType, eAlgorithm, pOptions,
1100
952
                OSRExportToWkt( hSRS, &pszOutputSRS );
1101
953
        }
1102
954
 
1103
 
        ProcessLayer( hLayer, hDstDS, poClipSrc, nXSize, nYSize,
 
955
        ProcessLayer( hLayer, hDstDS, poSpatialFilter, nXSize, nYSize,
1104
956
                      i + 1 + nBands - nLayerCount,
1105
957
                      bIsXExtentSet, bIsYExtentSet,
1106
958
                      dfXMin, dfXMax, dfYMin, dfYMax, pszBurnAttribute,
1135
987
    OGR_DS_Destroy( hSrcDS );
1136
988
    GDALClose( hDstDS );
1137
989
    OGRGeometryFactory::destroyGeometry( poSpatialFilter );
1138
 
    OGRGeometryFactory::destroyGeometry( poClipSrc );
1139
990
 
1140
991
    CPLFree( pOptions );
1141
992
    CSLDestroy( papszCreateOptions );