[gdal-dev] Kriging Interpolation

[Even Rouault]

Ian,

Looking at your script, for the IDW, I see you don't use the max_points, nor 
the radius1 and radius2 to limit the neighbourhoud into which the points are 
searched.

For the linear search, I've profiled the run and found that an enormous amout 
of time was spent processing the points outside of the triangulation. 
Investigating, I found there was a missed optimization that I've added per
https://github.com/OSGeo/gdal/commit/134b9c4e437e347aac48c459a18a933ed526be5b

With that optimization, the runtime with 8 threads is 16 sec, vs 48 sec before 
(on a non-optimized debug build)
And if you don't want the nearest neighbour search for points outside of the 
triangulation, you can also add :radius=0.0, and this will speed up things 
again. Down to 2 sec.

Regarding Kriging, I guess you specified some neighbourood value in ArcGIS ? 
Otherwise as mentionned by Thomas, Knudsen, the computation time and memory 
requirements would be enormous (between O(n^2) and O(n^3) where n is the 
number of points).
I gave a try to https://github.com/bsmurphy/PyKrige and used Ordinary Kriging. 
For the above reasons, it blows up with a MemoryError on the full layer. I 
also found that the extent of the output grid had also a strong influence on 
the memory requirements. Which is less expected, perhaps a side effect of that 
particular implementation.
Anyway, if arbitrarily restricting to the first 1000 points of your layer 
(which is roughly a strip in the middle of the dataset) and to their extent 
(producing a 1321x45 raster with 1m resolution), it gives a result (2 minutes 
single-threaded, 1.8 GB RAM), which looks consistent with the one from ArcGIS. 
This is not exactly the same due to a non-moving neighbouroud having been used 
due to my simplified methodology, but for the 'inner' of the computed raster, 
this is quite similar.
So for Kriging, a moving window approach limiting the number of input points 
would be absolutely needed for "large" point datasets, and the processing time 
would still be significantly slower than linear interpolation.

Even


> Hi Gdal,
> 
> I'd like to know if it is possible to make Kriging Interpolation part of
> gdal_grid.  We work with a good deal of scattered point elevation data and
> neither IDW or Linear interpolations(smooth or not smooth)  produce the
> desired visual results.  Nor are they fast enough for the density of data
> we commonly work with. Linear interpolation is a reliable and fast method,
> however the triangulated results are difficult to work with visually.
> 
> For some visual comparisons and reference, here is a link to a folder of
> sample points(shp) for a test area, two examples of kriged data for that
> test region, and a BASH script to reproduce IDW and Linear interpolations
> of the same region.
> 
> https://drive.google.com/open?id=1HKbK4wMmVN6JZSvhxBQDub1WXBLE710H
> 
> Here is an example of the Kriging method used by ArcGIS:
> http://desktop.arcgis.com/en/arcmap/latest/extensions/geostatistical-analyst
> /kriging-in-geostatistical-analyst.htm
> 
> Few graphs explaining best interpolation methods to use with given data. 
> Gridded vs. scattered.
> https://www.neonscience.org/spatial-interpolation-basics
> 
> If it is possible, we are curious to know what type of funding would be
> needed.
> 
> Cheers,
> 
> Ian
> 
> ________________________________
> 
> This message contains information, which may be in confidence and may be
> subject to legal privilege. If you are not the intended recipient, you must
> not peruse, use, disseminate, distribute or copy this message. If you have
> received this message in error, please notify us immediately (Phone 0800
> 665 463 or info at linz.govt.nz) and destroy the original message. LINZ
> accepts no responsibility for changes to this email, or for any
> attachments, after its transmission from LINZ. Thank You.


-- 
Spatialys - Geospatial professional services
http://www.spatialys.com