matplotlib-checkins Mailing List for matplotlib

Revision: 8333
 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=8333&view=rev
Author: jswhit
Date: 2010年05月23日 16:49:22 +0000 (2010年5月23日)
Log Message:
-----------
handle masked arrays for unstructured meshes.
Modified Paths:
--------------
 trunk/toolkits/basemap/examples/ploticos.py
 trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
Modified: trunk/toolkits/basemap/examples/ploticos.py
===================================================================
--- trunk/toolkits/basemap/examples/ploticos.py	2010年05月23日 12:31:00 UTC (rev 8332)
+++ trunk/toolkits/basemap/examples/ploticos.py	2010年05月23日 16:49:22 UTC (rev 8333)
@@ -1,6 +1,7 @@
 from mpl_toolkits.basemap import Basemap, NetCDFFile
 import matplotlib.pyplot as plt
 import numpy as np
+from numpy import ma
 # read in orography of icosahedral global grid.
 f = NetCDFFile('C02562.orog.nc')
 lons = (180./np.pi)*f.variables['grid_center_lon'][:]
@@ -11,5 +12,8 @@
 map.drawcoastlines()
 map.drawmapboundary()
 # tri=True forwards to axes.tripcolor
-map.pcolor(x,y,z,tri=True,shading='faceted')
+#z = ma.masked_where(z < 1.e-5, z) # for testing masked arrays.
+map.pcolor(x,y,z,tri=True,shading='faceted',vmin=0,vmax=3000)
+#map.contourf(x,y,z,np.arange(0,3000,150),tri=True)
+#map.contour(x,y,z,np.arange(0,3000,150),tri=True)
 plt.show()
Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
===================================================================
--- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py	2010年05月23日 12:31:00 UTC (rev 8332)
+++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py	2010年05月23日 16:49:22 UTC (rev 8333)
@@ -2757,11 +2757,24 @@
 if tri:
 # for unstructured grids, toss out points outside
 # projection limb (don't use those points in triangulation).
+ if hasattr(data,'mask'):
+ data = data.filled(fill_value=1.e30)
+ masked=True
+ else:
+ masked=False
 mask = np.logical_or(x<1.e20,y<1.e20)
 x = np.compress(mask,x)
 y = np.compress(mask,y)
 data = np.compress(mask,data)
- ret = ax.tripcolor(x,y,data,**kwargs)
+ if masked:
+ import matplotlib.tri as tri
+ triang = tri.Triangulation(x, y)
+ z = data[triang.triangles]
+ mask = (z > 1.e20).sum(axis=-1)
+ triang.set_mask(mask)
+ ret = ax.tripcolor(triang,data,**kwargs)
+ else:
+ ret = ax.tripcolor(x,y,data,**kwargs)
 else:
 # make x,y masked arrays
 # (masked where data is outside of projection limb)
@@ -2828,11 +2841,26 @@
 ax.hold(h)
 try:
 if kwargs.has_key('tri') and kwargs['tri']:
+ # for unstructured grids, toss out points outside
+ # projection limb (don't use those points in triangulation).
+ if hasattr(data,'mask'):
+ data = data.filled(fill_value=1.e30)
+ masked=True
+ else:
+ masked=False
 mask = np.logical_or(x<1.e20,y<1.e20)
 x = np.compress(mask,x)
 y = np.compress(mask,y)
 data = np.compress(mask,data)
- CS = ax.tricontour(x,y,data,*args,**kwargs)
+ if masked:
+ import matplotlib.tri as tri
+ triang = tri.Triangulation(x, y)
+ z = data[triang.triangles]
+ mask = (z > 1.e20).sum(axis=-1)
+ triang.set_mask(mask)
+ CS = ax.tricontour(triang,data,*args,**kwargs)
+ else:
+ CS = ax.tricontour(x,y,data,*args,**kwargs)
 else:
 # make sure x is monotonically increasing - if not,
 # print warning suggesting that the data be shifted in longitude
@@ -2894,11 +2922,26 @@
 ax.hold(h)
 try:
 if kwargs.has_key('tri') and kwargs['tri']:
+ # for unstructured grids, toss out points outside
+ # projection limb (don't use those points in triangulation).
+ if hasattr(data,'mask'):
+ data = data.filled(fill_value=1.e30)
+ masked=True
+ else:
+ masked=False
 mask = np.logical_or(x<1.e20,y<1.e20)
 x = np.compress(mask,x)
 y = np.compress(mask,y)
 data = np.compress(mask,data)
