Asked by Ivar Teeseling
on 16 Sep 2018

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 by jonas
on 16 Sep 2018

Edited by jonas
on 17 Sep 2018

Accepted Answer

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
on 17 Sep 2018

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
on 17 Sep 2018

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.