Thread Subject: Interpolation with NaNs

Hello,

Is there a way to get 2 dimensional interpolation to ignore certain
values?

I will outline my problem first as it's the easiest way to make myself
understood. I have a dataset of surface values for points in the world
ocean. Currently it is at 1 degree resolution but I want to
interpolate to 1/10 degree resolution using cubic interpolation. The
input array is 360x180 points (longitude x latitude) and points that
are on land are given as NaNs.

When I use the interpolation as specified the NaNs propagate out and
areas that were ocean are now being identified as land. I can't just
take ocean values, put them in a new 'ocean' array and then
interpolate, because this would cause interpolation between points
that would ordinarily be separated by land.

I need a way for matlab to consider NaNs as hard boundaries and to
ignore interpolation across them. Is there a way to do this? Or will I
have to code the interpolation algorithm from scratch?

I hope I've made myself clear. If not I'm happy to clarify further

Thanks,
Rob

CSM <robthevagrant@gmail.com> wrote in message
<e143d594-200b-41fb-b3d5-
ff856de1d93c@s8g2000prg.googlegroups.com>...
> Hello,
>
> Is there a way to get 2 dimensional interpolation to ignore certain
> values?
>
> I will outline my problem first as it's the easiest way to make myself
> understood. I have a dataset of surface values for points in the world
> ocean. Currently it is at 1 degree resolution but I want to
> interpolate to 1/10 degree resolution using cubic interpolation. The
> input array is 360x180 points (longitude x latitude) and points that
> are on land are given as NaNs.
>
> When I use the interpolation as specified the NaNs propagate out and
> areas that were ocean are now being identified as land. I can't just
> take ocean values, put them in a new 'ocean' array and then
> interpolate, because this would cause interpolation between points
> that would ordinarily be separated by land.
>
> I need a way for matlab to consider NaNs as hard boundaries and to
> ignore interpolation across them. Is there a way to do this? Or will I
> have to code the interpolation algorithm from scratch?
>
> I hope I've made myself clear. If not I'm happy to clarify further
>
> Thanks,
> Rob
------------
  The interpolation-extrapolation process depends heavily on an assumption
of continuity of some sort when increasing the resolution of an array of data.
However, when discontinuities are known to be present, such as the
transition from ocean to land in your case, an increase in resolution becomes
impossible in the vicinity of that discontinuity. The necessary information is
simply not there in the data.

  Specifically, if you have two array points a degree apart with one being an
ocean point and the other a land point, there is no way of determining to a
resolution of a tenth of a degree where the transition from ocean to land
occurs between those two points. The same difficulty is present in a two-
dimensional representation of a discontinuity boundary. Your shore line
could be drawn at a higher resolution in a great many different ways, all
consistent with the given lower resolution data. There is no way of
distinguishing which of these would be the most accurate without some
additional information being made available.

Roger Stafford

Roger Stafford <ellieandrogerxyzzy@mindspring.com.invalid> wrote:
> CSM <robthevagrant@gmail.com> wrote in message
> <e143d594-200b-41fb-b3d5-
> ff856de1d93c@s8g2000prg.googlegroups.com>...
>> Hello,
>>
>> Is there a way to get 2 dimensional interpolation to ignore certain
>> values?
>>
>> I will outline my problem first as it's the easiest way to make myself
>> understood. I have a dataset of surface values for points in the world
>> ocean. Currently it is at 1 degree resolution but I want to
>> interpolate to 1/10 degree resolution using cubic interpolation. The
>> input array is 360x180 points (longitude x latitude) and points that
>> are on land are given as NaNs.
>>
[snip]

> The interpolation-extrapolation process depends heavily on an assumption
> of continuity of some sort when increasing the resolution of an array of
data.
> However, when discontinuities are known to be present, such as the
> transition from ocean to land in your case, an increase in resolution
becomes
> impossible in the vicinity of that discontinuity. The necessary
information is
> simply not there in the data.
>
[snip]