- CS = ax.tricontourf(x,y,data,*args,**kwargs)
+ if masked:
+ import matplotlib.tri as tri
+ triang = tri.Triangulation(x, y)
+ z = data[triang.triangles]
+ mask = (z > 1.e20).sum(axis=-1)
+ triang.set_mask(mask)
+ CS = ax.tricontourf(triang,data,*args,**kwargs)
+ else:
+ CS = ax.tricontourf(x,y,data,*args,**kwargs)
 else:
 # make sure x is monotonically increasing - if not,
 # print warning suggesting that the data be shifted in longitude
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.

Revision: 8332
 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=8332&view=rev
Author: jswhit
Date: 2010年05月23日 12:31:00 +0000 (2010年5月23日)
Log Message:
-----------
add support for unstructured meshes in pcolor, contour, contourf methods.
Modified Paths:
--------------
 trunk/toolkits/basemap/Changelog
 trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
Added Paths:
-----------
 trunk/toolkits/basemap/examples/C02562.orog.nc
 trunk/toolkits/basemap/examples/ploticos.py
Modified: trunk/toolkits/basemap/Changelog
===================================================================
--- trunk/toolkits/basemap/Changelog	2010年05月21日 18:26:12 UTC (rev 8331)
+++ trunk/toolkits/basemap/Changelog	2010年05月23日 12:31:00 UTC (rev 8332)
@@ -1,7 +1,11 @@
 version 0.99.5 (not yet released)
+ * add support for plotting on unstructured grids using
+ keyword 'tri' in pcolor,contour,contourf methods (which
+ then forward to tripcolor, tricontour, tricontourf axes
+ methods).
 	 * let continents that fill the whole map be filled.
 * added option for cubic spline interpolation in interp function
- (order=3) using scipy.ndimage.
+ (order=3) using scipy.ndimage.
 * added "near-sided perspective" projection for a satellite
 view at an arbitrary altitude.
 * patch from Stephane Raynaud to pass format string to
Added: trunk/toolkits/basemap/examples/C02562.orog.nc
===================================================================
(Binary files differ)
Property changes on: trunk/toolkits/basemap/examples/C02562.orog.nc
___________________________________________________________________
Added: svn:mime-type
 + application/octet-stream
Added: trunk/toolkits/basemap/examples/ploticos.py
===================================================================
--- trunk/toolkits/basemap/examples/ploticos.py	 (rev 0)
+++ trunk/toolkits/basemap/examples/ploticos.py	2010年05月23日 12:31:00 UTC (rev 8332)
@@ -0,0 +1,15 @@
+from mpl_toolkits.basemap import Basemap, NetCDFFile
+import matplotlib.pyplot as plt
+import numpy as np
+# read in orography of icosahedral global grid.
+f = NetCDFFile('C02562.orog.nc')
+lons = (180./np.pi)*f.variables['grid_center_lon'][:]
+lats = (180./np.pi)*f.variables['grid_center_lat'][:]
+z = f.variables['zs'][:]
+map = Basemap(projection='ortho',lon_0=-105,lat_0=40)
+x,y = map(lons, lats)
+map.drawcoastlines()
+map.drawmapboundary()
+# tri=True forwards to axes.tripcolor
+map.pcolor(x,y,z,tri=True,shading='faceted')
+plt.show()
Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
===================================================================
--- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py	2010年05月21日 18:26:12 UTC (rev 8331)
+++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py	2010年05月23日 12:31:00 UTC (rev 8332)
@@ -2730,7 +2730,7 @@
 self.set_axes_limits(ax=ax)
 return ret
 
- def pcolor(self,x,y,data,**kwargs):
+ def pcolor(self,x,y,data,tri=False,**kwargs):
 """
 Make a pseudo-color plot over the map
 (see matplotlib.pyplot.pcolor documentation).
@@ -2739,22 +2739,35 @@
 they will be convert to masked arrays with those values masked.
 As a result, those values will not be plotted.
 
+ If ``tri`` is set to ``True``, an unstructured grid is assumed
+ (x,y,data must be 1-d) and matplotlib.pyplot.tricolor is used.
+
 Extra keyword ``ax`` can be used to override the default axis instance.
 
- Other \**kwargs passed on to matplotlib.pyplot.pcolor.
+ Other \**kwargs passed on to matplotlib.pyplot.pcolor (or tricolor if
+ ``tri=True``).
 """
 ax, plt = self._ax_plt_from_kw(kwargs)
