Filled contour with whiteout regions based on xyz data

9 views (last 30 days)
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 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;
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.
Raymond Lam
Raymond Lam on 29 Dec 2022
Perhaps others might have an alternative solution to this?...

Sign in to comment.

Answers (1)

J. Alex Lee
J. Alex Lee on 3 Jan 2023
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 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.
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
J. Alex Lee
J. Alex Lee on 6 Jan 2023
Edited: J. Alex Lee on 6 Jan 2023
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)

Sign in to comment.


Find more on Contour Plots in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!