>
> Roger Stafford
>

Roger's comments are spot on here, which is why you don't find a simple
answer with the standard versions of the interpolation routines.

You might find this easiest if done in a few steps.

First, make note of where your NaNs are (find(isnan(x)), then replace them
all with zeros, -Inf, or their nearest neighbours. In one dimension I
would do this by something like the following:

ii = (1:length(x))'; % An index vector
ff = find(isnan(x)); % These are NaNs.
ff2 = find(~isnan(x)); % These are not
xa = x(ff2); % Without the NaNs.
iia = ii(ff2);
xn = interp1(iia,xa,ii,'nearest');

Now do your increased resolution interp:
ii2 = (1:.1:length(x))';
x2 = interp1(ii,xn,ii2,'pchip');

Next you need to replace the land areas with NaNs again.
xnan = zeros(size(ii));
xnan(ff) = nan;
x3 = interp1(ii,xnan,ii2,'nearest');

x4 = x2 + x3;

--
Dr Tristram J. Scott | 44 Montague Road
Energy Consultant | Cambridge, CB4 1BX
                                   | England
tristram.scott@quantmodels.co.uk | ph (+44 1223) 526 255

> The interpolation-extrapolation process depends heavily
on an assumption
> of continuity of some sort when increasing the resolution
of an array of data.
> However, when discontinuities are known to be present,
such as the
> transition from ocean to land in your case, an increase in
resolution becomes
> impossible in the vicinity of that discontinuity. The
necessary information is
> simply not there in the data.
>
> Specifically, if you have two array points a degree
apart with one being an
> ocean point and the other a land point, there is no way of
determining to a
> resolution of a tenth of a degree where the transition
from ocean to land
> occurs between those two points. The same difficulty is
present in a two-
> dimensional representation of a discontinuity boundary.
Your shore line
> could be drawn at a higher resolution in a great many
different ways, all
> consistent with the given lower resolution data. There is
no way of
> distinguishing which of these would be the most accurate
without some
> additional information being made available.
>
> Roger Stafford


I agree with Roger's comments as well. What you might try
is using an interpolation on a window around the points of
interest. You could determine the region of geological
relevance (say 10 miles for example) and then apply the
interpolation considering only that region. It won't help
you on the land-adjacent regions, but its better than nothing.

One of Lauren's recent blogs described a vectorized way to
do something similar with brain data.

Hope that helps

CSM <robthevagrant@gmail.com> wrote in message
<e143d594-200b-41fb-b3d5-ff856de1d93c@s8g2000prg.googlegroups.com>...
> Hello,
>
> Is there a way to get 2 dimensional interpolation to
ignore certain
> values?
>
> I will outline my problem first as it's the easiest way to
make myself
> understood. I have a dataset of surface values for points
in the world
> ocean. Currently it is at 1 degree resolution but I want to
> interpolate to 1/10 degree resolution using cubic
interpolation. The
> input array is 360x180 points (longitude x latitude) and
points that
> are on land are given as NaNs.
>
> When I use the interpolation as specified the NaNs
propagate out and
> areas that were ocean are now being identified as land. I
can't just
> take ocean values, put them in a new 'ocean' array and then
> interpolate, because this would cause interpolation
between points
> that would ordinarily be separated by land.
>
> I need a way for matlab to consider NaNs as hard
boundaries and to
> ignore interpolation across them. Is there a way to do
this? Or will I
> have to code the interpolation algorithm from scratch?
>
> I hope I've made myself clear. If not I'm happy to clarify
further
>
If you reshape your 2-D structures to 1-D and remove the
nans, then you can use griddata to interpolate without
painting in the land areas. Maybe you'd have to cut up your
2-D structures into several smaller areas if griddata runs
out of memory, if so do it with some overlap at the borders.
Then you'd have to insert nans again where there should be
land. I implicitly assume that you're not perticularly
interested in areas near land.

HTH,
Bjeorn

Please try the bellow.

Bruno

%
% Generate data
%
[img]=peaks(20);
figure(1); clf;
imagesc(img);
[ny nx]=size(img);
[X Y]=meshgrid(1:nx,1:ny);
disp('Please enter a land (polygon), <CR> when finish');
[xpoly ypoly]=ginput;
img(inpolygon(X,Y,xpoly,ypoly))=NaN;
imagesc(img);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Expand image by 1 pixel
expand=[nan(1,nx+2); ...
        nan(ny,1) img nan(ny,1); ...
        nan(1,nx+2)];
mask=double(~isnan(expand));
expand(isnan(expand))=0; % do not at the NaN

imgconv=expand(1:end-2,2:end-1)+expand(3:end,2:end-1)+...
        expand(2:end-1,1:end-2)+expand(2:end-1,3:end);

maskconv=mask(1:end-2,2:end-1)+mask(3:end,2:end-1)+...
    mask(2:end-1,1:end-2)+mask(2:end-1,3:end);

% closure nan in the original image by 4 neighboring
imgclosure=imgconv./maskconv;

% Generate finer grid
rg=4; % resolution gain
[xi yi]=meshgrid(linspace(1,nx,nx*rg),linspace(1,nx,ny*rg));
% Interpolation
imginterp=interp2(1:nx,1:ny,imgclosure,xi,yi);

%
% Restore the land as it is
%
[inan jnan]=find(isnan(img));
[inan k]=ndgrid(rg*(inan-1),repmat(1:rg,1,rg));
inan=inan+k;
inan=inan(:);
[jnan k]=ndgrid(rg*(jnan-1),repmat(1:rg,rg,1));
jnan=jnan+k;
jnan=jnan(:);
imginterp(sub2ind(size(imginterp),inan,jnan))=NaN;

figure(2); clf;
imagesc(imginterp);
                         

Opps! Sorry I find a bug (forget to restore unconvoluted
data before interpolation). This should be better.

Bruno


%
% Generate data
%
[img]=peaks(20);
figure(1); clf;
imagesc(img);
[ny nx]=size(img);
[X Y]=meshgrid(1:nx,1:ny);
disp('Please enter a land (polygon), <CR> when finish');
[xpoly ypoly]=ginput;
img(inpolygon(X,Y,xpoly,ypoly))=NaN;
imagesc(img);
[inan jnan]=find(isnan(img));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Expand image by 1 pixel
expand=[nan(1,nx+2); ...
        nan(ny,1) img nan(ny,1); ...
        nan(1,nx+2)];
mask=double(~isnan(expand));
expand(isnan(expand))=0; % do not at the NaN

% neighbouring indices
i0=sub2ind(size(expand),inan+1,jnan+1);
i1=i0+1;
i2=i0-1;
i3=i0+size(expand,1);
i4=i0-size(expand,1);

imgconv=expand(i1)+expand(i2)+expand(i3)+expand(i4);
maskconv=mask(i1)+mask(i2)+mask(i3)+mask(i4);

% closure nan in the original image by 4 neighboring
imgclosure=imgconv./maskconv;
img(sub2ind(size(img),inan,jnan))=imgclosure;

% Generate finer grid
rg=8; % resolution gain
[xi yi]=meshgrid(linspace(1,nx,nx*rg),linspace(1,nx,ny*rg));
% Interpolation
imginterp=interp2(1:nx,1:ny,img,xi,yi,'spline');

%
% Restore the land as it is
%
[inan k]=ndgrid(rg*(inan-1),repmat(1:rg,1,rg));
inan=inan+k;
[jnan k]=ndgrid(rg*(jnan-1),repmat(1:rg,rg,1));
jnan=jnan+k;
imginterp(sub2ind(size(imginterp),inan(:),jnan(:)))=NaN;

figure(2); clf;
imagesc(imginterp);