contour and contourf produce conflicing results for certain projections

[johnhcc]

Originally posted here because I wasn't convinced it was a bug (I'm new to matplotlib)... but I'm increasingly convinced it is a bug, so here we are.

I'm plotting some NCEP temperature anomaly data using the stereographic projection. When I plot the data with filled contours, contourf, I see the first picture below. Note the red (positive values) over most of North America. When I change contourf to contour, making no other changes, I see the second picture below. Note the streak of blue (negative values) over central North America. The contour result with the blue is what is actually in the data, and therefore what I expect to see. The contourf result is not what I expect. Is this a bug, or am I using these tools improperly? (Update: The first plot above should look like the below plot, which was produced using the cylindrical projection:) The python code is reproduced below... you can obtain the data as a small binary file here if you want to try it yourself.

Some more notes... (1) shifting the longitude range to -180 to 180 using addcyclic and shiftgrid does not help. (2) the Lambert Conformal ('lcc') projection produces results that are less wrong, but still not correct, when using contourf. (3) the Alberts Equal Area ('aea') projection centered on North America looks correct, however when I use the North Pole-centric version ('nplaea') I only see correct results for certain values of lon_0. I am happy to provide an example if this turns out to be a different issue.

import numpy as np
import matplotlib.pyplot as plt

# read data

f = open('data5.bin', 'r')
lat = np.fromfile(f,dtype=np.float32,count=73)
lon = np.fromfile(f,dtype=np.float32,count=144)
data = np.reshape(np.fromfile(f,dtype=np.float32,count=-1),(73,144))
f.close()

# plot

m = Basemap(width=10000000,height=6000000,
        resolution='l',projection='stere',\
        lat_ts=50,lat_0=50,lon_0=253)
m.drawcoastlines()
lon2d, lat2d = np.meshgrid(lon,lat)
x, y = m(lon2d,lat2d)
mymap = plt.contourf(x,y,data,levels=np.arange(17)-8,cmap=plt.cm.bwr)
plt.colorbar(mymap,orientation='vertical',shrink=0.75)
plt.show()

[efiring]

What version of matplotlib are you using?

[johnhcc]

Am using version 1.5.0 of matplotlib.

[efiring]

In your call to contourf, try adding the kwarg corner_mask='legacy'. This will change to a different contouring algorithm.

[johnhcc]

I added corner_mask='legacy' to my contourf call and see no difference in the plot, just the warning about the deprecated attribute.

[guziy]

Is there any chance that those negative values are very small??

maybe z[abs(z) < 1e-2] = 0

could help?

Cheers

2015-11-23 12:43 GMT-05:00 jhcsu notifications@github.com:

I added corner_mask='legacy' and see no difference in the plot, just the warning about the deprecated attribute.

— Reply to this email directly or view it on GitHub https://github.com/matplotlib/basemap/issues/229#issuecomment-159007842.

Sasha

[johnhcc]

The negative values shown in the blue contours are not particularly small. The value at 50 degrees north, 255 degrees east (southern Saskatchewan), for example, is approximately -2.4:

print (lat[16], lon[102], data[16,102])
(50.0, 255.0, -2.4153743)

[johnhcc]

I've updated my original report to include the following... this was created using the cylindrical projection, and seems to be correct. Note that it looks much like the contour version of the stereographic plot (but not much like the contourf version).

[efiring]

On 2015/11/23 8:29 AM, jhcsu wrote:

I've updated my original report to include the following... this was created using the cylindrical projection, and seems to be correct. Note that it looks much like the contour version (but not much like the contourf version).

Try one more thing (another shot in the dark): use floating point for your levels.

[johnhcc]

Try one more thing (another shot in the dark): use floating point for your levels.

I changed the code to levels=np.arange(17)-8., forcing conversion to floating point, and no difference. I can appreciate shots in the dark, though!

[johnhcc]

This is odd. If I subset longitude starting at lon[30]=75.0 (which isn't even on the map!), it works. If I subset longitude at lon[29]=72.5 (or don't subset at all, as in my initial example), it doesn't. What is special about 75 degrees longitude?

Let's demo by making a slight change just to this segment of the code. The following works:

lon2d, lat2d = np.meshgrid(lon[30:],lat)
x, y = m(lon2d,lat2d)
mymap = plt.contourf(x,y,data[:,30:],levels=np.arange(17)-8,cmap=plt.cm.bwr)

The following does not (and looks exactly like my initial example):

lon2d, lat2d = np.meshgrid(lon[29:],lat)
x, y = m(lon2d,lat2d)
mymap = plt.contourf(x,y,data[:,29:],levels=np.arange(17)-8,cmap=plt.cm.bwr)

[tacaswell]

Do the collection objects created by the contourf calls have enough layers (and the rendering is failing) or is contourf is failing to create the correct collections?

[johnhcc]

Do the collection objects created by the contourf calls have enough layers (and the rendering is failing) or is contourf is failing to create the correct collections?

Although I can't answer that question... I will clarify that no error is produced. If I didn't already have some expectation of what the output looked like, I would have incorrectly assumed it was correct.