Interpolating 2D data including gaps (rawZ = 0) with GRIDDATA(), but without filling in the data gaps

Background:

• I am trying to project scattered 2D raw data (rawX, rawY, rawZ) onto a 2D grid (intX, intY) using GRIDDATA()
• The scattered 2D raw data has a data gap where no measurements have been made (rawZ = 0), as shown in the figure
• When I use GRIDDATA(), the data gap gets automatically filled in by the interpolation operation (see figure), but I want the data gap to exist in the interpolated data as well
• intX, intY are created using MESHGRID()

Tried Methods

• Setting rawX and/or rawY = NaN, but this does not help
• I do ignore the gapped data points (rawZ == 0) doing the interpolation as the zero values will locally distort the interpolated values

Question

• Is there method to directly have GRIDDATA() ignore the rawZ data gap when creating intZ?

Much appreciated!!

I have attached the code below:

% Remove data gaps for interpolation
mask = (rawZ==0);
rawX_adj = rawX(~mask);
rawY_adj = raw(~mask);
rawZ_adj = raw(~mask);

% Create the interpolation mesh grid (intX, intY)
intX = linspace(min(rawX_adj(:)), max(rawX_adj(:)), 100);
intY = linspace(min(rawY_adj(:)), max(rawY_adj(:)), 100);
[intX, intY] = meshgrid(intX, intY);

% Project the raw data onto the (intX, intY) grid
intZ = nangriddata(rawX_adj, rawY_adj, rawZ_adj, intX, intY);

Answer 1

[jonas]

The nature of your data is a bit unclear. You say scattered data but that the hole is associated with a value of zero. If you have coordinates for the hole, then Id just extract those, define a polygon around the hole using boundary and remove all points inside the polygon after interpolation using inpolygon .

Some untested code:

%%Find coordinates inside hole
xh=rawX((rawZ==0)):
yh=rawY((rawZ==0)):

%%Define bounding polygon
p=boundary(xh,yh);

%%Interpolate
[X,Y]=meshgrid(resx,resy);
Z=griddata(rawX,rawY,rawZ,X,Y)

%%Find and remove pts inside polygon and plot
in=inpolygon(X,Y,xh(p),yh(p))
Z(in)=NaN;
surf(X,Y,Z)

Thats simple. The following two solutions are for when the holes are defined by their lack of data:

Technically, there are gaps between every one of your data points and you want to remove large gaps. There is, to my knowledge, no option in griddata to define the "largest distance" over which it is allowed to interpolate. Here are however two options for you where you first interpolate and then remove the data where the holes were.

Method 1:

Create a grid for interpolation and find the closest distance between each point in your grid and a data point using dsearchn .

%%Make some data
x=rand(1,1e4)';
y=rand(1,1e4)';
z=x.*y;
xp=[0.4 0.6 0.6 0.4 0.4];
yp=[0.4 0.4 0.6 0.6 0.4];
in=inpolygon(x,y,xp,yp);
x(in)=[];
y(in)=[];
z(in)=[];

%%Plot data
subplot(1,3,1)
plot(x,y,'.')

%%Desired grid
[X,Y]=meshgrid(0:0.01:1,0:0.01:1);

%%Find holes
[a b]=dsearchn([x,y],[X(:),Y(:)])
subplot(1,3,2)
scatter(X(:),Y(:),[],b)
id=a(b>.02)

You clearly see the hole in the right figure. Now you can just grab the indices (using some threshold value) and set those points to NaN after interpolating onto your grid using griddata.

Method 2 (faster)

Use hist3 to find bins with missing data. Define a polygon around those points and remove points inside of the polygon.

%%Same data as previous example
%%Desired grid and interp
[X,Y]=meshgrid(0:0.01:1,0:0.01:1);
Z = griddata(x,y,z,X,Y);

%%Bin count, find empty bins and their bounds
[count,bins]=hist3([x,y],'ctrs',{0:0.03:1,0:0.03:1},'cdatamode','auto');
bins=[bins{1}(:) bins{2}(:)];
[idx,idy]=find(count==0);
xh=bins(idx,2);
yh=bins(idy,1);
p=boundary(xh,yh);

%%Plot data with boundary
subplot(1,4,2)
plot(x,y,'.',...
    xh(p),yh(p),'r');

%%Find grid points inside of polygon and remove
in=inpolygon(X,Y,xh(p),yh(p))
Z(in)=NaN;

%%Plot surface with ole
subplot(1,4,3)
surf(X,Y,Z,'edgecolor','none')

Note that the hole is not completely filled. You can fill a larger portion of the hole by reducing the bin size. However, with smaller bins you may also get other bins with holes, depending on the nature of your data. If this happens, then you have to remove outliers before running boundary, else you will get a very oddly shaped polygon.

[Ivar Teeseling]

Dear Jonas,

First of all; thank you for the very useful and exhaustive answer. I understand the confusion and should have been more clear:

The data 'gaps' are the locations where no measurements were made, which are marked by the value of zero for the measured data parameter (rawZ). However, I do know their coordinates in the raw data set (rawX, rawY). In that sense, and according to its strict definition, it is not a true data gap as I do have the corresponding locations associated with (rawZ = 0) values.

About your solutions:

Method 1

I was unaware of the functions Boundary() and inpolygon(). They look awesome and will be very easy to implement. Does exactly the trick!

Method 2

This was indeed what I was trying to do as I hoped there would be a way to tell griddata() when a gap was too large for interpolation. It is able to determine for itself as increasing the size of the hole will at some point stop griddata() from interpolating across the hole.

This method will work fine as well, but as I know the (rawX, rawY) locations of the 'holes', I can more easily work with method 1.

Great answer. Thanks!

I will make sure to upload the final section of code for future readers.

[jonas]

My pleasure! I've been thinking of ways of doing this myself, but never found the motivation to actually write the code. That's why the answer is a bit more comprehensive than necessary ;)