Contour and Vector Plots

Development matplotlib-devel
1

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)
gca().add_collection(coll)
savefig('kontour')

show()

2

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[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 <avle@...153...>
    
    """
    
    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 )