[gdal-dev] Kriging Interpolation

[Ian Reese]

Hi All,

Thank you so much everyone for your replies.  I should be clear, speed is important to our processes, but most important in my work is the visual outcome.  Kriging by far produces the best results for this. No amount of tweaking IDW or any other method in GDAL produces the same results. Kriging works best with scattered data and produces the nicest results across nodata locations.  I have plenty of methods for working around speed issues and have worked with any number of interpolations methods in GRASS, SAGA and Python. What I don't have is a Kriging interpolation in gdal_grid. I love using gdal for its ease, speed, and efficiency. I am merely asking if it is possible incorporate kriging interpolation into gdal_grid, is there interest in the wider community to do so, and what would it take to help make that happen.

It is my fault that I was not clear about the examples I provided.  I offered the most stripped down examples, but was not interested in the speed of the returns. I am more interested in the comparisons between the visual outcomes of the outputs.

Thanks again for all the responses.

Cheers,

Ian

-----Original Message-----
From: Even Rouault [mailto: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.