6

I am looking for some fast and robust open source interpolation library/script/tool.

The goal is to re-interpolate relatively big (1M+ points) unstructured (xyz) grids into high resolution cartesian grids (tif files or other format I can process with GDAL) to produce nice smooth contours.

So far the only usable way I found is to use triangulation in saga_cmd but the more points I have the longer it takes and from some higher number it is not finished even after hours while if I do the same in other GIS software through GUI it takes only seconds to maximum few minutes.

I wanted to use GDAL but though nearest is fast, it doesn't get me detail I need and other methods took ages even for lower numbers of points... I am not sure if/how could I use GRASS. One of my requirements is that the solution must be portable (so far I used portable python, gdal and saga).

Based on answers so far:

  1. Example of data enter image description here

  2. gdal_fillnodata is incredibly fast but leaving spikes no matter how much smoothed, same with SAGA close gaps, where the speed is highly dependent on tension threshold, the bigger the faster enter image description here

  3. all (2 and nearest) gdal interpolations are not really suited for this task - nearest just takes the nearest value (not really interpolation), moving average is one step better but creates big steps not suitable for contours, the IDW is very hard to set to some reasonable results and anyway it takes too long enter image description here

  4. so far the best results I got are still with triangulation. Bus as I mentioned it takes long time with SAGA for middle size grids and for big it fails enter image description here

There are still few things I haven't tried yet or not tested properly:

Both I am not sure about the speed and will be glad for any experience or better sugestions. Thanks.

Improve this question
8
  • Did you look into grass.osgeo.org/grass70/manuals/r.in.xyz.html ?
    – markusN
    Aug 18, 2014 at 23:23
  • Triangulation is difficult... been there, done that! It gave me more respect for the Esri developers. I did a web search for delaunay triangulation and found some C/C++ code that works. The problem with triangulation is that the processing time increases exponentially with the number of points; for the Esri TIN commands a rule-of-thumb is to keep the number of points below 1Meg, you may need to tile your points with overlap and then merge post rasterization.
    – Michael Stimson
    Aug 18, 2014 at 23:50
  • @markusN Is it possible to have GRASS portable on Windows? Should I start new question on that topic? All I found are installation packages. Also it seems to me that function does the same as gdal_rasterize - convert points to raster. As I mentioned I have unstructured (irregular) and/or rectilinear grids meaning I need to interpolate values for cells between the points of known values. For that there is mentioned v.surf.rst but after further reading that is too slow to process millions of points.
    – Miro
    Aug 19, 2014 at 0:26
  • @MichaelMiles-Stimson I can break points to tiles but there is this tricky issue with triangulation that it interpolates only between known points meaning breaking into tiles will cause holes with no data values on the edges of tiles... Any links for the compiled C/C++ library/tool I could use in command line?
    – Miro
    Aug 19, 2014 at 0:35
  • 1
    @MichaelMiles-Stimson it is plugin for QGIS called Qgis2threejs - plugins.qgis.org/plugins/Qgis2threejs - it can export whatever you have in the map (raster/vector) to 3D visualization based on nice option of settings you can use - all as html/javascript/images witch you can view in any modern web browser. It is great tool, I really like it. There are only few things I miss there but might be added in the future :)
    – Miro
    Apr 24, 2015 at 23:47

4 Answers 4

Reset to default
3

The legendary triangulation is available http://www.cs.cmu.edu/~quake/triangle.html

from the page:

I timed the Delaunay triangulation of 1,000,000 vertices uniformly randomly distributed in a square. The output contained 1,999,955 triangles. (I used the -I switch to suppress the rewriting of input vertices to another file, since the vertices written would be identical to the ones read. Hence, only the triangles were written to disk.)

Delaunay triangulation time: 20.761 sec
File output time: 28.794 sec
Total time (above plus file input): 56.115 sec

This has been translated to a python library according to http://dzhelil.info/triangle/index.html# which includes documentation, download https://pypi.python.org/pypi/triangle/2013.04.05 - note, these links may become dead as the version changes. I found them with google in under 10 minutes.