- # make x,y masked arrays
- # (masked where data is outside of projection limb)
- x = ma.masked_values(np.where(x > 1.e20,1.e20,x), 1.e20)
- y = ma.masked_values(np.where(y > 1.e20,1.e20,y), 1.e20)
 # allow callers to override the hold state by passing hold=True|False
 b = ax.ishold()
 h = kwargs.pop('hold',None)
 if h is not None:
 ax.hold(h)
 try:
- ret = ax.pcolor(x,y,data,**kwargs)
+ if tri:
+ # for unstructured grids, toss out points outside
+ # projection limb (don't use those points in triangulation).
+ mask = np.logical_or(x<1.e20,y<1.e20)
+ x = np.compress(mask,x)
+ y = np.compress(mask,y)
+ data = np.compress(mask,data)
+ ret = ax.tripcolor(x,y,data,**kwargs)
+ else:
+ # make x,y masked arrays
+ # (masked where data is outside of projection limb)
+ x = ma.masked_values(np.where(x > 1.e20,1.e20,x), 1.e20)
+ y = ma.masked_values(np.where(y > 1.e20,1.e20,y), 1.e20)
+ ret = ax.pcolor(x,y,data,**kwargs)
 except:
 ax.hold(b)
 raise
@@ -2801,40 +2814,51 @@
 
 Extra keyword ``ax`` can be used to override the default axis instance.
 
- Other \*args and \**kwargs passed on to matplotlib.pyplot.contour.
+ If ``tri`` is set to ``True``, an unstructured grid is assumed
+ (x,y,data must be 1-d) and matplotlib.pyplot.tricontour is used.
+
+ Other \*args and \**kwargs passed on to matplotlib.pyplot.contour
+ (or tricontour if ``tri=True``).
 """
 ax, plt = self._ax_plt_from_kw(kwargs)
- # make sure x is monotonically increasing - if not,
- # print warning suggesting that the data be shifted in longitude
- # with the shiftgrid function.
- # only do this check for global projections.
- if self.projection in _cylproj + _pseudocyl:
- xx = x[x.shape[0]/2,:]
- condition = (xx >= self.xmin) & (xx <= self.xmax)
- xl = xx.compress(condition).tolist()
- xs = xl[:]
- xs.sort()
- if xl != xs:
- print dedent("""
- WARNING: x coordinate not montonically increasing - contour plot
- may not be what you expect. If it looks odd, your can either
- adjust the map projection region to be consistent with your data, or
- (if your data is on a global lat/lon grid) use the shiftgrid
- function to adjust the data to be consistent with the map projection
- region (see examples/contour_demo.py).""")
- # mask for points outside projection limb.
- xymask = np.logical_or(np.greater(x,1.e20),np.greater(y,1.e20))
- data = ma.asarray(data)
- # combine with data mask.
- mask = np.logical_or(ma.getmaskarray(data),xymask)
- data = ma.masked_array(data,mask=mask)
 # allow callers to override the hold state by passing hold=True|False
 b = ax.ishold()
 h = kwargs.pop('hold',None)
 if h is not None:
 ax.hold(h)
 try:
- CS = ax.contour(x,y,data,*args,**kwargs)
+ if kwargs.has_key('tri') and kwargs['tri']:
+ mask = np.logical_or(x<1.e20,y<1.e20)
+ x = np.compress(mask,x)
+ y = np.compress(mask,y)
+ data = np.compress(mask,data)
+ CS = ax.tricontour(x,y,data,*args,**kwargs)
+ else:
+ # make sure x is monotonically increasing - if not,
+ # print warning suggesting that the data be shifted in longitude
+ # with the shiftgrid function.
+ # only do this check for global projections.
+ if self.projection in _cylproj + _pseudocyl:
+ xx = x[x.shape[0]/2,:]
+ condition = (xx >= self.xmin) & (xx <= self.xmax)
+ xl = xx.compress(condition).tolist()
+ xs = xl[:]
+ xs.sort()
+ if xl != xs:
+ print dedent("""
+ WARNING: x coordinate not montonically increasing - contour plot
+ may not be what you expect. If it looks odd, your can either
+ adjust the map projection region to be consistent with your data, or
+ (if your data is on a global lat/lon grid) use the shiftgrid
+ function to adjust the data to be consistent with the map projection
+ region (see examples/contour_demo.py).""")
+ # mask for points outside projection limb.
+ xymask = np.logical_or(np.greater(x,1.e20),np.greater(y,1.e20))
+ data = ma.asarray(data)
+ # combine with data mask.
+ mask = np.logical_or(ma.getmaskarray(data),xymask)
+ data = ma.masked_array(data,mask=mask)
+ CS = ax.contour(x,y,data,*args,**kwargs)
 except:
 ax.hold(b)
 raise
@@ -2856,48 +2880,59 @@
 
 Extra keyword 'ax' can be used to override the default axis instance.
 
- Other \*args and \**kwargs passed on to matplotlib.pyplot.scatter.
+ If ``tri`` is set to ``True``, an unstructured grid is assumed
+ (x,y,data must be 1-d) and matplotlib.pyplot.tricontourf is used.
+
+ Other \*args and \**kwargs passed on to matplotlib.pyplot.contourf
+ (or tricontourf if ``tri=True``).
 """
 ax, plt = self._ax_plt_from_kw(kwargs)
