New gdal_grid interpolation algorithm

[mckelvym]

Hello,

I would like to contribute the attached patch for a new algorithm to be included in gdal_grid. This algorithm option, "invdistnn", is a variation on the existing inverse distance weighting algorithm with the following features:

• Use a quadtree to search for points only in the neighborhood (within radius) of each grid cell. This dramatically improves on the time to produce an interpolated grid since it is not examining the entire point space like "invdist"
• Retains up to "max_points" nearest cells instead of arbitrarily the first "max_points + 1" cells

This implementation will produce interpolated grids that match the output of ESRI's ArcPy arcpy.gp.Idw_sa method.

The affected files are:

• alg/gdal_alg.h
• alg/gdalgrid.cpp
• alg/gdalgrid.h
• apps/gdal_grid.cpp

[mckelvym]

Implementation of "invdistnn" for gdal_grid

[Even Rouault]

I didn't actually try it, but it looks good to me. It would also be necessary to document the new alg in apps/gdal_utilities.dox. And having a new test case in autotest/utilities/test_gdal_grid.py would be great also.

[mckelvym]

Patch to document invdistnn in apps/gdal_utilities.dox

[mckelvym]

Patch to test_gdal_grid to add test 12 for invdistnn

[mckelvym]

utilities/ref_data/grid_invdistnn*tif reference testing data

[mckelvym]

utilities/ref_data/grid_invdistnn*tif reference testing data

[mckelvym]

utilities/ref_data/grid_invdistnn*tif reference testing data

[mckelvym]

Please let me know if there is anything else. Thanks.

[Even Rouault]

I'm middly enthousiastic in seeing computed fields in GDALGridInverseDistanceToAPowerNearestNeighborOptions such as dfPowerDiv2PreComp, dfRadiusPower2PreComp and dfRadiusPower4PreComp. Cannot they just be computed in GDALGridInverseDistanceToAPowerNearestNeighbor() ? If their precomputation really improve performance, couldn't those value be passed with the hExtraParamsIn structure ?

[mckelvym]

Yes, these could be computed in GDALGridInverseDistanceToAPowerNearestNeighbor(), but since they are grid cell invariant, it is smarter to compute them in advance. I think hExtraParamsIn may be a better place to introduce this, as you suggest.

The motivation for computing them in advance is a dataset I am working with that has 29 million grid cells -- these extra multiplications * 29+ million is wasted cpu time.

Do you suggest I make these modifications to move dfPowerDiv2PreComp, dfRadiusPower2PreComp, and dfRadiusPower4PreComp from the structure GDALGridInverseDistanceToAPowerNearestNeighborOptions to this structure (GDALGridExtraParameters) instead? If so, I can submit a new patch to use instead of the first one.

Thanks

[Even Rouault]

Do you suggest I make these modifications to move dfPowerDiv2PreComp, dfRadiusPower2PreComp, and dfRadiusPower4PreComp from the structure GDALGridInverseDistanceToAPowerNearestNeighborOptions to this structure (GDALGridExtraParameters) instead? If so, I can submit a new patch to use instead of the first one.

Yes, that would be better

[mckelvym]

Implementation of "invdistnn" for gdal_grid. This patch replaces gdal_grid_invdistnn_20150715.patch

[mckelvym]

rouault, please see attachment:gdal_grid_invdistnn_20150720.patch​ which contains the changes described above. Use this one instead of gdal_grid_invdistnn_20150715.patch

Note: This patch was generated against [29540].

[Even Rouault]

trunk r29541, r29542 "gdal_grid: add invdistnn algorithm, variation on the existing inverse distance weighting algorithm with quadtree to search for points only in the neighborhood (patch by mckelvym, #6038)"