How to use pchip to interpolate between data points in cartesian coordinate format

I am trying to use pchip to interpolate between a series of data points. My data are spatial (lat/long) but have been converted into meters in a projected cartesian coordinate system (UTM 31N- i.e. x= 431804m, y=4571664m).

The code that I am using is below, but for some reason I get error messages when trying to create objects 't' and 'p'. long_m is a 1xcolumn matrix of longitude in meters, lat_m is a 1xcolumn matrix of latitude in meters, 0.01667= 1/60th in decimal, specifying that I want to end up 60 times more data points than I have now via the interpolation:

x=long_m;
y=lat_m;
t=x(1,1):0.016667:x(211,1);
p=pchip(x,y,t);

Output message:

t =
     Empty matrix: 1-by-0

p =
   Empty matrix: 1-by-0

I thought that the function would return a string of numbers that relate to my interpolated points, but instead I end up with an empty 't' matrix and an empty 'p' matrix. What am I doing wrong?

Some example data is below:

x= [518666 521872 519984 519591 518800];
y= [4694989 4667173 4644884 4645622 4647104];

I am new to matlab (although I have knowledge of programming language from R), so any advice would be greatly appreciated.

Answer 1

[Kelly Kearney]

The all-increasing or all-decreasing issue you're encountering isn't a restriction of pchip, but rather of the : operator. If you try to define a vector with a positive step size but decreasing values, you'll get an empty vector. For example:

 >> 5:1:2

 ans =

   Empty matrix: 1-by-0

I'm guessing that what you're really after in this case is a parametric interpolation, where you "fill in" points in both the x- and y-direction.

 % Some fake data

 theta = linspace(0, 2*pi, 10);
 x = cos(theta);
 y = sin(theta);

 % Interpolate

 told = 1:10;
 tnew = linspace(1,10,50);
 xnew = pchip(told, x, t);
 ynew = pchip(told, y, t);

 % Plot

 plot(x,y,'-b.', xnew, ynew, '-ro');

You can get more sophisticated to evenly space your final points, but this gives you there general idea.

[Rhiannon]

@ Kelly thanks for your suggestion. Yes, I am after a curvilinear interpolation method to 'fill in' points in both x and y. linspace() seems to be closer to doing what I want. Thanks! However, the code above returns an equal number of interpolated points between each original xy (i.e 5 new points between each original point). Is there an easy way to get pchip to space interpolated points equally within the entire point sequence (i.e. the first and last point in the sequence)?

[Kelly Kearney]

To approximate that, use the same idea as above, but define told based on distance along the curve:

 d = sqrt((x(2:end) - x(1:end-1)).^2 + (y(2:end) - y(1:end-1)).^2);
 told = [0 cumsum(d)];
 tnew = linspace(0, max(told), 50);

Note that the key to this method isn't the use of linspace as opposed to :; I could have gotten approximately the same results by saying

 tnew = 0:0.125:1;

The key is that I'm using a third variable ( t ) to represent "distance along the curve," and then interpolating both x and y against that.

[Rhiannon]

fantastic, many thanks!

Answer 2

[Shashank Prasanna]

This works perfectly fine for me:

>> x= [518666 521872 519984 519591 518800];
>> y= [4694989 4667173 4644884 4645622 4647104];
>> p=pchip(x,y,x)
p =

       4694989     4667173     4644884     4645622     4647104

p=pchip(x,y,x(1):0.01:x(end));

[Rhiannon]

@ Shashank thanks for your response. This example just gives me a 1x13401 matrix for p. Shouldn't using a step of 0.01 give a matrix with 400 rows if interpolating between 5 points? What I am after is an output matrix of interpolated x and y points (i.e. a 2xmatrix with x and y columns), so that I can then use these points to plot my interpolated data spatially. Can pchip not do this?

[Shashank Prasanna]

Your understanding of how the function works is not correct.

pchip first builds the model and lets you interpolate based on the pchip model.

If you have new data which is more finely defined:

x_new = x(1):0.01:x(end);
length(x_new)
ans =

         13401

p = pchip(x,y,x_new);

pchip will find all the corresponding points for y based on your n_new

Which means you x_new and p are the new x and new y which you can use for plotting.

The output matrix you are referring to is simply:

out = [x_new p];

The doc does a very good job at explaining this as well as provides you with good examples:

http://www.mathworks.com/help/matlab/ref/pchip.html

HTH

[Rhiannon]

Thanks for the clarification. If the subinterval (i.e. 0.01 in x_new) relates to the number of interpolated points between each pair of locations (as it seems to in the matlab example) then I still don't understand why length of x_new is 13401 rather than 400? Also, while the example that I've given above seems to work, other subsets of data that I have don't work. i.e.

x= [488461 487899 487486 485761 485760];
y= [4611769 4611769 4610523 4610173 4609059];
t = x(1):0.1:x(end);
p = pchip(x,y,t);

and give the same output message that I get with larger datasets:

t =

     Empty matrix: 1-by-0

p =

     Empty matrix: 1-by-0

Apologies if I am missing something very obvious here. Thanks.

[Shashank Prasanna]

x is decreasing in value.

x(1):0.1:x(end) makes no sense for decreasing values.

this does:

x(1):-0.1:x(end)

[Rhiannon]

Thanks. Using the minus value does work for some of the points, however I have a dataset that contains a mixture of decreasing and increasing values within the point series (the sequence of numbers in example 1 above shows this) and this seems to be the problem. Can pchip only handle a sequence of either all increasing, or all decreasing, values? If so, is there a different function for interpolation that can handle both. Thank you so much for your help, greatly appreciated!

[Shashank Prasanna]

Rhiannon, your Y can be increasing or decreasing. You can always sort your X which is meant to be increasing on the X axis. Your Y could be anything and PCHIP will happily interpolate for new X.

[sortX, idx] = sort(X);
sortY = Y(idx);

This way you can always take the sorted X to interpolate

newY = pchip(sortX,sortY,sortX(1):0.01:sortX(end));

And this doesn't change your data or results at all.

[Rhiannon]

Thanks Shashank for responding again. Unfortunately, I cannot sort my data into another order as data points are temporally related, which I now realise I didn't specify in my original question (apologies). However, linspace() seems to be a solution to the problem, so I'm giving that a go. Thanks for taking the time to help!