Menu
▾
▴
matplotlib-devel
|
From: Curtis C. <cu...@hi...> - 2004-09-30 23:54:46
|
Hi, My research is in computational fluid dynamics (specifically, the meteorologies of planetary atmospheres). Working contour and vector plots in matplotlib would make it possible for me to make 2D meteorological maps of atmospheric layers, etc. I noticed for the first time in the goals page that contour plots are being worked on, apparently by STSci. I have been considering implementing these two plot types as sets of line collections, but now that I know contour plots are being worked on, and vector plots are simpler to implement (in 2D), I will work on making vector plots. The mathematics is fairly straightforward. I just need to learn how to use the class library. About contour plots, however, I have a couple of questions. How is it being implemented? I was about to try to write a marching squares contouring routine, although it might have been painfully slow in Python. Does anyone have experience with this? Cheers, Curtis |
|
From: Perry G. <pe...@st...> - 2004-10-01 02:17:03
|
Curtis Cooper writes: > My research is in computational fluid dynamics (specifically, the > meteorologies of planetary atmospheres). Working contour and vector plots > in matplotlib would make it possible for me to make 2D meteorological maps > of atmospheric layers, etc. > > I noticed for the first time in the goals page that contour plots are > being worked on, apparently by STSci. I have been considering > implementing these two plot types as sets of line collections, but now > that I know contour plots are being worked on, and vector plots are > simpler to implement (in 2D), I will work on making vector plots. The > mathematics is fairly straightforward. I just need to learn how to use > the class library. > > About contour plots, however, I have a couple of questions. How is it > being implemented? I was about to try to write a marching squares > contouring routine, although it might have been painfully slow in Python. > Does anyone have experience with this? > We are trying to adapt the C contour program that is used by gist (and can be found in the contour routine used by xplt in scipy). It would be best to look at the source for the precise description of the algorithm it uses (note though that gist apparently uses two different pieces of contour code for its contour tasks. The one we are looking to adapt, mainly because it appears much easier to isolate from the gist environment is the gcntr.c version). I would be amazed if one could find a pure Python algorithm to do contouring that was fast enough. Our current plan is to use these C routines to generate the contour segments, and do the plotting from within Python (as well as any contour labeling). If you have expertise in this area you may be able to do it better and faster than we can. Currently it is being worked on part time so we aren't able to do it as fast as we would like. I'm hoping that we will have at least a basic version (e.g., no labeling) in a couple weeks. If you want me to send or point you to the source code we are using as the basis, let me know. Perry Greenfield |
|
From: Curtis C. <cu...@hi...> - 2004-10-01 03:29:35
|
> We are trying to adapt the C contour program that is used by gist > (and can be found in the contour routine used by xplt in scipy). > It would be best to look at the source for the precise description > of the algorithm it uses (note though that gist apparently uses > two different pieces of contour code for its contour tasks. The > one we are looking to adapt, mainly because it appears much easier > to isolate from the gist environment is the gcntr.c version). > I would be amazed if one could find a pure Python algorithm to do > contouring that was fast enough. Our current plan is to use these > C routines to generate the contour segments, and do the plotting > from within Python (as well as any contour labeling). > > If you have expertise in this area you may be able to do it better > and faster than we can. Currently it is being worked on part time > so we aren't able to do it as fast as we would like. I'm hoping that we > will have at least a basic version (e.g., no labeling) in a couple > weeks. > > If you want me to send or point you to the source code we are > using as the basis, let me know. Hi again, Yes, I had the thought that using C for the algorithm would be easier as well. There are actually some very well-written marching squares contouring algorithms in C already out there. I will try to find such an implementation and point you to it or send you the source code. The second half is just the drawing, which should be implemented in matplotlib using the line collections class. Since vector plotting is not that hard, I will try to get that working first. Then, someone can take my source code and adapt it easily to the contouring problem, once an effective and sufficiently high-performance algorithm implementation can be found. Cheers, Curtis |
|
From: Perry G. <pe...@st...> - 2004-10-01 03:32:48
|
> Hi again, > > Yes, I had the thought that using C for the algorithm would be easier as > well. There are actually some very well-written marching squares > contouring algorithms in C already out there. I will try to find such an > implementation and point you to it or send you the source code. > Thanks, that would be helpful. In my search I didn't come across many. Keep in mind the license needs to be compatible with that of matplotlib. > The second half is just the drawing, which should be implemented in > matplotlib using the line collections class. Since vector plotting is not Yeah, that's what we have in mind. > that hard, I will try to get that working first. Then, someone can take > my source code and adapt it easily to the contouring problem, once an > effective and sufficiently high-performance algorithm implementation can > be found. > > Cheers, > Curtis > OK, Perry |
|
From: Helge A. <av...@ii...> - 2004-10-01 15:44:49
|
John Hunter <jdh...@ac...> writes: | I was concerned by the fact that the lines were not smooth - if you | plot a connected line they line jumps from side to side. But it | does get the contour right, and is implemented in pure numeric, and | so it occurs to me that it might be easier to fix this problem than | start from scratch. Perhaps Helge or one of you has some insight | into how to fix this. | | I'm attaching a modified version of the tarfile Helge initially sent | me. I've included a script testkont_mpl.py that calls Helge's lib. | Change the '.' linestyle to '-' to see the problem I discussed. Hi, not sure if I have matplotlib 100% correctly installed, but this is what I see using your example script: http://www.ii.uib.no/~avle/mpl/c0.png (and with the current algorithm, more or less what I would expect...) to get straight lines you must plot segments one by one since they are not ordered. if I use gist for this(see the script at the end) I get http://www.ii.uib.no/~avle/mpl/c1.png the first points of the segments are given by the vectors (x1,y1) the second (x2,y2). you can get pretty lines in matplotlib as well, but only by using the scattered line drawing methods of gtk. (something like self.area.window.draw_segments(self.gc, zip( x1,y1,x2,y2)?) if you want do do it "right" in matplotlib, you should implement a contour following algorithm (in C) - with this I mean an routine that returns the linesegments defining each countour in bundles. the current alg. is sort of marching cubes in 2D, a simplified version of CONREC http://astronomy.swin.edu.au/~pbourke/projection/conrec/ but only using 2 triangles per square. doing contour following alg. it is also much easier to implement automatic contour labelling. I suspect python loops are too slow for such algorithms - it may perhaps be possible to do them in Numeric, but it will still be much slower than my simple library. I think you may use the GPL'ed PLPLOT (C) for an example of contour following alg. Helge from matplotlib.matlab import * import hutil delta = 0.05 x = y = arange(-3.0, 3.0, delta) X, Y = meshgrid(x, y) Z1 = bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0) Z2 = bivariate_normal(X, Y, 1.5, 0.5, 1, 1) Z = Z2-Z1 print Z.shape #fsm = ones(Z.shape, Z.typecode()) fsm = ones(Z.shape, 'l') zmax, zmin = hutil.maxmin(Z) depths=linspace(zmin, zmax, 10) x1,y1,x2,y2 = hutil.contour2(Z, fsm, depths ) #imshow(Z, origin='lower', interpolation='nearest') #plot(y2,x2,'-') #show() import gist gist.pldefault(dpi=100,style='framed.gs') gist.palette('rainbow.gp') gist.pli(transpose(Z)) gist.pldj(x1,y1,x2,y2) # draw disjoint segments |
|
From: John H. <jdh...@ac...> - 2004-10-01 16:04:34
|
>>>>> "Helge" == Helge Avlesen <av...@ii...> writes:
Helge> http://www.ii.uib.no/~avle/mpl/c1.png
Helge> the first points of the segments are given by the vectors
Helge> (x1,y1) the second (x2,y2). you can get pretty lines in
Helge> matplotlib as well, but only by using the scattered line
Helge> drawing methods of gtk. (something like
Helge> self.area.window.draw_segments(self.gc, zip( x1,y1,x2,y2)?)
OK, I see. I didn't fully understand that x1,x2,y1,y2 were the verts
of unordered line segments. Then one can easily use a LineCollection
to draw these efficiently in matplotlib - script below and screenshot
http://nitace.bsd.uchicago.edu:8080/files/share/kontour.png. Jeez, I
feel bad for sitting on this since February!
Helge> if you want do do it "right" in matplotlib, you should
Helge> implement a contour following algorithm (in C) - with this
Helge> I mean an routine that returns the linesegments defining
Helge> each countour in bundles. the current alg. is sort of
Helge> marching cubes in 2D, a simplified version of CONREC
Helge> http://astronomy.swin.edu.au/~pbourke/projection/conrec/
Helge> but only using 2 triangles per square.
Do you have any thoughts on how we might do labels with your code?
Helge> doing contour following alg. it is also much easier to
Helge> implement automatic contour labelling. I suspect python
Helge> loops are too slow for such algorithms - it may perhaps be
Helge> possible to do them in Numeric, but it will still be much
Helge> slower than my simple library. I think you may use the
Helge> GPL'ed PLPLOT (C) for an example of contour following alg.
We have a problem in that we cannot use GPL'd code in matplotlib
because the GPL does not allow redistribution of closed code, which
the matplotlib (and python license) do.
If we decide to go with your routines, at least for the time being
until we can "do it right", would you be willing to contribute your
code to matplotlib under the matplotlib license (PSF inspired, free
for commercial and noncommercial reuse)?
Thanks!
JDH
from matplotlib.matlab import *
from matplotlib.collections import LineCollection
import hutil
delta = 0.05
x = y = arange(-3.0, 3.0, delta)
X, Y = meshgrid(x, y)
Z1 = bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
Z2 = bivariate_normal(X, Y, 1.5, 0.5, 1, 1)
Z = Z2-Z1
print Z.shape
fsm = ones(Z.shape, Z.typecode())
zmax, zmin = hutil.maxmin(Z)
depths=linspace(zmin, zmax, 10)
x1,y1,x2,y2 = hutil.contour2(Z, fsm, depths )
imshow(Z, origin='lower', interpolation='nearest')
segments = [ ( (thisy1, thisx1), (thisy2, thisx2) )
for thisx1, thisy1, thisx2, thisy2 in zip( x1,y1,x2,y2)]
coll = LineCollection(segments)
gca().add_collection(coll)
savefig('kontour')
show()
|
|
From: Helge A. <av...@ii...> - 2004-10-01 16:52:23
|
John Hunter <jdh...@ac...> writes:
|
| Do you have any thoughts on how we might do labels with your code?
2 ways:
1) automatically: define a coarser(user defined coarseness) mesh on
the function to be labelled, and add the labels in the vertices of
this mesh. the angle of each label can easily be found from the
closest contour line segment. this method avoids clusters of labels
effectively, but will not be good in areas with high variability.
2) manually: let the user point and click on contours, I implemented
this for use with gist that you could take a look at (clabel below).
matlab also does this, and I think this is the best way for real
publication quality. I have never seen automatic routines that do
labelling well. TECPLOT is close but not quite there.
I also attached a routine to do this for 2D stretched coordinates (see
the contour plots with labels on
http://www.ii.uib.no/~avle/python.html
for examples) which is called vclabel
|
| If we decide to go with your routines, at least for the time being
| until we can "do it right", would you be willing to contribute your
| code to matplotlib under the matplotlib license (PSF inspired, free
| for commercial and noncommercial reuse)?
sure, no problem.
Helge
def clabel(z,clevels,opa=1,col='black',meth='std',digits=1):
"""
clabel(z,clevels,opa=1,col='black')
At the point where the mouse is clicked, print the contour level
from clevels that are closest to the interpolated value in this
point. Meant to be useful for labeling contour lines manually...
optional arguments:
opa=0 : Transparent text.
opa=1 : Erase background under label
col='white' : text color.
meth= 'std' bilinear interpo on a std grid, might give error near coast
'grid' values not given cellcentered but on corners of cells.
'bigrid' values taken from a bilinearly interpol fine mesh from
futil.contour
'cell' takes value from cell directly
Helge Avlesen <av...@ii...>
"""
print """
Insert contour levels by left clicking, middle button display
values, right button finishes.
"""
button=0
while button<>3:
mus=gist.mouse()
button=mus[9]
if meth=='bigrid':
i=int(2*mus[0])+1 ; j=int(2*mus[1])+1
x=2*mus[0]-i+1 ; y=2*mus[1]-j+1
elif meth=='std':
i=int( mus[0] ); j=int( mus[1] )
x=mus[0]-i ; y=mus[1]-j
elif meth=='grid':
xm=mus[0]+0.5 ; ym=mus[1]+0.5
i=int( xm ); j=int( ym )
x=xm-i ; y=ym-j
elif meth=='cell':
print mus[0], mus[1]
i=int(round(mus[0])) ; j=int(round(mus[1]))
if meth=='cell':
val=z[i,j]
print val,i,j
else:
# bilinear interpolation to find value
a00=z[i,j]
a10=z[i+1,j]-a00
a01=z[i,j+1]-a00
a11=z[i+1,j+1]-(a00+a10+a01)
val=a00 + a10*x + a01*y + a11*x*y
print val,i,j,x,y
if button==1:
# compare this value to the selected levels
diff= abs( clevels-val )
# use the closest
label=fpformat.fix( clevels[ Numeric.argmin(diff) ], digits )
gist.plt(label, mus[0], mus[1], opaque=opa, tosys=1, \
height=8, justify="CH", color=col )
def vclabel(z,sx,sy,clevels,opa=1,col='black',digits=1):
"""
manual(z,clevels,opa=1,col='black')
At the point where the mouse is clicked, print the contour level
from clevels that are closest to the interpolated value in this
point. Meant to be useful for labeling contour lines manually...
sx[i,j],sy[i,j] is the x,z coordinate of point z[i,j]. if [:,1]
denotes the top layer, [:,kb-1] the bottom (common in
oceanography) j_is_down will be true. (z is always positive in
the upward direction, but the indice j may go downwards)
optional arguments:
opa=0 : Transparent text.
opa=1 : Use background color for text.
col='white' : text color.
digits: number of decimals in label
Helge Avlesen <av...@ii...>
"""
print """
Insert contour levels by left clicking, middle button display
depth, right button finishes.
"""
kb=z.shape[1]
im=z.shape[0]
j_is_down=0
if sy[0,0]>sy[0,1]:
j_is_down=1
button=0
while button<>3:
mus=gist.mouse()
button=mus[9]
# bisection search for the indices
i=hbisect( sx[:,0], mus[0] )
if i<0 or i>im-1:
print 'outside:',i
continue
x=(mus[0]-sx[i,0])/(sx[i+1,0]-sx[i,0])
if j_is_down:
finn=hbisect( (1.-x)*sy[i,::-1] + x*sy[i+1,::-1] , mus[1] )
j=kb-2-finn
if finn<0 or finn>kb-1:
print 'outside',i,finn
continue
xa=Numeric.array((sx[i,j+1]+x*(sx[i+1,j+1]-sx[i,j+1]),\
sy[i,j+1]+x*(sy[i+1,j+1]-sy[i,j+1])))
if mus[1]-xa[1]<0:
print 'below'
continue
xb=Numeric.array((sx[i,j]+x*(sx[i+1,j]-sx[i,j]),\
sy[i,j]+x*(sy[i+1,j]-sy[i,j])))
y=(((mus[0]-xa[0])**2 + (mus[1]-xa[1])**2 )\
/((xb[0]-xa[0])**2 + (xb[1]-xa[1])**2 ))**0.5
# bilinear interpolation to find value
a1=z[i,j+1]
a2=z[i+1,j+1]-a1
a3=z[i,j]-a1
a4=z[i+1,j]-(a1+a2+a3)
val=a1 + a2*x + a3*y + a4*x*y
else:
print 'increasing j upwards not yet implemented'
continue
print 'x,y=',x,y,' i,j=',i,j, 'val=',val
if button==1:
# compare this value to the selected levels
diff= abs( clevels-val )
# use the closest
label=fpformat.fix( clevels[ Numeric.argmin(diff) ], 1)
if opa==1:
label=' '+label+' '
gist.plt(label, mus[0], mus[1], opaque=opa, tosys=1, \
height=8, justify="CH", color=col )
|
|
From: John H. <jdh...@ac...> - 2004-10-01 16:37:29
|
>>>>> "Perry" == Perry Greenfield <pe...@st...> writes:
>> Hi again,
>>
>> Yes, I had the thought that using C for the algorithm would be
>> easier as well. There are actually some very well-written
>> marching squares contouring algorithms in C already out there.
>> I will try to find such an implementation and point you to it
>> or send you the source code.
>>
Perry> Thanks, that would be helpful. In my search I didn't come
Perry> across many. Keep in mind the license needs to be
Perry> compatible with that of matplotlib.
Of course, in addition to the license, there is the patent issue. I
believe marching squares is patented. I know marching cubes is.
http://patft.uspto.gov/netacgi/nph-Parser?Sect1=PTO1&Sect2=HITOFF&d=PALL&p=1&u=/netahtml/srchnum.htm&r=1&f=G&l=50&s1=4,710,876.WKU.&OS=PN/4,710,876&RS=PN/4,710,876
I checked the header of vtkMarchingSquares.cxx which states
Program: Visualization Toolkit
Module: $RCSfile: vtkMarchingSquares.cxx,v $
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
THIS CLASS IS PATENTED UNDER UNITED STATES PATENT NUMBER 4,710,876
"System and Method for the Display of Surface Structures Contained
Within the Interior Region of a Solid Body".
Application of this software for commercial purposes requires
a license grant from GE. Contact:
but the patent number they reference which is linked above begins with
A method and apparatus for displaying *three dimensional surface
images* includes the utilization of a case table for rapid retrieval
of surface approximation information.
emphasis mine. So I don't know for sure what the patent status of the
2D algorithm is.
JDH
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X