Improve this answer
8
  • Hi, thanks, I found that one before but was frightened by it's age (newest version from 2005) and moreover with no clue how to operate with it to produce cartesian grid or contours from that...
    – Miro
    Aug 19, 2014 at 0:56
  • I was scared by the original triangle code too, I read the source and couldn't make any sense from it; neither could I understand the required format for files. The good thing is the python library - read the docs, it's easy as! According to the statement it should triangulate your points in less than one minute (don't know the specs of what it was benchmarked on).
    – Michael Stimson
    Aug 19, 2014 at 1:00
  • I mean let's say I will be able to put together all the python bindings, make it portable and run it on my points to create the triangulation mesh. So I will get triangulation mesh in some unspecified text format. How should I use that to get to my goal which is interpolated cartesian grid or contours?
    – Miro
    Aug 19, 2014 at 1:20
  • I'm not sure what the output of the function is, it looks like either a list or a dictionary object. It looks like you would just iterate over the result and create points using OGR or similar. The only way to tell is install it, do a triangulation on an arbitrary set and see what object is returned.
    – Michael Stimson
    Aug 19, 2014 at 1:33
  • I am sure that the first part of triangulation interpolation is creating the triangulation mesh and this tool seems to be super robust and fast at that - create 2D triangles from points. But it ends there. I omit it doesn't create triangles based on XYZ but only on XY which is fine with me. But there is still completely missing the part which will interpolate Z values based on these triangles for new points into Cartesian grid.
    – Miro
    Aug 19, 2014 at 2:56
3

In saga, make sure you read your data as a pointcloud "Import Point Cloud from Text File" and process it as a pointcloud: Point Cloud to Grid. This should be much faster than using shape. Last step may be close gaps to fill any non data values.

I would not recommend using triangulation if you just want a grid. the algorithm takes long and is not better than others.

Improve this answer
4
  • 1
    Thank you. In that case I would probably stick with GDAL which has gdal_rasterize and gdal_fillnodata. Actually I am not sure why I left that idea, should try it again because it could work very well in this case and it is definitely fastest and very robust option.
    – Miro
    Aug 19, 2014 at 8:14
  • 1
    I should add that triangulation would actually be a bad idea if you want to have smooth contour lines (they will have sharp edges all over if you have little data). I think that indeed GDAL is probably the best option (I was not aware of gdal_fillnodata). As a tip, you could consider working with in memory rasters in gdal so you don't have to write/read to disk all the time.
    – johanvdw
    Aug 19, 2014 at 8:51
  • Thanks for tips. I am going to try both methods - compare SAGA vs GDAL and write result here but it will take me few days to get to it now, thanks again.
    – Miro
    Aug 20, 2014 at 23:09
  • gdal_fillnodata is incredibly fast compared to any interpolation but it doesn't work well for too sparse points which is unfortunately my case. SAGA has more methods and I am still through it... Also I found this post sugesting using Scipy: gis.stackexchange.com/questions/17432/…
    – Miro
    Sep 9, 2014 at 1:22
2

I've never used the tool myself, but points2grid was designed to be a lightweight tool to interpolate large point data sets to a structured grid. It is found in a few places:

Finally, to create contours, try using gdal_contour with the interpolated raster.

Improve this answer
4
  • Installation of GEON-points2grid-Utility-Setup-v1_3 on my Windows 7 machine fails with error so I can't test it. Anyway that is accesible only through GUI which is not my point. The more recent version on github looks promising but no portable version for Windows, not even installer. Well if there is no easier way I will challenge myself with it.
    – Miro
    Aug 19, 2014 at 3:09
  • @Miro Do you know how to install programs using OSGeo4W? If so, look there for a precompiled version for Windows. If not, I can explain further how to install this to an OSGeo4W shell. It looks like it is only available for the 64-bit version.
    – Mike T
    Aug 19, 2014 at 4:27
  • Thank you. I have installed point2grid through advanced install of OSGeo4W and have now pts2grd.lib in C:\OSGeo4W64\lib. Now the lame question, how do I use it?
    – Miro
    Aug 19, 2014 at 5:21
  • From an OSGeo4W shell, first see points2grid --help. For a real-world example, download example.las to the working directory and run points2grid -i example.las -o output. The result is a mess of files with the prefix output.
    – Mike T
    Aug 19, 2014 at 5:35
1

I have ended up using Python with GDAL etc as it comes as easily installable and after that portable package when installing QGIS OSGeo4W. Matplotlib - matplotlib.mlab.griddata (interp='linear') for superfast interpolation into regular grid. On that can by applied matplotlib.pyplot.contour to very fast contouring. Package osgeo.ogr and shapely to save to GIS formats. The core of my code was hugely inspired by python QGIS plugin called "Contour plugin" - https://github.com/ccrook/QGIS-Contour-Plugin (Plugin built by Lionel Roubeyrie. Contributions from Chris Crook)

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.