[gdal-dev] Kriging Interpolation

[Ian Reese]

Hi Evan,

I probably wasn't clear in my response to you.  I did try optimizing IDW  using max points and various radii.  The results were mixed, but overall by using radius and max points, the process is slowed considerably.  If it is helpful, I can give you some time differences when I get back to testing.  Usually I can work around speed issues by reducing the area I am interpolating, but there does reach a point where the area is too small for interpolating and artifacts near the edges become an issue for visualization.

As far as kriging goes, for testing and using the exact same dataset, I used the default settings in Arc.  So:

Kriging Method: Ordinary
SemiVariogram: Spherical
Search Radius: Variable
Number of Points: 12

For test examples, I tired to keep things as simple as possible since the more complex solutions were providing similar visual result as the basic examples.

I really appreciate you looking into the processes I provided.  I plan to upgrade my Gdal in few days too.  You suggested in a previous thread that doing so should solve some of the speed issues.  I let you know how it turned out.

Thanks again,

Ian
________________________________________
From: Even Rouault [even.rouault at spatialys.com]
Sent: Tuesday, 26 June 2018 12:21 a.m.
To: gdal-dev at lists.osgeo.org
Cc: Ian Reese; Jeremy Palmer
Subject: Re: [gdal-dev] Kriging Interpolation

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

________________________________

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.