Manifold currently offers two algorithms to create gridded datasets. Below is a list of the algorithms most commonly used to create gridded datasets, and which I would like to see included within a future version of Manifold.

• Nearest neighbour
• Moving average
• Natural neighbour
• Minimum curvature
• Radial basis function
• Modified Shepard's method
• Triangulation with linear interpolation
• Point Kriging
• Block kriging

Also, two more recently developed algorithms/methodologies

• Franklin Approximation (c.f. http://www.terrainmap.com/rm32.html and http://www.terrainmap.com/downloads/finalreport.pdf
• Univariate curve interpolation (c.f. C¹-continuous terrain resconstrurction from sparse contours http://www.multires.caltech.edu/pubs/TerrRecon.pdf )

Does anyone know of any other methods, that should be added to make a definitive list?

Notes:

1. Manifold currently offers a choice of kriging and median Polish kriging. However, the documentation with respect to what is termed "kriging" is not very clear as to what form is being used. So, for completness, I have listed above the two commonest forms of "plain" kriging, one or neither of which may correspond to what is currently available within Manifold.
2. The last two algorithms listed are state-of-the art and information suggests they would seem to be ideal for artifact-free DEM production from digitised or vectorised scans of paper maps, if they can be successfully implemented.

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

kernel density - with and without regression

A relevant reference:

Bowman & Azzalini (1997) 'Applied Smoothing Techniques for Data Analysis: The Kernel Approach with S-Plus Illustrations' Oxford University Press.

http://mirror.aarnet.edu.au/pub/CRAN/src/contrib/Descriptions/sm.html

Thanks Mike , the more the better...

Of related interest may be some of the ideas presented in:

http://www.ecse.rpi.edu/Homepages/wrf/research/ica99/

Particularly the ones concerning lossless compression of gridded data and the automated creation of auxillery contours from a vector drawing of contours without gridding.

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

I jotted down a few ideas from the list you gave, and went and talked to a local guru - I was thinking to organize those methods within a true math perspective and give some structure to the ideas. Sure enough, he concurred that it's just such a big topic *even he* can't do that - it's just so open-ended. (And he would know, believe me 8)

A couple of points, and thanks for those links - that contouring reference is fascinating:

I assume you are referring to methods to create surfaces that represent topography?

(of course, it's a much bigger topic if you extend it from those domains, which is what I'm thinking about - for example, even within maths the methods really tend to be application-specific - think of 3-D reconstructions from 2-D scans in medicine for one massive field of work, and then multidimensional analysis where the surface is just a visual aid for slices across some informative dimension)

I've been following your discussion of surfaces as "cells" vs. a matrix of points closely, I do tend to think of them always as matrices that hold data that fills a "grid" that we specify - but the difference is very important analytically, conceptually, programmatically, and visually. etc.

I'm used to dealing with "surfaces" as matrices within a code context (mostly R), that and the grid structure is imposed upon that. My linear algebra is very limited, but I've identified it as an absolute key area I need to develop (for example). I sometimes think I'd like Manifold to allow me to manipulate surfaces with commands that forget the grid structure (i.e. as a matrix that I can transpose, solve, ...) - but R does that so welll there's no point in asking for it - and it already can do much of it it with transforms, just not how I like - I should point out. Manifold is about spatial data primarily, I suppose.

You first 3 methods are really ways of creating gridded point data - *before* you've even begun the task of constructing a surface from which those data are taken. No?

Can you elaborate on each method some more? For example, "nearest neighbour" has many definitions in many different contexts, and many of the others are names that are application specific.

The big divisions I see are regular vs. irregular point data - everything else stems from that first issue.

Hi Mike ,

Thanks for the comments.

First off I am not expert, "I just fly 'em - I dont build 'em". So If I suggest something that you think is erroneous, speak out.

I am would like to see Manifold have interpolators that are designed for both general gridding and ones that are are more suited, but not exclusively so, for the special case of the gridding of topographic elevation data.

The way I would classify interpolation or gridding algoriths is on the basis of whether they are exact or smoothing/averaging interpolators, or can be either. See below (mainly from Surfer manual). All these methods are gridding algorithms in that they produce arrays of evenly spaced z-values in the xy plane (grids)

Exact interpolators

Inverse distance to a power, if no smoothing factor is specified

Krigging, point or block, when used without specifying an error nugget effect

Nearest neighbour under all circumstances

Radial basis function when an R² value is not specified

Modified Shepard's method when a smoothing factor is not specified

Triangulation with linear interpolation

Natural neighbour

and from information online:

Franklin approxiamtion

Univariate curve interpolation

Smoothing or averaging interpolators

Inverse distance to a power when a smoothing factor is specified

Krigging, point or block, when an error nugget effect is specified

Radial basis function when an R2 values is specified

Modified Shepard's method when a smoothing factor is specified

Moving average

I think we can safely limit interpolators to those that are used in geoinfomatics in its broadest sense. So we can exclude the more obscure ones. The above list is based on what is currently offered by Golden Software's Surfer 8 (which BTW is a good model for Manifold to aspire to) with a few additions.

The additions are Franklin Approximation( a la John Chile's work) and Univariate curve interpolation. I have included these because, although somewhat new and experimental, they have the potential to make DEM creation a lot less time consuming and labour intensive, by removing the need to manually remove artifacts.

You say:

" I sometimes think I'd like Manifold to allow me to manipulate surfaces with commands that forget the grid structure (i.e. as a matrix that I can transpose, solve, ...) - but R does that so welll there's no point in asking for it - and it already can do much of it it with transforms, just not how I like - I should point out. Manifold is about spatial data primarily, I suppose."

Performing grid math is not something I have to much experience of. However, I think there needs to tbe the ability to perform mathematical operations on single entire grids, multiple colocated entire grids and on a point by point basis within these grids according to whether the z-values at grid points meet a particular condition. If R does this well, please put in your ha' penny's worth.

Hope this clarifies a few issues, not complicates matters further.:-)

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

Hey Nick, if someone could add

1. pros and cons of the methods and typical situations to use them and

2. the current (and the coming) methods Manifold uses

to your list would - decorated with a mathematician - fill a gap in the background articels in manifolds help system.

mdsumner - 2004-10-23 4:44 AM

The big divisions I see are regular vs. irregular point data - everything else stems from that first issue.

I think thats true. As far as I can see - but I'm no expert - regular grids like STRM data are determine or processed as regular data points before they reach us. So there is no reality in a grid cell but the area the data point is representativ for.

Thus it should be allowed to use whatever methode suites to interpolate or smoothen sub-grids - as long as we know what happens. So I would like to know which method is used in Manifold to show height in the Status Bar. It changes before the cursor leaves the indicated cell. And more important: Which methode is used to transform different surfaces with different grid size.

To keep in mind the implications I use to carry the irregular points that made a surface throught the whole project as an all time companion of the surface. Surfaces based on irregular points in my projects would better be reflected as varitable 3D vector drawings. But of course I would like to be able to paste them as interpolated grid and level out some outliers or combine 3D-drawings with surfaces in a terrain. What makes up this new component is still to determine before it compiles to a wishlist item. Whereas surfaces based on irregular data points are a typical issue thought has to be given to other aspects like vertical walls or subsurfaces.

Klaus

KlausDE - 2004-10-23 5:26 PM ...if someone could add

1. pros and cons of the methods and typical situations to use them

To address this point in general terms:

The type of gridding method you should use depends on your objectives, what you want to produce, the data set being gridded and the time you are prepared to wait for the results (!). In general, you would use averaging interpolators when you are attempting to analyse longwave trends in z-value and where there may be significant shortwave z-value "noise" obscuring these trends. Conversely, exact interplators would normally be used where it is important to faithfully reproduce z-values at the locations of the original datapoints, althoug these values will not be recorded in the grid unless a gridpoint perfectly coincides with the location of a datapoint. For obvious reasons exact interpolators can be much more computationally intensive and time consuming than averaging interpolators.

An exception to the above is the nearest neighbour method, this is generally used to resample existing grids to create grids of different gridpoint spacing, orientation etc, when it is important to retain the exact values of the original gridpoints, without "bluring" the values through interpolation. This method is commonly used to resample imagery in remote sensing, as classification using NN resampled images is easier.

Another reason to have different methods is provide a choice of localised or global (=gride-wide) interpolation to determine gridpoint z-values. The former uses only the values around a gridpoint, the latter uses all the datapoints and produces a set of gridpoint z-values that are unique and dependant on the original datapoint values and their locations. Away from a gridpoint, the former is largely insensitive to the addition, subtraction of datapoints, change in their z-value or location. However, wrt to the global methods, adding, subtracting, or altering the z-value of a datapoint or its location, alters the interpolated z-values of ALL the grid points.

In general, it is a good policy to have a variety of gridding algorithms at ones disposal, so as to be able to compare the results from gridding using several. Then, and this is particulalry important for abstract datasets, one can assess whether a particular z-value trend or pattern is an artifact of the algorithm or real. Also, and this is especially important, where the object is to interpolate between unevenly distributed data, clumpy data, different algorithms can interpolate different values and trends across regions where there are no or few datapoints, it is then up to you to assess whether the results produced are acceptable or realistic.

In short wrt to which gridding algorithm you use, it is a case of "horses for courses". You should grid the data according to your objectives, using different methods, and then reject the results you like least.

KlausDE - 2004-10-23 5:26 PM As far as I can see - but I'm no expert - regular grids like STRM data are determine or processed as regular data points before they reach us. So there is no reality in a grid cell but the area the data point is representative for.

I think this goes beyond gridding methods. It is more an issue of how already-gridded datasets are interpreted. xy coordinates of SRTM gridpoints, for instance, within Manifold, are presently considered to define the location of the SW corner of a pixel 3 arcsends on a side (or 1 arcsecond on aside for the US). I suspect that it would be just as valid to consider the xy coordinates for a z-value within SRTM data as defining the centre point of a grid cell with dimensions 1- or 3-arseconds on a side. It is rather a matter or knowing one's data and how it was created.

I suspect Manifold designers chose to interpret grid point xy coordinates the way they have, because it makes for an easier life when coding as images can be created of the gridded data without the need for the representation of half gridcells at the edges of the image or quarter gridcells at image corners. For most uses with small grid cells this may be splitting proverbial hairs, but if the grid cells are large then it becomes non-trivial.

Thus it should be allowed to use whatever methode suites to interpolate or smoothen sub-grids - as long as we know what happens. So I would like to know which method is used in Manifold to show height in the Status Bar. It changes before the cursor leaves the indicated cell. And more important: Which methode is used to transform different surfaces with different grid size. To keep in mind the implications I use to carry the irregular points that made a surface throught the whole project as an all time companion of the surface. Surfaces based on irregular points in my projects would better be reflected as varitable 3D vector drawings. But of course I would like to be able to paste them as interpolated grid and level out some outliers or combine 3D-drawings with surfaces in a terrain. What makes up this new component is still to determine before it compiles to a wishlist item. Whereas surfaces based on irregular data points are a typical issue thought has to be given to other aspects like vertical walls or subsurfaces.

I dont understand this.

Cheers

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

Tomorrow, I will see if i can find some mathematical definitions, guidances on usage, and refrences about most of the algorithms listed above.

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

The Surfer 7 manual has this to say, by way of general recommendations, about the usage of the interpolation algorithms it supports:

The following list gives you a quick overview of each gridding method and some advantages and disadvantages in selecting one method over another.

Inverse Distance to a Power is fast but has the tendency to generate "bull's-eye" patterns of concentric contours around the data points. Inverse Distance to a Power does not extrapolate Z values beyond the range of data.

Kriging is one of the more flexible methods and is useful for gridding almost any type of data set. With most data sets, Kriging with the default linear variogram is quite effective. In general, we would most often recommend this method. Kriging is the default gridding method because it generates a good map for most data sets. For larger data sets, Kriging can be rather slow. Kriging can extrapolate grid values beyond your data's Z range.

Minimum Curvature generates smooth surfaces and is fast for most data sets but it can create high magnitude artifacts in areas of no data. The internal tension and boundary tension allow you control over the amount of smoothing. Minimum Curvature can extrapolate values beyond your data's Z range.

Natural Neighbor generates good contours from data sets containing dense data in some areas and sparse data in other areas. It does not generate data in areas without data. Natural Neighbor does not extrapolate Z grid values beyond the range of data.

Nearest Neighbor is useful for converting regularly spaced (or almost regularly spaced) XYZ data files to grid files. When your observations lie on a nearly complete grid with few missing holes, this method is useful for filling in the holes, or creating a grid file with the blanking value assigned to those locations where no data are present. Nearest Neighbor does not extrapolate Z grid values beyond the range of data.

Radial Basis Function is quite flexible. It compares to Kriging since it generates the best overall interpretations of most data sets. This method produces a result quite similar to Kriging.

Modified Shepard's Method is similar to Inverse Distance to a Power but does not tend to generate "bull's eye" patterns, especially when a smoothing factor is used. Modified Shepard's Method can extrapolate values beyond your data's Z range.

Triangulation with Linear Interpolation is fast. When you use small data sets, Triangulation with Linear Interpolation generates distinct triangular faces between data points. Triangulation with Linear Interpolation does not extrapolate Z values beyond the range of data.

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

If forum members want to see the results of gridding with the algorithms supported by Surfer 8, they can download the Surfer 8 demo (save and export disabled) from:

http://www.goldensoftware.com/demo.shtml

Further information on these methods should also be available via the Surfer 8 Help file accompanying the demo.

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

Nick Verge - 2004-10-23 10:04 PM

....I dont understand this.

Sorry, it seems my english needs more words than I used and then I was tricked by the unexpected appearance of a known bug.

Have you noticed that the status bar shows height of surfaces not in the raster that is displayed as image of the surface? I interpretet this as Manifold interpolating height for some sub-grid raster at least for display in the status bar if not for transforms and more. But your call back made me start a more thorough research: The status bar simply shows height shifted for halve a grid cell on the display. Same bug we have discussed in http://www.georeference.org/Forums/forums/thread-view.asp?tid=1143&posts=10

Sorry for adding confusion.

For all that have created a surface from known points it's important to notice, that readings of height from the status bar have to be corrected for a position (X-[local scale X]/2, Y-[local scale Y]/2) in v6.00 but the Transfer Height funcion returns the correct height of the grid point in the displayed grid cell without any interpolation or smoothing.

If we combine drawings and surfaces we have to correct the local offset. But for instance for mass calculation - based on the differenc between surface before and after some construction work - the height of the original grid point representativ for the area of a grid cell is used. The result is correct even in corse grids (as long as you keep away from the margins of our worlds slice).

So far the first point that was not clear.

Let me skip the next part that more or less deals with what you called "knowing one's data".

But at last and derived from the nature of the irregular grids I know I announced a risky thesis that I would like to discuss and may be it will breed a well founded contribution to the wishlist:

The thesis: The appropriate object model for an irregular grid of points in the three dimensions is a 3D vector drawing. This is able to keep the precision of these data in scaling and numeric processing and is able to show sharp edges, walls and subsurfaces.

Is the 3D drawing the right model for all types of irregular grided data? If we could transform them to the type of surfaces that we know this question could be left open.

hope it's clearer now and thanks for your contribution to methodology that especially non-geographers like me wellcome.

Klaus

Have sent a new feature request to Sales, for the gridding methods discussed above to be incorporated into a future version of Manifold, plus for any methods to permit control of datapoint search (search rules and asymetric searching) and datapoint weighting (anisotropy) supported by these algorithms, to be included also.

I encourage others to do the same.

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

Another web page worth looking at on this subject is

http://skagit.meas.ncsu.edu/~helena/gmslab/

One of the methods discussed here is RST (Relaxed Spline with Tension). This is not a new method and I've seen it used in code from years ago, but I suspect it shows up under different names. I don't think it is represented in your list as such, but apparently special cases of RST are mathematically equivalent to Krigging and (I think) minimum curvature. Perhaps the method should be added to the list since it is well documented and a generalized method.

I notice a few of the methods listed may already be supported in Manifold and would require a script to implement them more conveniently rather than new functionality. In particular I'm thinking of the nearest and natural neighbour methods which sound like spatial overlays of regularly spaced points on a Voroni tessalation. (Its been a while since I used Surfer - I don't recall the difference between nearest neighbour and TIN with linear interpolation.)

Cheers,

John Harrop

John

Funnily enough i was reading about RST recently and was about to put in a suggestion for it.

http://skagit.meas.ncsu.edu/~helena/gmslab/papers/MG-I-93.pdf

http://skagit.meas.ncsu.edu/~helena/gmslab/papers/MG-II-93.pdf

Regularised (not Relaxed) Spline with Tension

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

Nice link, lots of material for new feature requests there :-)

