[gdal-dev] Kriging Interpolation

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

[gdal-dev] Kriging Interpolation

Ian Reese
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 [hidden email]) 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.

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Kriging Interpolation

Chris Marsh
For what is worth, as an end-user of GDAL, I would be very interested in having such capacity in GDAL.
Cheers
Chris


On 24 June 2018 at 21:19, Ian Reese <[hidden email]> wrote:
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 [hidden email]) 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.

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Kriging Interpolation

Thomas Knudsen
In reply to this post by Ian Reese
Hi Ian,

If you find IDW on N points slow, you will find Kriging on the same number of points unbearably slow, since for Kriging you will have to solve a system of N linear equations (essentially inverting an NxN matrix) to compute the Kriging weights. To set up the matrix you will even need to compute the same distance elements you need for IDW. So Kriging will never become faster than IDW.

Hence, the only realistic way to use Kriging on large data sets is to base each estimate on local subsets - typically by working in a TIN model, just like in the linear interpolation case. And if your data set is dense, the difference between Kriging and linear interpolation will be miniscule. Kriging will, however, provide a variance estimate for the prediction, which may be realistic assuming you provide a realistic covariance function for the phenomenon under investigation.

So to make any difference wrt the linear interpolation, you will need to include a larger number of data, e.g. by including neighbouring triangles of the TIN, or using another data structure to subset your observations.

Some years ago I wrote an informal note about this subject. The title is “Data Assimilation for Updates of Digital Terrain Models” and it is available over at http://sdfe.dk/media/2916626/kms_technical_report_14.pdf

/Thomas




2018-06-25 5:19 GMT+02:00 Ian Reese <[hidden email]>:
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 [hidden email]) 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.

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Kriging Interpolation

Even Rouault-2
In reply to this post by Ian Reese
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 [hidden email]) 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
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: 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:[hidden email]]
Sent: Tuesday, 26 June 2018 12:21 a.m.
To: [hidden email]
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 [hidden email]) 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 [hidden email]) 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.
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Kriging Interpolation

Joaquim Luis


I have plenty of methods for working around speed issues and have
|>worked with any number of interpolations methods in GRASS, SAGA and
|>Python.


You have also the minimum curvature in GMT. Not lightning fast but quicker than krigging and normally better looking, specially on the extrapolation zones.

Joaquim

 
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Kriging Interpolation

Ian Reese
Cheer Joakim!  That I am not familiar with and will have a look
________________________________________
From: Joaquim Manuel Freire Luís [[hidden email]]
Sent: Tuesday, 26 June 2018 10:58 a.m.
To: Ian Reese; Even Rouault; [hidden email]
Cc: Jeremy Palmer
Subject: RE: [gdal-dev] Kriging Interpolation

I have plenty of methods for working around speed issues and have
|>worked with any number of interpolations methods in GRASS, SAGA and
|>Python.


You have also the minimum curvature in GMT. Not lightning fast but quicker than krigging and normally better looking, specially on the extrapolation zones.

Joaquim



________________________________

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 [hidden email]) 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.
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Kriging Interpolation

Ian Reese
In reply to this post by Even Rouault-2
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 [[hidden email]]
Sent: Tuesday, 26 June 2018 12:21 a.m.
To: [hidden email]
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 [hidden email]) 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 [hidden email]) 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.
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev