contour not closed across 0 meridian in polar projection

[Xunius]

Hi all,
I'm trying to fetch contour line coordinates from a npaeqd projection plot and I noticed that the contour lines will be broken into 2 parts whenever it crosses the 0-degree longitude, even though they form a closed contour and after calling addcyclic(). Below is minimal working example:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import addcyclic
from mpl_toolkits.basemap import Basemap

lats=np.linspace(0,90,90)
lons=np.linspace(0,360,360)

# make some toy data
xx,yy=np.meshgrid(lons,lats)
z=np.cos(xx*np.pi/180)*np.sin(yy*np.pi/180)

# add cyclic
z,lons=addcyclic(z,lons)
xx,yy=np.meshgrid(lons,lats)

# get contours
bmap=Basemap(projection='npaeqd',boundinglat=0,lon_0=0,
        resolution='l')

fig=plt.figure(figsize=(12,6))
ax1=fig.add_subplot(1,2,1)

contours=bmap.contour(xx,yy,z,[-0.6,0.6],latlon=True,ax=ax1)
bmap.drawcoastlines(ax=ax1)

clines1=contours.collections[0].get_paths()
clines2=contours.collections[1].get_paths()

print 'len(clines1), num of contours across 180', len(clines1)
print 'len(clines2), num of contours across 0', len(clines2)

# plot contours
ax2=fig.add_subplot(1,2,2)

xs=clines1[0].vertices[:,0]
ys=clines1[0].vertices[:,1]
ax2.plot(xs,ys,'b-',label='Contour across 180')

xs=clines2[0].vertices[:,0]
ys=clines2[0].vertices[:,1]
ax2.plot(xs,ys,'r-',label='Half contour across 0')

xs=clines2[1].vertices[:,0]
ys=clines2[1].vertices[:,1]
ax2.plot(xs,ys,'g-',label='Half contour across 0')

ax2.legend()

plt.show(block=False)

Figure output here
The yellow contour on the left are made up by 2 lines (red+green) on the right. This makes it difficult when I try to detect and track some features that move across the 0-meridian.
Is it intended or a bug?

Some specs:
basemap 1.0.7
matplotlib 2.2.2, both installed via conda install

[WeatherGod]

What is getting contoured here spans all of 0 to 360. A discontinuity has to happen somewhere. In this case, it is at 0/360. Now, if you projected the data ahead of time and contoured that, then you wouldn't have a discontinuity at any longitude.

…

On Wed, Aug 1, 2018 at 4:18 AM, Xunius ***@***.***> wrote: Hi all, I'm trying to fetch contour line coordinates from a npaeqd projection plot and I noticed that the contour lines will be broken into 2 parts whenever it crosses the 0-degree longitude, even though they form a closed contour and after calling addcyclic(). Below is minimal working example: import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.basemap import addcyclic from mpl_toolkits.basemap import Basemap lats=np.linspace(0,90,90) lons=np.linspace(0,360,360) # make some toy data xx,yy=np.meshgrid(lons,lats) z=np.cos(xx*np.pi/180)*np.sin(yy*np.pi/180) # add cyclic z,lons=addcyclic(z,lons) xx,yy=np.meshgrid(lons,lats) # get contours bmap=Basemap(projection='npaeqd',boundinglat=0,lon_0=0, resolution='l') fig=plt.figure(figsize=(12,6)) ax1=fig.add_subplot(1,2,1) contours=bmap.contour(xx,yy,z,[-0.6,0.6],latlon=True,ax=ax1) bmap.drawcoastlines(ax=ax1) clines1=contours.collections[0].get_paths() clines2=contours.collections[1].get_paths() print 'len(clines1), num of contours across 180', len(clines1) print 'len(clines2), num of contours across 0', len(clines2) # plot contours ax2=fig.add_subplot(1,2,2) xs=clines1[0].vertices[:,0] ys=clines1[0].vertices[:,1] ax2.plot(xs,ys,'b-',label='Contour across 180') xs=clines2[0].vertices[:,0] ys=clines2[0].vertices[:,1] ax2.plot(xs,ys,'r-',label='Half contour across 0') xs=clines2[1].vertices[:,0] ys=clines2[1].vertices[:,1] ax2.plot(xs,ys,'g-',label='Half contour across 0') ax2.legend() plt.show(block=False) Figure output here <https://imgur.com/a/XxgcUCS> The yellow contour on the left are made up by 2 lines (red+green) on the right. This makes it difficult when I try to detect and track some features that move across the 0-meridian. Is it intended or a bug? Some specs: basemap 1.0.7 matplotlib 2.2.2, both installed via conda install — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#416>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AARy-M5dLzT2fywQ8oV9b7mpbIdpyPYgks5uMWQ_gaJpZM4VqDLA> .

[Xunius]

Cheers @WeatherGod. I assume by projecting the data you meant calling the m.transform_scalar() function? I thought basemap did that before doing the contour under the hood so the contours should be closed.
I'm currently trying to manually link those discontinuities as their coordinates do overlap. I'm bit concerned with possible accuracy lost in the projection approach. I've to project the data to do contouring, and also project the longitude/latitude mesh so I could get lon/lat coordinates from the contours, do I?

[WeatherGod]

Huh, you are right, it does project the stuff ahead of time for you, So now I am confused why there is that discontinuity...

…

On Wed, Aug 1, 2018 at 5:25 PM, Xunius ***@***.***> wrote: Cheers @WeatherGod <https://github.com/WeatherGod>. I assume by projecting the data you meant calling the m.transform_scalar() function? I thought basemap did that before doing the contour under the hood so the contours should be closed. I'm currently trying to manually link those discontinuities as their coordinates do overlap. I'm bit concerned with possible accuracy lost in the projection approach. I've to project the data to do contouring, and also project the longitude/latitude mesh so I could get lon/lat coordinates from the contours, do I? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#416 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AARy-B9FliZa6vA3Dep__r6Rckbc5Z29ks5uMhzigaJpZM4VqDLA> .

[Xunius]

Or it computes the contour and only project the contour?

[chayanroyc]

I still get this error, File "C:\Users\acer\Anaconda3\lib\site-packages\mpl_toolkits\basemap\__init__.py", line 5111, in addcyclic return list(map(_addcyclic,arr[:-1]) + [_addcyclic_lon(arr[-1])]) TypeError: unsupported operand type(s) for +: 'map' and 'list'

[Xunius]

@chayanroyc See this issue #440