Asked by Ivar Teeseling
on 16 Sep 2018 at 17:38

**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 at 23:26

Edited by jonas
on 17 Sep 2018 at 7:25

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 at 7:24

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 at 8:11

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.