Re Relaxed Spline with Tension.

I think RST may be another name for the Minimum Curvature algorithm available in Surfer 7 and later. The Surfer 7 manual says that there are three factors that can be used to control MC interpolation, internal tension, boundary tension and relaxation. The relaxation factor controls how fast the algorithm converges (or whether it will converge at all!)

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

I just wanted to push this out of the archive again. It has been over 3 years since Nick tried to push for more algorithms to create gridded datasets. Altough Nick seems to be primarly concerned with the creation of terrain models, I also have an interest in gridded datasets, but related to point pattern/density analysis.

I am still wishing for Kernel Density and other spatial analysis methods, but development in this respect doesnt seem to have gotten much focus in the past 3 years.

More interestingly, a couple of people from the NHS (National Health Service) that attended the recent Manifold introduction training course have asked me about the capabilities of Manifold as a spatial point density analysis tool, and specifcally noticed the lack of relevant interpolators for this work.

Any more thoughts from the community on this subject, and maybe if we get lucky even some thoughts from the forum gods?

Check out the Manifold Wiki with SQL and scripting examples at http://www.manipedia.eu/

Spatial Knowledge, my personal blog.

Patrick,

it is good you have brought this thread to the fore again.

So there is no confusion, i am for as including the widest possible variety of interpolators, not just special ones for DEM creation.

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

Just wanted to add another great resource for spatial analysis in general, from which we can develop more focused feature requests:

http://www.spatialanalysisonline.com/

Check out the Manifold Wiki with SQL and scripting examples at http://www.manipedia.eu/

Spatial Knowledge, my personal blog.

Patrick'

and another...

http://spatial-analyst.net/index.php

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe

Thus far all the algorithms suggested above, are suitable for 2D interpolation of scaler fields, only. However, not all quantities that vary in space and time are scalars. Some are vectors (eg velocity, forces, topographic slope, dip of rock strata etc).

So what you might ask? why not simply independantly interpolate the magnitude and directional components of scattered data points of the vector field to create two sets of arrayed information describing the magnitude and direction of the field. Unfortunately, this cannot be done because if it is done, the resulting magnitude and direction values of the vector field at any point in the arraywould not be an unique solution!

The attached paper explains why, this is so. Although wind is the vector field under discussion therein, the principles apply equally to interpolation of scattered data point of other types of vector fields.

http://www.flame.org/~cdoswell/publications/Schaefer&Doswell_79.pdf

Dream no small dreams for they have no power to change the minds of men. - After Johann Wolfgang von Goethe