Contour plot ignores a set of data points

Hello,

I'm constructing a contour plot using location (longitude and latitude) and respective depth data. However, the output contour plot seems to ignore (or interpolate in a way that is not intended) a set of data points. This makes the result undesirable. Is there any thing that I can do to fix this?

Following is the piece of code that I use to plot the contour plot.

% reading the bathy data
data = xlsread('datafile.xlsx');
lon = data(:,1);                    % Longitude
lat = data(:,2);                    % Latitude
alt = data(:,3);                    % Altitude

[xq,yq] = meshgrid(min(lon):.001:max(lon), min(lat):.001:max(lat));
[x,y,z] = griddata(lon,lat,alt,xq,yq, 'cubic');

contourf(x,y,z);

I have attached the figues for reference.

• Data_Points.png (Locations of the data points)
• Overlay.png (Result contour plot is overlaid on Google Earth)
• Overlay2.png (Result contourf plot)

[Subhadeep Koley]

@ James Can you share the datafile.xlsx file?

[Jake]

I know this is going to sound bad, but is it possible to have any pointers or advice without the datafile.xlsx?

I'm sorry, but the bathymetry data in this location is scarce and the project that I'm on will not allow me to share the data :( I'm really sorry.

[John D'Errico]

This is just some random conjecture on my part, but it seems from the pictures you post as if you may not want the contours to cross over irregularly shaped land regions that intrude into the water. As one might expect, such a boundary will be highly complex and irregular, certainly not convex in any respect. And all you seem to have is the point set, which indicates only information sampled where there is water.

If I am correct, then it is not that the contour plot ignores the points provided, but that the contour plot has no hint that the points exist. You are not contouring data in fact. The contour plot is applied to the result of a griddata operation ON your data. contourf sees only what griddata produces.

The problem of course is then you used meshgrid to populate the region with a grid, then interpolated using griddata. The real problem is then griddata, which has no clue that a void between two points is not valid to interpolate across, whereas another void is perfectly valid to cross. Note that EVERY pair of points has a void between them. Points are just points after all, scattered data. If some of the data is slightly more scattered, what is the griddata algoithm to know?

As such, I think (think - but this is only my perhaps weak conjecture) is you are asking for a way to force the contours to live ONLY inside a polygonal region describing the shoreline. Griddata has no such concept of course, and that makes what you tried to do fail.

If so, it is theoretically possible to solve such a problem. The scheme I would follow would have me triangulate the polygonal shoreline region, then using the data which lives only inside that polygonal region, use that data to infer a surface on that triangulation. Once you now have z(x,y) over a triangulated region, you can interpolate to produce contours that live only within the shoreline polygon. At least that would be my algorithmic approach. However, I'd not be amazed if perhaps the mapping toolbox already has some existing capability to solve this problem too.

At the very least, one simpler option might be to use what griddata produces, then force the grid points outside the polygonal region of interest to take on NaN values. inpolygon or polyshapes can offer the necessary tests. This will give you a slightly coarse approximation to the shoreline, and the contours will not cross through a region of NaNs.

[Jake]

Absolutely, yes! That is the exact approach I need and my apologies for not putting it out with the right words. I'm still learning MATLAB (Alone, I must add) and the approach that you have mentioned seems quite complex.

"The scheme I would follow would have me triangulate the polygonal shoreline region, then using the data which lives only inside that polygonal region, use that data to infer a surface on that triangulation."

How can I approach this?

I'd really appreciate any advice or suggestion. TIA!

[John D'Errico]

Oh geez. I had to stick my neck out. I had to stick my nose into this. ;-) Well, first, I would not even consider suggesting a re-creation of the wheel. Does the mapping toolbox have such capability, or something that would be sufficient? If so, then you never want to do work twice. And I lack the mapping TB, so I have no clue in this respect.

So hoping to avoid some effort here, does the mapping TB have what you can use?

[Jake]

I'm sorry for dragging you back after your explaination, but I'm stuck and so far I've had no luck :(

And to make things worse, I don't have the mapping toolbox either.

I'd search for more options and would totally understand if you avoid the problem altogether :) Your explaination alone gave me a fair idea of what is to be done, although I don't have a solid idea how to.

[Bjorn Gustavsson]

The easiest way (if there is no mapping-toolbox methods that allow you to define a perimeter on the tri-scatter-interpolation), would be for you to determine the points that lie along the perimeter of your data-set. Then you can "simply" use the coordinates of those points to determine which points are inside and outside the boundary using inpolygon as John suggested. The inpolygon call migh be a bit time-consuming, but for any set of coordinates from meshgrid you only need to do this once. Then you get away from the spurious contours by setting everything outside to nan.

[John D'Errico]

A very simple example of what I described is shown here:

px = [0, 1, .5, .49, 0];
py = [0, 0, sqrt(3), sqrt(3)-.5, 0];
ps = polyshape(px,py);
[x,y] = meshgrid(linspace(0,1,100),linspace(0,1.8,200));
z = x.^2 + y.^2;
inp = ps.isinterior(x(:),y(:));
z(~inp) = NaN;
H = contourf(x,y,z);
hold on
HB = plot(ps);
HB.FaceColor = 'none';

The edges are coarsely stairstepped, since I used only a fairly coarse mesh. Also note that the result lives ENTIRELY inside the polygonal region. So the steps intrude into the polyshape, because contourf will not contour any part of the grid that even touches a NaN. A finer mesh would reduce the stairstepping artifacts.

I defined z there using a simple computation, but nothing would have prevented me from using a tool like griddata or scatteredInterpolant if I had real data.

Is this a bit of a hack? Well, sadly, yes. That is because griddata still is using data from acrsss a hole to infer the shape of that surface in some region. I just deleted what griddata did in those regions, so contourf did not see the stuff I wanted hidden. The point is, griddata is still doing what is technically a bad thing, because it does not understand what you want. And that MAY introduce artifacts into the result. (Carefully consider this issue, as it may easily produce non-obvious artifacts. I can give an example of where this would be a problem.)

However, if you can live with this simplistic solution, then it may be a good choice for you.

[Jake]

@ Bjorn Thanks a lot! However, my problem would be identifying the points that lie along the perimeter of my data-set. Because I have thousands of Longitude/Latitude and corresponding Altitude data :( I will have to find the points manually on a Google Earth scatter map, I think.

This would be one of the best approaches, if I can't find another solution :) Thank you so much for your time.

[Jake]

@ John: This looks perfect! But as I commented above, I have to find the coordinates of px and py manuallu from a data set of thousands of lon/lat values :) But, like both of you have mentioned, this would be an easier approach and I'm definitely trying this now.

I won't drag you more into this because that would be a huge effort on your part. I'd accept the answers if you comment them seperately :)

You've been a huge help. Thanks again!

[Bjorn Gustavsson]

James, you might get away with generating a perimeter manually selecting points along a perimeter. Once you've plotted your points like in the Data_Points.PNG image you can manually select points along a perimeter (with the accuracy of your chosing) using ginput. If that becomes too tedious you might find this file-exchange contribution useful: snake-algorithm, some of the alpha-shape contributions, or active contour contributions.

[Jake]

Honestly, I had no idea about these approaches. Thanks a lot, and this is going to take some time for me to learn and play around because I have not tried these before :)

I will comment here again with my progress. Thanks again, Bjorn! You made the day!

Answer 1

[John D'Errico]

I've moved what was a comment into an answer, then added a bit...

A very simple example of what I described is shown here:

px = [0, 1, .5, .49, 0];
py = [0, 0, sqrt(3), sqrt(3)-.5, 0];
ps = polyshape(px,py);
[x,y] = meshgrid(linspace(0,1,100),linspace(0,1.8,200));
z = x.^2 + y.^2;
inp = ps.isinterior(x(:),y(:));
z(~inp) = NaN;
H = contourf(x,y,z);
hold on
HB = plot(ps);
HB.FaceColor = 'none';

The edges are coarsely stairstepped, since I used only a fairly coarse mesh. Also note that the result lives ENTIRELY inside the polygonal region. So the steps intrude into the polyshape, because contourf will not contour any part of the grid that even touches a NaN. A finer mesh would reduce the stairstepping artifacts.

I defined z there using a simple computation, but nothing would have prevented me from using a tool like griddata or scatteredInterpolant if I had real data.

Is this a bit of a hack? Well, sadly, yes. That is because griddata still is using data from acrsss a hole to infer the shape of that surface in some region. I just deleted what griddata did in those regions, so contourf did not see the stuff I wanted hidden. The point is, griddata is still doing what is technically a bad thing, because it does not understand what you want. And that MAY introduce artifacts into the result.

As I said, artifacts might arise. Consider a randomly u-shaped region. In this case, the region was just based on some quick clicks of the mouse using ginput. Then I created a polyshape, based on the perimeter polygon, with some randomly generated data inside of it.

My fear is that what happens in one wing of that u-shaped region can cross the peninsula between those wings to influence the shape in the other wing. This could be made more significant due to the use of cubic interpolation in griddata. That would not be true in the case where the surface is constructed from a triangulation that lives inside only the polygon. At the same time, the only tools I have to build that surface on a general non-convex triangulated region would require some tools that are now fairly outdated.

The final issue mentioned is how to extract the polygon itself, if all you have are the scattered points. This is somewhat difficult to do via an automatic scheme, because all scattered data has holes between EVERY pair of points. It seems as if the data points around the perimeter were fairly closely spaced. So you might try to use some sort of intelligent scheme to build the polygon. Personally, unless you are going to do this more than a few times, I might just carefully draw it in using a mouse/tablet and ginput.

If you can live with this simplistic solution, then it may be a good choice for you.

[Jake]

Perfect! I will use this approach, together with what Bjorn suggested in his comment. (That will be required to get the coordinates needed for the polyshape)

Thank you so much for the detailed explaination. I'm really grateful! :)