Filled contour with whiteout regions based on xyz data

I would like to produce a filled 2D contour plot with my xyz data, where x and y are the coordinates in metres with size of (1 x 61299) each and z is the variable. The scatter plot of x and y gives the following plot, (18000 x 18000) metres,

As recommended in https://uk.mathworks.com/matlabcentral/answers/373008-whiteout-regions-on-contour-plot and I have a fair bit of data, I have used delaunayTriangulation and nearestNeighbour to determine which points are too far from the scattered data and so should be set as NaN. My code for gridding is this:

[X,Y] = meshgrid(linspace(min(x), max(x),1000), linspace(min(y), max(y),1000)); % this should give a resolution of 18m ?
Z = griddata(x, y, z, X, Y);

DT = delaunayTriangulation(x, y);
[vi, d] = nearestNeighbor(DT, X, Y);

cutoff_distance = 500;
mask = d > cutoff_distance;
Z(mask) = nan;

contourf(X,Y,Z)
colorbar()

This gives the following contour plot.

I believe the problem is that the cutoff distance specified is shorter than the distances between the some of the points, so that the contour shows some incorrectly white out regions as compared to the scatter plot. Not to mention there is a very small white region (island) in the middle of the dense region. See this,

What are your recommendations to solve this problem? I have attached the data for your reference, where the columns are x, y and z.

[Image Analyst]

For the first problem, I'm not sure what the problem is. What would you expect if the problem weren't there? Maybe that will help me determine the source of the first problem.

For the second problem, the white islands/holes, that would be easy to solve if the data were in the form of a digital image. But if the region is just some coordinates that are not in your xyz data, then it's harder because we can't use functions like imfill or bwareaopen

[Raymond Lam]

Thanks for your reply, I would like to produce a contour plot while showing the boundaries/coastlines of the island, as shown in the large white regions of the scatter plot. The problem with my code is that there were regions in the contour plot that were incorrectly white out. Because the cut-off distance was too low for some of the distances between points (or nodes if you prefer to call it this way). This was shown in comparison between the contour and scatter plot.

Those islands (whether they are big or small) are coordinates that are not in my xyz data.

[Image Analyst]

Instead of delaunay would you consider a scatteredInterpolant? A demo is attached.

[Raymond Lam]

I have tried using scatteredInterpolant like this, which is similar to griddata. But this would produce a contour plot filling the white islands/regions because the z values of those coordinates at white regions are 0. That is why I tried to identify which points should be set to NaN using the above method (by setting a cutoff distance) Yet, it failed because the density of the nodes are not uniform.

[X,Y] = meshgrid(linspace(min(x), min(y), 1000), linspace(min(y), max(y), 1000));
F = scatteredInterpolant(x,y,z)
Z = F(X,Y);
contourf(X,Y,Z)

[Raymond Lam]

Perhaps others might have an alternative solution to this?...

Answer 1

[J. Alex Lee]

