;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; plot_binned.pro ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Author: wd (Wolfgang.Dobler@ncl.ac.uk)
;;; $Date: 2007-08-16 23:22:43 $
;;; $Revision: 1.9 $
;;; Description:
;;; Scatter-plot data, but bin them first in order to reduce the
;;; size of thusly created PostScript files. For 60^3 data points,
;;; this can easily result in a 5 times smaller file.
;;; Alternatively, the density of the binned data can be displayed
;;; using PLOTIMAGE (recommended) or CONTOUR
;;;
;;; Usage: like `plot', but with the additional key words
;;; NX=nx, NY=ny --- set vertical and horizontal number of points
;;; /CONTOUR --- do contour plot (by default filled)
;;; /FILL --- /CONTOUR, FILL=0 does contour lines
;;; NLEVELS=nlevels --- set NLEVELS for /CONTOUR (default=60)
;;; /PLOTIMAGE --- use PLOTIMAGE to display density
;;; (faster than CONTOUR and PS files are
;;; much smaller)
;;; /INVERT --- invert data for /CONTOUR (default is
;;; white for no data, dark for dense data
;;; with standard colour maps)
;;; /OVERPLOT --- do overplotting (kind of `oplot_binned')
;;; /EQUX, /EQUY --- equalize in x/y, i.e. normalize
;;; vertical/horizontal sum of counts to
;;; unity (except where there is no
;;; point). Useful with /PLOTIMAGE or
;;; /CONTOUR when plotting over radius
;;; data in cubic box
;;; ENHANCE=enhance --- multiply density map with ENHANCE,
;;; which must be a non-negative 2d array
;;; of size NX x NY
;;; ASINHSCALE=scale --- plot asinh(scale*density) instead of density
;;; DENSITY=density --- return density map (only values 0 an 1
;;; unless /CONTOUR was sepcified
;;; XCOORDS=xcoords --- corresponding x coordinates
;;; YCOORDS=ycoords --- corresponding x coordinates
;;;
;;; Some key words like PSYM, SYMSIZE and COLOR are intecepted, such as
;;; to work in overplot mode as well as in normal plot mode. For the
;;; rest, draw the frame first explicitly and then call plot_binned
;;; with /OVERPLOT.
;;;
;;; Examples:
;;; N=8000 & xx=1+5.*randomu(seed,N) & ff=2.5+cos(4*xx)+0.2*randomn(seed,N)
;;; plot_binned, xx, ff
;;; plot_binned, xx, ff, NX=60, NY=40, /PLOTIMAGE, /XLOG, /YLOG
;;; plot_binned, xx, ff, NX=60, NY=40, /CONTOUR
;;; plot_binned, xx, ff, NX=60, NY=40, /CONTOUR, /INVERT
;;; plot_binned, xx, ff, NX=60, NY=40, /PLOTIMAGE, /EQUX
;;; plot_binned, xx, ff, /PLOTIMAGE, /EQUX, ASINHSCALE=5.
;;;
;;; To do:
;;; Reduce memory usage: xx,yy,good,xind,yind all use the same
;;; amount of space as the original data xvar, yvar
;;; [Partly accomplished]
pro plot_binned, xvar, yvar, $
NX=nmx, NY=nmy, $
CONTOUR=contour, $
FILL=fill, NLEVELS=nlevels, $
PLOTIMAGE=plotimage, $
INVERT=invert, $
OVERPLOT=overplot, $
PSYM=psym, SYMSIZE=symsize, $
COLOR=color, $
XRANGE=xrange, YRANGE=yrange, $
XSTYLE=xstyle, YSTYLE=ystyle, $
EQUX=equx, EQUY=equy, $
ENHANCE=enhance, $
ASINHSCALE=asinhscale, $
DENSITY=map, XCOORDS=mapx, YCOORDS=mapy, $
HELP=help, $
_EXTRA=extra
if (keyword_set(help)) then extract_help, 'plot_binned'
if (n_elements(nmx) eq 0) then nmx = 500/(!p.multi[1] > 1)
if (n_elements(nmy) eq 0) then nmy = 350/(!p.multi[2] > 1)
if (n_elements(psym) eq 0) then psym = 3
if (n_elements(contour) eq 0) then contour = 0
if (n_elements(plotimage) eq 0) then plotimage = 0
if (n_elements(invert) eq 0) then invert = 0
if (n_elements(fill) eq 0) then fill = 1
if (n_elements(nlevels) eq 0) then nlevels = 60
if (n_elements(equx) eq 0) then equx = 0
if (n_elements(equy) eq 0) then equy = 0
if (n_elements(enhance) eq 0) then enhance = -1
if (n_elements(asinhscale) eq 0) then asinhscale = 0
if (n_elements(xrange) eq 0) then xrange = [0,0]
if (n_elements(yrange) eq 0) then yrange = [0,0]
if (n_elements(xstyle) eq 0) then xstyle = !x.style
if (n_elements(ystyle) eq 0) then ystyle = !y.style
;; Sanitize keywords
if (plotimage ne 0) then contour=0
if ((plotimage ne 0) or (contour ne 0)) then scatter=0 else scatter=1
if ((equx ne 0) and (equy ne 0)) then $
message, '/EQUX and /EQUY keywords are mutualy exclusive'
;; Save so we can reset later to avoid plotting to a new frame
pmulti1 = !p.multi
;; Plot the frame (necessary to get x- and y-range)
if (n_elements(overplot) eq 0) then begin
if (scatter eq 0) then begin ; don't plot and avoid extra PS frame
xstyle1 = (xstyle or 4)
ystyle1 = (ystyle or 4)
endif else begin
xstyle1 = xstyle
ystyle1 = ystyle
endelse
plot, xvar, yvar, /NODATA, $
XRANGE=xrange, YRANGE=yrange, $
XSTYLE=xstyle1, YSTYLE=ystyle1, $
COLOR=color, _EXTRA=extra
pmulti2 = !p.multi ; save so we can exit with correct !p.multi
endif
xr = !x.crange ; NB: for /xlog, this is alog10(data_range)
yr = !y.crange
dmx = (xr[1]-xr[0])*1./Nmx
dmy = (yr[1]-yr[0])*1./Nmy
;; Handle /xlog and /ylog
if (!x.type ne 1) then begin
xx=xvar
mapx = linspace(xr,Nmx)
xbin = xr[0] + (findgen(Nmx)+0.5)*dmx
imgxrange = xr
endif else begin
xx=alog10(xvar)
mapx = logspace(xr,Nmx)
xbin = 10^(xr[0] + (findgen(Nmx)+0.5)*dmx)
imgxrange = 10^xr
endelse
;
if (!y.type ne 1) then begin
yy=yvar
mapy = linspace(yr,Nmy)
ybin = yr[0] + (findgen(Nmy)+0.5)*dmy
imgyrange = yr
endif else begin
yy=alog10(yvar)
mapy = logspace(yr,Nmy)
ybin = 10^(yr[0] + (findgen(Nmy)+0.5)*dmy)
imgyrange = 10^yr
endelse
map = hist_2d(xx, yy, $
MIN1=xr[0], MAX1=xr[1], BIN1=dmx, $
MIN2=yr[0], MAX2=yr[1], BIN2=dmy) * 1.D0
;; hist and hist_2d do really bizarre things to the last bin. Here
;; this just means we should remove the last entry in each
;; direction:
map = map[0:NMx-1,0:NMy-1]
;; Equalize if requested
if (equx ne 0) then begin
for ix=0,NMx-1 do begin
sum = total(1.D0*abs(map[ix,*]))
if (sum gt 0) then map[ix,*] = map[ix,*]/sum
endfor
endif
;
if (equy ne 0) then begin
for iy=0,NMy-1 do begin
sum = total(1.D0*abs(map[*,iy]))
if (sum gt 0) then map[*,iy] = map[*,iy]/sum
endfor
endif
;; Multiply with ENHANCE, if requested
if (enhance[0] ge 0) then begin
s = size(enhance)
if ((s[0] ne 2) or (s[1] ne NMx) or (s[2] ne NMy)) then $
message, 'ENHANCE must be of size '+strtrim(NMx,2)+'x'+strtrim(NMy,2)
map = map*enhance
endif
;; Apply asinh if requested
if (asinhscale ne 0) then begin
map = map / max(abs(map))
map = asinh(asinhscale*map)
endif
if (scatter eq 0) then !p.multi = pmulti1
;; Plot
if (scatter ne 0) then begin ; scatter plot
;; Re-use xx, yy here to save memory
;; 2-d x- and y-coordinates for the map
xx = spread(indgen(NMx),[1],[NMy])
yy = spread(indgen(NMy),[0],[NMx])
goodm = where(map ne 0)
;
if (goodm[0] ge 0) then begin
oplot, xbin[xx[goodm]], ybin[yy[goodm]], $
PSYM=psym, SYMSIZE=symsize, COLOR=color, $
_EXTRA=extra
endif else begin
message, /INFO, 'No valid points in domain' $
+ '[' + wdtrim(xr[0]) + ', '+ wdtrim(xr[1]) +']'$
+ ' x ' $
+ '[' + wdtrim(yr[0]) + ', '+ wdtrim(yr[1]) +']'
endelse
endif else begin
if (invert ne 0) then begin
levels = linspace(0.,max(map),NLEVELS,GHOST=0.5)
fact = 1.
endif else begin
levels = linspace(-max(map),0,NLEVELS,GHOST=0.5)
fact = -1.
endelse
if (contour ne 0) then begin ; contour plot
contour, fact*map, mapx, mapy, $
/OVERPLOT, $
FILL=fill, LEVELS=levels, $
_EXTRA=extra
;; Replot the frame as CONTOUR will have erased most of it
if (n_elements(overplot) eq 0) then begin
plot, xvar, yvar, /NODATA, COLOR=color, /NOERASE, $
XRANGE=xr, YRANGE=yr, $
XSTYLE=((xstyle or 1) and 13), $ ; set 1, unset 2
YSTYLE=((ystyle or 1) and 13), $
_EXTRA=extra
!p.multi = pmulti2
endif
endif else begin ; plotimage
plotimage, fact*map, RANGE=minmax(fact*map), $
IMGXRANGE=imgxrange, IMGYRANGE=imgyrange, $
XSTYLE=((xstyle or 1) and 13), $ ; set 1, unset 2
YSTYLE=((ystyle or 1) and 13), $
_EXTRA=extra
endelse
endelse
end
; End of file plot_binned.pro