contourm fails to identify closed contour lines in some cases

What I have is a 2-D geolocated data grid F for hundreds of timesteps. Though it is a global dataset, I am interested in the coordinates of a certain contour line (value=2), which encircles the main parts of the north polar region.

To get these coordinates I have to make a contourm plot, there is no low-level function in the mapping toolbox (as contourc in the x/y-plane). I get a contour matrix C which contains the lat/lon coordinates of every closed contour line (or not closed if on the boundary).

The data of each contourline are separated in C, so they can be processed individually.

The code looks like that:

 lat=90:0.75:-90;
 lon=0:0.75:359.25;
 figure 
 axesm('MapProjection','stereo','Origin',[90 0 0],'MapLatLimit',[27 90])
 C = contourm(lat,lon,F,[2 2]);

In certain cases the result is as expected, shown here. Note that the contour encircles the northpole.

But in other cases (here 12 hours later) a normally closed contour line is splitted in two separated sequences in C, for illustration I have plotted these two parts.

As it seems this happens if the contour line do not encircle the pole. In these cases the first point of the green line starts with lon=0 and as it seems the algorithm to identify a single contour stops when a point reaches lon=0 again, not taking into account the points that follow. However it should be one contour line.

As I said I am mainly interested in the coordinates not in the plot. Since I am dealing with many fields I cannot effort to check in every case whether the contour is splitted and if so which sequences in C I need to join.

If a developer of the mapping toolbox read this: Is this a known bug?

Or is there a workaround for this problem?

I am grateful for any hints or suggestions.

Alex

Answer 1

[Chad Greene]

I've had this problem with contourm on stereographic projections. After you've initialized the map, try this:

[X,Y] = mfwdtran(lat,lon); 
[C,h] = contour(X,Y,F,[2 2]);

[Alexander]

Thanks Chad. Unfortunately this doesn't work for me.

I did as you said, creating the Lat/Lon matrix, then after initializing the axes projected this to map coordinates. After this I plotted the contour.

 [Lon,Lat]= meshgrid(0 :0.75:360 , 90:-0.75:-90);
 figure
 axesm( .... )
 [X,Y] = mfwdtran(Lat,Lon);
 [C,h] = contour(X,Y,F,[2 2]);

But in the resulting contour matrix C the contour line is splitted again in two parts.

EDIT: I also tried this for my first case (closed line). And there in C the contour line is also closed as for contourm.

You said you had a similar problem. Do you happen to remember if your contour lines always encircled the pole?

[Chad Greene]

Oh interesting. I'm not sure what the problem is. Perhaps an issue that Matlab doesn't know 360 degrees longitude is the same thing as 0 degrees? Most of my polar datasets are equally spaced in polar stereographic coordinates, not equally spaced in geo coordinates, so that may be the issue.

A clunky fix, but perhaps you can convert C to x,y,z coordinates with C2xyz and color the contour lines manually by plotting with plot?

[Alexander]

Meanwhile I tried to find the code for the algorithm of finding the contour lines in contourm. I didn't manage but I think no matter if the data are projected or not everything is based on the contourc function. A short description of the Contouring Algorithm can be found in the Matlab Help.

So I think you're right. It seems that Matlab cannot deal with Longitudes in a right manner if one has a global data set.

The fact that it doesn't work neither with contourm nor contour (after re-mapping on x/y) seems to confirm this.

 A clunky fix, but perhaps you can convert C to x,y,z coordinates with C2xyz and color the contour lines manually by plotting with plot?

The problem is I need the coordinates of the line rather than the plots.

But thanks anyway.