calculate conservative interpolation of two vectors in matlab

Question

G'day

Firstly, apologies for poor wording - I'm at a bit of a loss of how to describe this problem. I'm trying to calculate the conservative interpolation between two different vertical coordinate systems.

I have a vector of ocean transport values Ts, that describe the amount of transport at different depth values S. These depths are unevenly spaced (and size(S) is equal to size(Ts)+1 as the values in S are the depths at the top and bottom over which the transport value applies). I want to interpolate(/project?) this onto a vector of regularly spaced depths Z, where each new transport value Tz is formed from the values of Ts but weighted by the amount of overlap.

I've drawn a picture of what I mean (sorry for the bad quality webcam picture) I want to go from Ts1,Ts2.Ts3...TsN (bottom lines) to Tz1,Tz2,...TzN (top lines). The locations in the x direction for these are s0,s1,s2,...sN and z0,z1,z2,...zN. An example of the 'weighted overlap' would be:

Tz1 = a/(s1-s0) Ts1 + b/(s2-s1) Ts2 + c/(s3-s2) Ts3

where a, b and c are shown in the image as the length of overlap.

Some more details:

Example of z and s follow:

z = 0:5:720;
s = [222.69;...
   223.74
   225.67
   228.53
   232.39
   237.35
   243.56
   251.17
   260.41
    271.5
   284.73
   300.42
    318.9
   340.54
   365.69
   394.69
   427.78
   465.11
   506.62
   551.98
   600.54
    651.2];

Note that I'm free to define z, but not s. Typically, z will be bigger than s (i.e. the smallest value in z will be smaller than in s, while the largest value in z will be larger than in s).

Help or tips greatly appreciated. Cheers, Dave

Answer 1

I don't think there is an easy solution, as stated in the comments. I'll give it a go though :

One hypothesis first : We assume z0>s0 in order for your problem to be defined.

The idea (for your example) would be to get to the array below :

1  (s1-z0)  s1-s0   Ts1
1  (s2-s1)  s2-s1   Ts2
1  (z1-s2)  s3-s2   Ts3
2  (s3-z1)  s3-s2   Ts3
2  (z2-s3)  s4-s3   Ts4
3  (z3-z2)  s4-s3   Ts4
      ......

Then we would be able to compute, for each row : column1*column3/column2 and then use accumarray to sum the results with respect to the indexes in the first column.

Now the hardest part is to get this array :

Suppose you have :

• A Nx1 vectors Ts

• 2 (N+1)x1 vectors s and z, with z(1)>s(1).

  Vectsz=sort([s(2:end);z]); % Sorted vector of s and z values

In your case this vector should look like :

z0
s1
s2
z1
s3
z2
z3
...

The first column will serve as a subscript to apply accumarray, so we'll want it to increase each time there is a z value in our vector Vectsz

First=interp1(z,1:length(z),Vectsz,'previous');

Second=[diff(Vectsz);0]; % Padded with a 0 to keep the right size

Temp=diff(s);
Third=interp1(s(1:end-1),Temp,Vectsz,'previous');

This will just repeat the diff value everytime you have a z value in your vector Vectsz.

The last column is built exactly like the third one

Fourth=interp1(s(1:end-1),Ts,Vectsz,'previous');

Now that the array is built, a call to accumarray is enough to get the final result :

Res=accumarray(First,Second.*Fourth./Third);

EDIT : There is actually no need for the use of interp1 with the previous option :

Vectsz=sort([s(2:end);z]);

First=cumsum(ismember(Vectsz,z));

Second=[diff(Vectsz);0];

idx=cumsum(ismember(Vectsz,s(2:end)))+1;

Diffs=[diff(s);0];

Third=Diffs(idx);

Fourth=Ts(idx);

Res=accumarray(First,Second.*Fourth./Third);