If you can tolerate manually identifying white-out regions (just by a single set of coordinates somewhere in the region, such as by ginput), then here is one way to do it, but it is not a very universal solution. Your data happens to look nice in that sampling is more dense near boundaries of the ostensible white-out areas. This suggests that if you make your markers large enough, you can "close" the white-out region boundaries as they squeeze togethe (as OP's images do).

So here's a script which attempts that, after acquiring white-out region seed coordinates by ginput and jsut hard-coding them as a list.

The steps are roughly:

1. find a scale/translation to move coordinates into rational image coordinates. I'm using ~1800x1800.
2. create an image where the pixels corresponding to coordinates where there is xyz data is 1 and 0 otherwise
3. dilate the image to "simulate" scatter plot with large marker size - dilation radius should be about half the max distance between points on the boundary of ostensible white-out region. I chose 15, which seems to work.
4. The complement of this mask are candidates of white-out regions - those regions that are not, are in between data points in regions that should not be whited out.
5. This is where you have to manually tag the regions (as by ginput+bwselect) that you think should be whited out. I have no idea where is supposed to be whited out and where is not, so i just chose wherever i thought.
6. Take the bwboundary of white-out regions on the image, and convert from image coordinates back to real coordinates. This is the "polygon" inside which any interpolated data from original data needs to be nan'd out using inpolygon.
7. define a meshgrid and use scatteredInterpolant to interpolate the xyz data onto the grid
8. find where grid coordinates are inpolygon() of the white-out regions defined in step 6, and nan the z data

I'm not sure if the inpolygon is the best way to go about this...it is super slow when the boundary to check is long, which is why I use a crude decimation technique to keep the script running in under the 55 second limit.

Also, this was a fun exercise, but the z data is 0 in most areas where you have ostensible white-out regions, so very uninteresting. If this is just an example data set out of many more where you have white-out regions in more interesting regions of data, this method might make more sense.

script_contour_whiteout

generating image from data took 0.39 sec
interpolating data on grid took 0.41 sec
ID-ing white-out areas from seeds took 2.38 sec
interpolating data on grid took 2.22 sec
whiting out data on interp'd grid took 9.35 sec

[Raymond Lam]

Thanks very much for your answer and help! I have tweaked the script a little bit to include the boundary of a small island in the middle. I am very satisfied the way it works although it did not capture/search some of the boundary edges (you will see when you zoom in the graph).

I also wonder if using contourf would be the best way to illustrate the variations of z-data with respect to x,y coordinates , e.g. my case. This dataset belongs to one of the timesteps (2881 steps). And it took me (and MATLAB) quite a bit of time to visualise the variations in time-series.

[J. Alex Lee]

@Raymond Lam, i think the "Best way to illustrate" depends...are the white-out regions always in the same place? do you actually even care about the regions where you white-out? (in this frame, your data is very boring near the white-out regions)

As a first pass, can you just use scatter and use the z-data for coloring and watch as the frames go by?

load("xyz.mat")
scatter(D(:,1),D(:,2),5,D(:,3),"filled")

[Raymond Lam]

Cheers Alex, the locations of white-out regions vary depending on the size of the domain such as below, where white regions represents the islands and black represents the scattered xy data with non-uniform density. White-out regions are important because they represents the coastlines of the islands.

I think the issue here is the difficulty in identifying the boundaries. Because (1) there is a large number of islands, (2) some islands are really small and (3) density of scattered xy data are non-uniform (some black-coloured areas further offshore would have a lower local data density than that near the coast).

I believe that would lead to inconsistencies in searching and identifying the boundaries with @J. Alex Lee 's script, and perhaps with other methods like delaunayTrigulation + nearestneighbour (https://uk.mathworks.com/matlabcentral/answers/373008-whiteout-regions-on-contour-plot) and boundary function (https://uk.mathworks.com/matlabcentral/answers/499716-generate-x-y-grid-coordinates-from-within-contour-boundary).

The results from using filled contour is what I would like my z-data to be. But I am yet to find a good way to do it and properly reflect the boundaries. Using scatter plot with large dot size would distort the actual view in high density regions.

What's your thought about this issue?

I have attached the x-y coordinates here for the above difgure. Feel free to make up the z-data if necessary.

Thanks a lot!

[J. Alex Lee]

For any data sets that share a common domain (island locations), you only need to detect boundaries once and you can re-use. Assuming you keep the same meshgrid query points, you can also just reuse the nan-masks.

Where did that B&W map come from? Note that my strategy is to approximate that B&W map by manually tagging regions and relying on a closed/dilated image from the dots to fill the regions. If you already have B&W maps, this part may be redundant and you can use your maps directly.

You would need to register the B&W map pixel coordinates onto your real xy coordinate locations somehow. Maybe you already know those coordinates, or others there is some manual work there, or some automated work if you have information in digital form somewhere.

You can use your maps directly as a mask by re-sampling your data exactly on the map's pixel grid and vice-versa - so the higher resolution map you have the better. In fact you could save a lot of time in my script by taking that strategy instead of using inpolygon (which could be useful for example if island positions were fixed, but for whatever reason you wanted to change your meshgrid query points every time)