Re: Q: contour levels from IDL

[hcp(at)newsread.ed.ac.uk (H C Pumphrey)]

In article <8b5pfs$m7e$1@nnrp1.deja.com>, Mirko Vukovic <mvukovic@taz.telusa.com> writes:

|> My problem is how to choose intelligent levels.  The only way I could
|> think of is to plot (2dplot) the z component of the data in a temp
|> (ram?) window,  and observe the z ticks used by IDL, and use those as
|> contour levels.

Faced with this, I wrote my own contour chooser which I enclose below. 
You give it a max, min and a number of contours and you get back an array 
with values whose last digits are like 1,2,3,4 or 2,4,6,8 or 
5,10,15,20,25. I like it. YMMV. Sorry about the commented-out
debugging statements.


Enjoy

Hugh


function conpick,minval,maxval,ndesired
; This function returns a set of contours suitable for the range of
; values beween max and min. on input ndesired is the number of contours
; you want; you actually get between ndesired and 2*ndesired
; contours. 


;*** Work out a nice spacing for the contours *** 

if maxval lt minval then begin
    tmp=maxval
    maxval=minval
    minval=tmp
endif
range=double(maxval)-double(minval) 
;print,'range=',range
initspace=range / (ndesired-1)
;print,'initspace=',initspace
reduspace=alog10(initspace)
;print,'reduspace=',reduspace
pwroften=long(reduspace)
if reduspace lt 0.0 then pwroften=pwroften-1
;print,'pwroften=',pwroften
reduspace=reduspace-pwroften
;print,'reduspace=',reduspace
reduspace=10.^reduspace
;print,'reduspace=',reduspace,' Choosing spacing'

if(reduspace ge 5.0) then spacing=5.0 $
;else if(reduspace ge 2.5) then spacing=2.5$
else if(reduspace ge 2.0) then spacing=2.0$ 

else spacing=1.0
;print,'spacing=',spacing
spacing=spacing*10.^pwroften
;print,'spacing=',spacing

if( (not finite(spacing)) or (not finite(reduspace))) then begin
    print,'Error in routine conpick: failed to pick suitable contours'
    return,findgen(ndesired)
endif

; *** Work out a minimum value that fits in with the spacing *** 
gak=(abs(minval)/(10.0^(pwroften+1)))
igak=long(gak)
if gak lt 0.0 then igak=igak-1 
newminval=(10.0^(pwroften+1))* igak
;print,'newminval=',newminval

while(newminval lt abs(minval)) do newminval=newminval+spacing
;print,'newminval=',newminval

if(minval lt 0.0) then newminval=-newminval else newminval=newminval-spacing
;print,'newminval=',newminval,' Making contours'

upper=newminval
ncontours=0
while ( upper lt maxval ) do begin
    upper=newminval+ncontours*spacing
    ncontours=ncontours+1
endwhile
;print,' made',ncontours,' Contours'
clevs=findgen(ncontours)*spacing + newminval

return,clevs

end




-- 

==========================================================================
Hugh C. Pumphrey              |  Telephone 0131-650-6026
Department of Meteorology     |  FAX       0131-650-5780
The University of Edinburgh   |  Replace 0131 with +44-131 if outside U.K.
EDINBURGH EH9 3JZ, Scotland   |  Email     hcp@met.ed.ac.uk
OBDisclaimer: The views expressed herein are mine, not those of UofE.
==========================================================================