How do I interpolate?
If one has irregular data scattered over a 2D plane e.g. x, y, z1,
z2, ... then there are two functions which can be used to interpolate
onto a regular grid - bilinear() and bivariate(). The main reason for
interpolating irregular data is so that we can contour it.
N.B. Interpolation is essentially an educated guess about what values lie between the actual data points. For this reason, points generated by interpolation may or may not reflect reality.. With Gsharp it would be fairly simple to make a nice contour plot of the average height of females in the UK based on just four sample points, but it would of course be close to useless. Gsharp includes kriging functions for those who need to know the amount of error introduced.
N.B. If your data is located in a regular fashion then there is no need to interpolate your data. The following plot shows an example of a contour of a complex grid which can be contoured without interpolation. Please read the article "How do I make a contour?".
So I am assuming that we have a number of datapoints (x1,y1), (x2,y2), with one or more z values for each location. Our datasets are called x, y, z1, z2, etc.
Step 1 - Select Grid
The first step is to decide the location of the grid we shall interpolate onto. If no grid is specified then Gsharp will choose a grid of 33 by 33 which has the same extremes as the input data (x and y). If you require these values to use as X Data and Y Data of a contour plot then you can use the following commands:
gridx = range(x,33); gridy = range(y,33);
N.B. It is possible to make a 2D contour of a plot without specifying the position of the grid using the X Data and Y Data resources, but I recommend specifying the grid location if you know it as this will allow you to set the domain limits to zoom in or out on your data and also to accurately overlay your contour plot with other isolines or scatter data.
N.B. The X Data and Y Data resources are used to specify the location of the grid (gridx and gridy). Specifying the input datasets (x and y) will just produce an error.
There are a number of commands which can be used to create datasets containing the location of the grid we are to interpolate onto.
gridx = range(x,20)
# 20 points, equally spaced, with the same limits as x
gridx = nicerange(x,20)
# as above but the limits of gridx are rounded to "nice" numbers
gridx = 1:20
# 1,2,3, ... ,20
There is no need for the grid to be evenly spaced. e.g.:
gridx = 1//2//2.25//2.5//2.75//3//4;
Step 2 - Select Method
The second decision is decide whether to use bilinear interpolation of bivariate.
- If you have more than 250 points or if your data contains faults then you must use bilinear.
- If you want to interpolate within a region such as a circle you can
either
- use bilinear() and pass the region as arguments or
- interpolate over the whole grid (with either method) and apply the region when you plot.
The advantage of only applying the region when we plot is that we can then contour right up to the edge of our region by setting the resource XuNregionBorderCells to true.
It is not possible to say which interpolation method is best. They are both algorithms and their suitability will depend on the form of the input data.
- Bivariate interpolation performs a bivariate interpolation and produces a grid (two-dimensional) dataset with regularly spaced points. The new surface goes through all the original data points.
- Bilinear interpolation approximates existing points onto a grid, interpolates by linear interpolation, and then improves the estimate by computing gradients at the points using quadratic interpolation. Finally, the values are refined by distance weighting methods. The new gridded surface will smooth the original data points, but will not necessarily pass through them.
N.B. If you are not happy with either of these methods then you can either
- Link in your own interpolation methods using UserFunc.c (the interpolation method does not have to be written in C).
- Write your own interpolation routine in GSL - this will execute more slowly than code written in C or FORTRAN.
Step 3 - Perform interpolation
grid = bilinear(x,y,z1,gridx,gridy); grid2 = bivariate(x,y,z1,gridx,gridy);
Step 4 - Plot Grid
A contour plot can be created like so:
create Viewport page_1.view; create Domain page_1.view.dom; create Graph page_1.view.dom.graph ( XuNgraphType = "2DContour", XuNcolorDataGrid = "grid", XuNxData = "gridx", XuNyData = "gridy" );
Interpolation Examples
The following data is used for all of these examples:
x = rnd(51); y = rnd(51); z1 = rnd(51);
gridx = nicerange(x,21); gridy = nicerange(y,21);
xbar = .1//.2//.3//.5//.7//.9; ybar = .1//.2//.5//.6//.8//.9;
th = (0:360)*pi/180; xreg = cos(th)/2+.5; yreg = sin(th)/2+.5;
Example 1 - Basic 2D Contour plot
|
Either method is used to interpolate from our datasets onto our interpolation grid. There are many options within Gsharp to configure this basic plot into something more appealing to the eye, but we avoid them here to keep the example simple.
|
 |
Example 2 - Basic 3D contour plot
|
|
Example 3 - Interpolating within a region
|
|
Example 4 - Interpolate within region and across fault (barrier)
|
|
Example 5 - Overlay isoline over contour
|
|
Example 6 - Overlay original points
 |
Gsharp allows you to have multiple graphs within a single domain, provided the graphs are of the same type. Contour plots and scatter plots are of a different tpye, so we must put each plot into a separate domain. We want to ensure that a position and colour in our second domain corresponds to a position and colour in the first domain and so we must ensure that:
|
grid1 = bivariate(x, y, z1, gridx, gridy); classlimits = range(z1,100); reset;
create Viewport page_1.view
( XuNyRatio = 1,
XuNxRatio = 1
);
create Domain page_1.view.dom
( XuNclassType = "limits",
XuNclassData = "classlimits"
);
create Graph page_1.view.dom.graph
( XuNgraphType = "2DContour",
XuNcolorDataGrid = "grid1",
XuNxData = "gridx",
XuNyData = "gridy",
XuNregionXData = "xreg",
XuNregionYData = "yreg",
XuNregionBorderCells = "contoured"
);
recalc; create Domain page_1.view.dom2 ( XuNclassType = "limits", XuNclassData = "classlimits", XuNxMinimum = page_1.view.dom.XuNxMinimum, XuNxMaximum = page_1.view.dom.XuNxMaximum, XuNyMinimum = page_1.view.dom.XuNyMinimum, XuNyMaximum = page_1.view.dom.XuNyMaximum ); create Graph page_1.view.dom2.graph ( XuNgraphType = "scatter", XuNxData = "x", XuNyData = "y", XuNcolorData = "z1" ); |