What controls curvature in contour, contourm, contourfm,... ?

Hi all! I've been digging through the documentation and google all day but I can't figure out what controls the curvature of contourlines in MATLAB. I am mostly using contourm (mapping toolbox). I have produced a plot

and now I'd like to know, why did MATLAB connect all the areas of darkest red (instead of giving me multiple blobs)? It seems to me that the reason might be a maximum allowed curvature?

Any ideas? Thanks!!

[Chad Greene]

What data did you give contourfm to create the plot above?

[LMS]

vectors x,y,z where x is latitude, y is longitude, and z is elevation of the black points.

Answer 1

[John D'Errico]

I was hoping someone who had the mapping toolbox would answer this, so I could avoid stepping in. But the answer is clear enough, and nobody else has. :)

Contouring can be different, depending on the specific case. So the code for contour (which presumes a gridded array of point on which to build the contour) can differ from that of a tool like contourm, which apparently can handle an unstructured set of scattered points.

Given an unstructured set, the only viable approach for such automatic contouring is a simple one. (Yes, one could use more complex schemes, but they will tend to fail far more easily. So for AUTOMATIC, ROBUST contouring from scattered data, I'll argue there is only one simple approach that I would choose.) In fact, the approach for contouring scattered data is surely quite closely related to how one would contour a regular lattice, but for one crucial step.

So in order to explain it, I'll first need to explain how one does contouring on a regular lattice. (This is why I tried to avoid stepping in initially, since this post will be lengthy. And most of the time on a lengthy response like this, the OP will step in with a comment like "I know all that. But why...")

Given a regular lattice, so an n by m array of points, where the spacing along each axis may be uniform or not, such as can be used by contour, the code would typically employ a simple approach. In fact, that scheme is easily generalized into higher dimensions for iso-surface tools.

1. Break each rectangle in the lattice into a pair of triangles, using one arbitrarily chosen orientation for the diagonal. For simplicity, one would choose the same orientation for that dissection for all rectangles in the lattice, although that is not technically required. These dissections create a simple triangulation of the lattice, which is easy to build, although that triangulation is never explicitly constructed, using a call to tools like delaunay.

2. Test every triangle s created. If the contour crosses any edge of the triangle, then it MUST cross exactly two triangles edges. (Here, I'm going to ignore degenerate cases.)

3. Find those two points where the contour crosses those two indicated edges, using simple linear interpolation along the edges. Two points determine a line.

4. If we assume a linear interpolant across the triangle, then the contour would be a straight line in that local area, and we have just identified the two points where the contour crosses that triangle. So again, two points determine a line, so we have now generated ONE line segment of a contour line.

5. We can now repeat the above scheme for every triangle, testing to see if the contour crossed it, or we can employ a marching scheme, wherein we might follow the contour from triangle to triangle. Either method works, as long as one is careful in the code.

Ok, so that is a rough outline of how a general contouring scheme would work. I call it rough because I am not now looking at the specific code used in a tool like contour, but because I have written complete tools for this purpose myself. A contouring scheme that works on a regular lattice CAN (if the person writing the code wants to do so) get tricky in the dissection process for each square, deciding which orientation one might employ to dissect any small rectangle. That will resolve cases where the contour could be broken down into smaller pieces. Anyway, filling in all the details takes some work, because you need to worry about degeneracies, about all the things to make it robust. Since the basic MATAB contouring tools return a nicely structured set of contours, they probably employ some sort of marching triangles scheme, but that is not necessary, and is irrelevant to the case at hand.

A contouring scheme that works on scattered data, generating a plot like the one shown, will employ a similar scheme, but for one crucial step in the beginning. Here we start out using a Delaunay triangulation of the points in the (x,y) plane. Given that triangulation, is it trivial to fall into the general outline I describe above to generate contours.

To convince yourself of that fact, look at the coloring employed. See that it shows what are essentially large triangles in regions where the data is sparse. Also, see that the domain of the contour is a convex one, even though the data has holes in it along the edges. There are triangles that completely cover the convex hull of the data. A Delaunay triangulation was surely employed here.

Given that Delaunay triangulation, using linear interpolation across each triangle is surely employed, again, one can tell this simply from looking at the coloring, but also from understanding how such a code would be written.

However, you ask about curvature. There is NO question about local curvature. That will not be employed or measured anywhere in the code. Just simple linear interpolants across each triangle in the tessellation. While one could get trickier, worrying about the shape of a contour line for a regular lattice of points, that will NOT be done for a scattered set of points.

[LMS]

Hi there, and thank you for your answer! My bad, I didn't understand the comment above for what went into contourm so I answered incorrectly. Let me give you my code:

plot3m(y,x,z,'.','markersize',12,'color','black')
hold on
% set up mesh grid for x/y coords, step-size=0.5
[yi,xi] = meshgrid(min(y):0.5:max(y), min(x):0.5:max(x));
% interpolate z-values:
zi = griddata(y,x,z,yi,xi);
% plot the interpolated surface:
contourfm(yi,xi,zi) %hard contour lines, or alternatively:
%surfm(yi,xi,zi) %smooth surface with gradients

so it is in fact structured data I'm feeding into it, my apologies!

Thank you so much for your extensive explanation, I did in fact not know any of that, and I'm grateful that you explained it. So if I got you right, I cannot control the contour curvature in any way?

[John D'Errico]

I don't have MATLAB on this machine, and even then, I don't have the mapping toolbox, so I cannot be specific about contourfm or tools like surfm.

Yes. You generated a nicely structured grid when you called griddata there. griddata used a delaunay triangulation of your data. Outside of the convex hull, it will have produced NaNs, so nothing is plotted there.

But when you generated the interpolated result using griddata, you explicitly controlled how the interpolation was done, and therefore controlled the shape of those contours. The resulting surface will be piecewise locally linear, because that was what you created.

I looked briefly at the online docs for contourm. There is no control there over how the contours are generated. Again, they will be created using a locally linear interpolant, merely sets of line segments, and depending on how actively the code was written, there MAY be some effort done internally to resolve ambiguities, but that is the limit of what would be done. And when you feed it locally linearly interpolated data, there will be no effort made to improve things!

So IF you pass in locally linear interpolated data to a tool, how can it do more with the data than what you have already done? Answer: it cannot.

One option could be to try to use a tool like griddata to create a smoother interpolant. For example, you might have employed the 'natural' or 'cubic' methods for griddata. These will create what are in theory at least smoother interpolants, but ONLY WHEN THE DATA IS ITSELF PERFECTLY SMOOTH, taken from a nice, well-behaved function. So if you have noisy data, then the cubic interpolant result from gridfit will be a mess, creating all sorts of artifacts. Don't expect a well behaved result from noisy data, since when you interpolate noisy data, you also interpolate the noise. Note that when one interpolates noise, interpolants that claim to be "smooth", such as cubic splines, actually tend to amplify any noise.

So griddata is NOT a tool to use for smoothing noisy data!!!!!!!!! If your data has noise in z, that you somehow hope to smooth out, I think you will be unhappy with the results from those other methods in griddata.

You could have used a tool like my own gridfit to build that structured grid, but gridfit will extrapolate! You don't want to use it here, since those extrapolated values will be meaningless for your problem. (Yes, you could then restrict the final results from gridfit after the fact to lie entirely inside the convex hull.)

In the end, in order to gain any kind of control over the shape of the result, you will need to do more than you have done. Control is always possible, but it will take some work AND study AND code on your part to gain that control. All of that would be specific to your data and the noise structure in your data.

[LMS]

thank you. I don't expect my data be noisy, I don't think we need to worry about that. you're right though, the 'cubic' griddata resulted in a smoother plot.

that's only partially what I was looking for, mainly I would like to find out why the dots in the North and those in the South are connected, and if I can do anything to not have MATLAB connect clusters. Intuitively, I'd say that I need to limit curvature, hence the title of my original question.

To illustrate what I'm looking for, have this very crude ms paint drawing...

[John D'Errico]

See? These questions never do stop at one answer, mainly because the person asking the question usually had some underlying problem that was actually different from what they originally posed. Oh well, as long as the question stays interesting ... so I'll take another shot at this before I go off to play bridge. :)

Noise can be pertinent, even a tiny am when you have data points very near each other. A tiny amount of noise will force an interpolant to generate large gradients in regions with points very close. This is necessary in order to interpolate that data. I do see in your plot one region where there appears to be very high gradient, that may or may not be appropriate right there. It MAY be an artifact, introduced by the interpolation method.

Ignoring that issue, lets talk about your question

"why the dots in the North and those in the South are connected".

They are connected because in order to do ANY contouring, any surface modeling, you need to use some sort of interpolation or at least an implicit surface fitting scheme. Were you to not have used griddata there, but used a code that is able to work with scattered data to build contours, etc., it would still have used some scheme as I described. Where you have large holes in the data, a characteristic of a Delaunay triangulation is it will cover that hole with a large triangle. No data? No problem. Just use a large triangle to span the hole. Inside that triangle, the tool must then interpolate across the triangle, often linearly, but using a higher order scheme where you direct it to do so. In any event, the scheme will create a prediction at all points inside the convex hull of the data. This is why the dots in the North and South are connected. You could not build that contour map otherwise.

It sounds as if you really are not interested in having a global contour map as you are creating, but seem to want something that is purely local. Sadly, there are no tools that will do exactly as you desire. (At least, none that I know of.)

Were one to write something for this purpose, I might try an ad hoc method like:

1. Use a clustering tool to identify those clusters of points.

2. Build a low order LOCAL polynomial model, valid only in the immediate vicinity of each cluster so identified.

3. Build a contour map for each of those low order models, valid only in a small region, perhaps a circular disk of a fixed, known radius around that cluster.

All of this is doable. One could use a tool for explicit clustering, perhaps from the stats toolbox. Or, you can actually use a tool that builds and works with alpha shapes to do the clustering. You could use a tool like my polyfitn (on the file exchange) to build a low order polynomial model for a local surface from any cluster. Then just evaluate the polynomial over a region to then plot that local surface. Doable? Yes, very much so, with some moderate effort.

Or maybe I've misunderstood your goal here. Off to play bridge now. I'll look back in this evening.

[Chad Greene]

John, I'm glad no one else chimed in on this because I sure learned a lot from your responses. Thanks for that.

LMS--this is sort of a side note, but it's related to your question: Are you using xi and yi as lon and lat? Be careful when you do things like 2D gridding, contouring, interpolating, etc. with lats and lons because 1 degree of latitude does not have the same length as 1 degree of longiutde, yet interpolation schemes usually assume equality in both dimensions. I recommend projecting your lats and lons with projwd or mfwdtran, then grid your data.

[John D'Errico]

I had a few more thoughts to add. You are thinking about limiting curvature. This is the wrong way to look at it, as any surface that is locally fairly flat in some regions (but at differing levels, so effectively piecewise constant) is actually highly curved where it transitions between those various levels.

So you need to focus, as I suggested, on how to segregate the various populations. The classical way to do so is to use a clustering tool, if this must be done automatically. An alternative might be something like an alpha shape, although I think such a tool must fail. An idea from it might work though. The idea is if you can slide a circle of diameter alpha entirely between two clusters of points, then those clusters form a distinct population. Here you as the user must choose the value of alpha. Think of alpha as the minimum distance between clusters, so in terms of kilometers. This computation is easily done merely from an inter-point distance matrix. A nice thing about the general approach is it would be efficiently done, although as Chad points out, there could be issues with the fact that your points lie on the surface of a sphere (the Earth), so on a curved surface. The other virtue of such an approach is it does not rely on knowing in advance how many clusters of points one must expect. Many clustering tools need this information to run.

Regardless of how it is done, I would suggest that you need to extract these clusters of points as distinct populations, then model the local behavior within each cluster. You will now be able to populate the global plot from those local models.

[LMS]

Thank you, both of you, for your answers and advice. I have enough information now to proceed, and I've learned a lot. I appreciate your input and time!

All the best!

[LMS]

I got caught up in some other work, but I've now had the time to really digest everything you guys said in reply to my question. Once again thank you both for your patience and for explaining it so well.

Chad, I wanted to ask you again about your suggestion to project the coordinates to map coordinates.

If I convert my data to map coordinates, then grid the data, then interpolate - how do I plot? The mapping toolbox requires geographic coordinates to plot!

[Chad Greene]

Hey LMS,

You have two options. In some cases, the simplest option is to use regular non-mapping Matlab functions plot, pcolor, contour, etc because if you've projected all your data to (x,y) coordinates you can simply plot the data on (x,y) axes. Plotting data in (x,y) axes is common in the geosciences and when they do this they often name the axes "eastings (meters)" and "northings (meters)" and specify which projection they're using. Careful though, because eastings and northings don't always perfectly align with true east and west.

Your second option is to transform back to geo coordinates after you've done all your analysis in projected coordinates. If you've projected forward (from geo to xy) with projfwd you can get your analyzed data back into geo coordinates with projinv.

[Chad Greene]

Instead of projfwd and projinv, you may also consider these.

[LMS]

thank you so much! :)