Problem doing plots using contour with NaN elements in a matrix

Dear all,

I have 3 matrices: lon, lat, depth.

Than I want to remove all the depth values that are bigger than 0. To do so I use:

depth(depth>1)=NaN.

By doing this I remove the land part, which is ok.

The problem then is that I want to create a contour arround the sea level, which means at the depth =0. To do so I'm doing:

[C,H]=contour(lon,lat,depth',[0 0]);

As a result I obtain an empty C matrix and I don't understand why!!!

Hope you can help me!

Thanks in advance.

[dpb]

That would appear to indicate your resulting data don't have any zero values...not sure anybody here could do anything w/o the Z data.

[Ricardo Duarte]

The thing is:

If I do

[C,H]=contour(lon,lat,depth', [0 0]) 

before the depth(depth>1)=NaN it gives me a very nice contour plot. However, if I do it after the C matrix is empty.

I attached my matrices.

Thank you in advance.

[dpb]

I think that clearly shows you that there's nothing left for contour to have interpolated over to find the 0 level contour with NaN on the boundariess. But, only looking at the data itself will explain...I'll try to take a look later on.

[dpb]

So, what is it you are trying to do here, anyways? Why are you setting to NaN, anyway?

What's wrong with the returned contour for level==0 from the original array?

[Ricardo Duarte]

I want to remove land from my data in order to perform further calculations where the land should not be present.

What I want is to obtain a [C H], where the C is a matrix that has the information about the area selected and that will be used to perform later a delauny triangulation.

[dpb]

Well, that's still the same contour.

One can't interpolate across a NaN.

[Ricardo Duarte]

Yes I understand that, but supposedly the NaN are assigned to levels highers than 1 and not 0.

[dpb]

Yes, but there isn't necessarily a crossing with a value between 0 and 1 everywhere for contour to be able to find a closed path.

Again, you want the best estimate of what the zero-level contour is; you don't need to remove data to do that; in fact, you've demonstrated that doing so prevents you from finding that contour level.

[Ricardo Duarte]

The problem is that if I don't remove it, when I perform the delauny triangulation it will take into account also the land and I don't want that. I just want triangles on the water and not on land.

[dpb]

Find the contour first, THEN triangulate. If you have to remove the positions that are >0 then, so be it, but still the countour locations will have been determined and returned in the countour object.

Again, the countour is a boundary; you've got to have something that can be interpolated across the zero plane for countour() to be able to figure out where that is. When you have something that is like [-1.234 NaN] in one direction; what's it supposed to do? There simply is no way to infer where zero is between those two values.

As is, you're not removing all land, anyways, by using the >1 level; you probably will be better served to leave completely alone and find the zero contour object, then remove everything above 0.

I'm not sure what you're doing with the contour object, anyhows???

I've never messed with the guts of contour much other than recollect digging into the convoluted storage pattern of the returned contour object array for somebody here once upon a time to help them retrieve information from it, but having done is all that I do remember.

[Ricardo Duarte]

Ok, thank you very much for your help and time invested in this problem. :)

[dpb]

What are you trying to triangulate over? How do you try to combine the returned contour path with the remainder of the image that doesn't have anything excepting the elevations in it?

It's just not clear what your next step you're trying to do is -- I'm sure there's a way forward, just need more input on what your end objective really is here.

Doesn't help I don't have any of the mapping toolsets, sorry...

[Ricardo Duarte]

Ok, I will try to explain the complete problem.

I have the lon, lat, depth. With that I want to calculate the transmission loss (just at the sea) of a certain region represented by my lon, lat, depth matrices. To do so I need to define the sea level (depth=0) and remove land in order obtain my X and Y vectors and after performe the delauny triangulation. The problem here is that ,if I do not remove land before the contour, when I obtain the C matrix returned by contour, it will triangulate also in the land part and later (in the remaining part of my code) the TL will be calculated also in the land part, which is not correct.

So the code I have is:

[C,H]=contour(lon,lat,depth', [0 0]); I will use C to assing X and Y (at this stage I want to have just sea information)

X=C(1,:); Y=C(2,:);

Then I use the delauny doing this:

tri=delaunay(X,Y); (if I have land, the triangulation will be performed also in the land, which is what I don't want)

and then triplot(tri,X,Y);

And then the rest of the program...

[dpb]

You have to determine from the contour what is inside and what is outside the sea/land boundary. Removing values >0 from the depth matrix doesn't affect that; in fact by removing those you'll just make estimates of where those are less precise than will be otherwise.

I don't think the basic MATLAB delaunay function knows about how to handle holes automagically and I've not used it enough to know just how one would go about representing the connection points to do this. The routines are for convex polygons.

There is one example of setting constraints for a nonconvex polygon, but even it doesn't actually have a hole in the middle of it; it is roughly a J-shaped polygon.

All in all, I'm not so sure MATLAB is equipped to handle the problem...then again, I didn't DAGS for other examples, either.

[dpb]

A search uncovered a response from Bruno that pointed another user to <mesh2d-delaunay-based-unstructured-mesh-generation> on File Exchange

[Ricardo Duarte]

Thank you!