- # make sure x is monotonically increasing - if not,
- # print warning suggesting that the data be shifted in longitude
- # with the shiftgrid function.
- # only do this check for global projections.
- if self.projection in _cylproj + _pseudocyl:
- xx = x[x.shape[0]/2,:]
- condition = (xx >= self.xmin) & (xx <= self.xmax)
- xl = xx.compress(condition).tolist()
- xs = xl[:]
- xs.sort()
- if xl != xs:
- print dedent("""
- WARNING: x coordinate not montonically increasing - contour plot
- may not be what you expect. If it looks odd, your can either
- adjust the map projection region to be consistent with your data, or
- (if your data is on a global lat/lon grid) use the shiftgrid
- function to adjust the data to be consistent with the map projection
- region (see examples/contour_demo.py).""")
- # mask for points outside projection limb.
- xymask = np.logical_or(np.greater(x,1.e20),np.greater(y,1.e20))
- # mask outside projection region (workaround for contourf bug?)
- epsx = 0.1*(self.xmax-self.xmin)
- epsy = 0.1*(self.ymax-self.ymin)
- outsidemask = np.logical_or(np.logical_or(x > self.xmax+epsx,\
- x < self.xmin-epsy),\
- np.logical_or(y > self.ymax+epsy,\
- y < self.ymin-epsy))
- data = ma.asarray(data)
- # combine masks.
- mask = \
- np.logical_or(outsidemask,np.logical_or(ma.getmaskarray(data),xymask))
- data = ma.masked_array(data,mask=mask)
 # allow callers to override the hold state by passing hold=True|False
 b = ax.ishold()
 h = kwargs.pop('hold',None)
 if h is not None:
 ax.hold(h)
 try:
- CS = ax.contourf(x,y,data,*args,**kwargs)
+ if kwargs.has_key('tri') and kwargs['tri']:
+ mask = np.logical_or(x<1.e20,y<1.e20)
+ x = np.compress(mask,x)
+ y = np.compress(mask,y)
+ data = np.compress(mask,data)
+ CS = ax.tricontourf(x,y,data,*args,**kwargs)
+ else:
+ # make sure x is monotonically increasing - if not,
+ # print warning suggesting that the data be shifted in longitude
+ # with the shiftgrid function.
+ # only do this check for global projections.
+ if self.projection in _cylproj + _pseudocyl:
+ xx = x[x.shape[0]/2,:]
+ condition = (xx >= self.xmin) & (xx <= self.xmax)
+ xl = xx.compress(condition).tolist()
+ xs = xl[:]
+ xs.sort()
+ if xl != xs:
+ print dedent("""
+ WARNING: x coordinate not montonically increasing - contour plot
+ may not be what you expect. If it looks odd, your can either
+ adjust the map projection region to be consistent with your data, or
+ (if your data is on a global lat/lon grid) use the shiftgrid
+ function to adjust the data to be consistent with the map projection
+ region (see examples/contour_demo.py).""")
+ # mask for points outside projection limb.
+ xymask = np.logical_or(np.greater(x,1.e20),np.greater(y,1.e20))
+ # mask outside projection region (workaround for contourf bug?)
+ epsx = 0.1*(self.xmax-self.xmin)
+ epsy = 0.1*(self.ymax-self.ymin)
+ outsidemask = np.logical_or(np.logical_or(x > self.xmax+epsx,\
+ x < self.xmin-epsy),\
+ np.logical_or(y > self.ymax+epsy,\
+ y < self.ymin-epsy))
+ data = ma.asarray(data)
+ # combine masks.
+ mask = \
+ np.logical_or(outsidemask,np.logical_or(ma.getmaskarray(data),xymask))
+ data = ma.masked_array(data,mask=mask)
+ CS = ax.contourf(x,y,data,*args,**kwargs)
 except:
 ax.hold(b)
 raise
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.