basemap
basemap copied to clipboard
matplotlib
contour not closed across 0 meridian in polar projection
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
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 [email protected] 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(yynp.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 https://github.com/matplotlib/basemap/issues/416, or mute the thread https://github.com/notifications/unsubscribe-auth/AARy-M5dLzT2fywQ8oV9b7mpbIdpyPYgks5uMWQ_gaJpZM4VqDLA .
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?
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 [email protected] 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 https://github.com/matplotlib/basemap/issues/416#issuecomment-409729095, or mute the thread https://github.com/notifications/unsubscribe-auth/AARy-B9FliZa6vA3Dep__r6Rckbc5Z29ks5uMhzigaJpZM4VqDLA .
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'
