 # Contour and Vector Plots

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)?)

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!

> 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.

Do you have any thoughts on how we might do labels with your code?

> 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.

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)
savefig('kontour')

show()

John Hunter <jdhunter@...5...> 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 <avle@...153...>

"""

print """
Insert contour levels by left clicking, middle button display
values, right button finishes.
"""

button=0
while button<>3:
mus=gist.mouse()
button=mus

if meth=='bigrid':
i=int(2*mus)+1 ; j=int(2*mus)+1
x=2*mus-i+1 ; y=2*mus-j+1
elif meth=='std':
i=int( mus ); j=int( mus )
x=mus-i ; y=mus-j
elif meth=='grid':
xm=mus+0.5 ; ym=mus+0.5
i=int( xm ); j=int( ym )
x=xm-i ; y=ym-j
elif meth=='cell':
print mus, mus
i=int(round(mus)) ; j=int(round(mus))

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, mus, 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 <avle@...153...>

"""

print """
Insert contour levels by left clicking, middle button display
depth, right button finishes.
"""
kb=z.shape
im=z.shape

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

# bisection search for the indices
i=hbisect( sx[:,0], mus )
if i<0 or i>im-1:
print 'outside:',i
continue

x=(mus-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 )
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-xa<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-xa)**2 + (mus-xa)**2 )\
/((xb-xa)**2 + (xb-xa)**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, mus, opaque=opa, tosys=1, \
height=8, justify="CH", color=col )