# Re: Polar Plots

```Hi Chris,

I have written some code to plot what I was calling a wind rose
(It's not a real wind rose, just a routine to plot a value against a
direction, which in this case is wind direction). Anyway, you might
want to have a look at the code and modify it for your needs. Note
that 0 degrees in this plot is at "north" or top y-axis whereas in
true polar coords 0 degrees is right hand x-axis. The code does
a lot of other things that you are probably not interested in eg
can plot coastlines behind the polar plot.
Find the code plot_wind_rose.pro attached. Hope this helps!

Cheers, Paul

"Chris Bull" <cjbull@another.com> wrote in message
sCmF6.8777\$Ln6.1258097@news2-win.server.ntlworld.com">news:sCmF6.8777\$Ln6.1258097@news2-win.server.ntlworld.com...
> Hi all,
>
>
> I was just wondering if anyone had extended (and is willing to share) the
> IDL basic polar plot routine to include useful things like circular axes
> etc (before I slave over doing it myself! :*)
>
> Regards
>
>
> Chris Bull
>
>

FUNCTION CIRCLE, xcenter, ycenter, radius, nPts=nPts

IF N_Elements(nPts) EQ 0 THEN nPts = 100
points = (2 * !PI / (nPts-1)) * FIndGen(nPts)
x = xcenter + radius * Cos(points)
y = ycenter + radius * Sin(points)
RETURN, Transpose([[x],[y]])

END
;+
; NAME:
; PLOT_WIND_ROSE
;
; PURPOSE:
; This procedure will plot a given variable against
; wind direction in a polar plot. This type of plot
; is commonly known as a wind rose in meteorology.
; The routine has several options including the option
; of plotting the rose centred over a particular map
; location. Other options inlcude placing baseline
; indicators, any number of circles and cross hairs
; onto the plot.
;
; CATEGORY:
; Plotting/meteorology.
;
; CALLING SEQUENCE:
; PLOT_WIND_ROSE, Var, Wind_Dirn
;
; INPUTS:
; Var: This is the variable to plot on the wind rose.
;  This is an array of the same lenght as Wind_Dirn
;  corresponding to each element in Wind_Dirn.
; Wind_Dirn: This is the wind direction in degrees.
;
; KEYWORD PARAMETERS:
; TITLE: A string that contains a title for the plot.
; COLOR: This controls the colour of the wind rose line.
;  Default color is 0.
; THICK: This controls the thickness of the wind rose
;  line. Default thickness is 4.
; LINESTYLE: This controls the linestyle of the wind rose
;  line. Default linestyle is 0.
;
; MAXV: Set this keyword to a scalar of the maximum
;  value of the plot axis. Positive scalar only.
;
; VAR2: Set this to an array to plot a second wind rose.
;  This must be an array of the same lenght as WIND_DIRN2
;  corresponding to each element in WIND_DIRN2.
; WIND_DIRN2: This is the second wind direction in degrees.
; V2_COLOR: This controls the colour of the 2nd wind rose line.
;  Default color is 0.
; V2_THICK: This controls the thickness of the 2nd wind rose
;  line. Default thickness is 4.
; V2_LINESTYLE: This controls the linestyle of the 2nd wind rose
;  line. Default linestyle is 1.
;
; YTITLE: A string that contains a title for the y-axis.
; YTHICK: This controls the thickness of the y-axis.
;  Default thickness is 2.
; YMINOR: This controls the number of minor tickmarks.
;  Default is 5.
;
; CROSS_HAIRS: Set this keyword to place tick marks on
;  the horizontal and vertical radial lines through
;  the centre of the plot.
; E_CROSS: Set this keyword to a structure that will
;  contain any keywords used in the AXIS procedure.
;  The AXIS procedure is used to draw the cross hairs.
;
;  lines ("spokes") that will be plotted and labelled.
;  If not set at all then the default is to plot 12 radial
;  lines (ie every 30 degrees). Set N_RAD_LINES to zero to
;  plot NO radial lines. If set to zero and CROSS_HAIRS is
;  selected then the cross hair angles are lebelled (ie 0,
;  90, 180 and 270 degrees).
; E_RAD_LINES: Set this keyword to a structure that will
;  contain any keywords (linestyle, color, thickness etc)
;  used in the OPLOT procedure for plotting the radial
;  lines. Default is same as for OPLOT and THICK=1.0.
;  NOTE: If CROSS_HAIRS is selected and the angles
;  for the radial lines coincide with the cross hairs
;  angles then the radial lines are plotted OVER the
;  cross hair lines!
;
; BASELINE: Set this keyword if baseline indicator lines
;  are required. The default is plot the baseline
;  indicators for Cape Grim, Tasmania ie lines at 190
;  and 280 degrees. To change this, set this keyword to
;  an array of the angles required ie BASELINE=[150.,290.]
;  This array can contain any number of angles as long
;  as there is more than 1 angle.
; E_BASE: Set this keyword to a structure that will contain
;  any keywords used in the OPLOT procedure ie linestyle,
;  color, thickness etc. Defaults are the same as used
;  in OPLOT.
;
; MAP: Set this keyword to place a map on the wind rose.
;  The default is to plot a high resolution map with the
;  COASTS keywrod set centred on Cape Grim, Tasmania
;  (~-40.68,~144.68). To centre the map on a different
;  location set the keyword to a two element array
;  containing the latitude and longitude of the desired
;  point ie MAP=[-38.03,145.10]. The map is setup with
;  MAP_SET and drawn with MAP_CONTINENTS.
; E_MAP: Set this keyword to a structure that will contain
;  any keywords used in the MAP_SET procedure. Use this
;  to set the map scale or projection. The map scale
;  controls what the map region is. The default is 5.e6
;  so a map of scale 1:5e6 is plotted centred on the given
;  lat and lon. DO NOT use this to set actual plotting
;  features ie color, thickness, fill etc Use the E_CONT
;  keyword described below.
; E_CONT: Set this keyword to a structure that will contain
;  any keywords used in the MAP_CONTINENTS procedure. Use
;  this to set plotting variables of the map such as line
;  thickness, color, turn high res/coasts off, etc.
; GRID: Set this keyword to plot a grid on the map using
;  MAP_GRID.
; E_GRID: Set this keyword to a structure that will contain
;  any keywords used in the MAP_GRID procedure. Use this
;  keyword to set the linestyle, thickness, color, labels
;  etc for the grid. Refer to MAP_GRID.
;
; N_CIRC: Set this keyword to an array containing values
;  between zero and one indicating where to plot the
;  circles on the wind rose. The default if this is not
;  set is to plot two circles, one at 0.5 of the maximum
;  value, the other at 1. If you wanted to plot 4 circles
;  set this keyword to where the circles are required ie
;  N_CIRC=[0.25, 0.5, 0.75, 1.0].
; E_CIRC: Set this keyword to a structure that will contain
;  any keywords used in the PLOTS procedure which draws
;  the circles. Use this to control the color, thicknes
;  and linestyle of the circles. The defaults for these
;  are the same as for the PLOTS procedure.
;
; COPY_FIRST: Set this keyword to copy the the first array
;  element of the input arrays VAR and WIND_DIRN to the
;  end of these arrays. This is to "close" the wind
;  rose line.
;
; NODATA: Set this keyword to plot no data. This is useful
;  for visualising wind directions for a particuar
;  location by plotting a map in the background with the
;  windrose angles or coordinates over the top.
;
; EXAMPLE:
; To plot a wind rose of the variable ch4conc with a map
; centred over Cape Grim, with scale of 10e6, no grid lines
; on the map, 4 circles, cross hairs and the default
; baseline indicators:
; PLOT_WIND_ROSE, ch4conc, wind_dir, \$
;  color=pen(2), thick=6, linestyle=0,  \$
;  title='Methane at Cape Grim', \$
;  ytitle='CH!D4!N (ppb)', /cross_hairs, \$
;  map=1, E_MAP={scale:10e6}, \$
;  E_cont={thick:4,color:pen(4)}, \$
;  /baseline, \$
;  E_BASE={thick:4, color:pen(3), linestyle:2},  \$
;  E_CIRC={thick:2}, n_circ=[0.25,0.5,0.75,1.0]
;
; To plot the above with filled continents and a grid
; on the map:
; PLOT_WIND_ROSE, ch4conc, wind_dir, \$
;  color=pen(2), thick=6, linestyle=0,  \$
;  title='Methane at Cape Grim', \$
;  ytitle='CH!D4!N (ppb)', /cross_hairs, \$
;  map=1, E_MAP={scale:10e6}, \$
;  E_cont={thick:4,color:pen(4),fill_continents:1}, \$
;  /baseline, \$
;  E_BASE={thick:4, color:pen(3), linestyle:2},  \$
;  E_CIRC={thick:2}, n_circ=[0.25,0.5,0.75,1.0], \$
;  /GRID, E_GRID={thick:0.5}
;
; To plot just a simple wind rose:
; PLOT_WIND_ROSE, ch4conc, wind_dir, ytitle='CH!D4!N (ppb)'
;
; MODIFICATION HISTORY:
; Written by: Bronwyn Dunse and Paul Krummel,
;  27 October 1998, CSIRO Atmospheric Research, Australia.
; Modified by: Paul Krummel, 8 November 1998. Added map
;  functionality.
; Modified by: Paul Krummel, 15 January 1999. Added many
;  extra keywords (E_???), cleaned up the routine, added
;  help and proper header information, also other small
;  improvements and bug fixes.
; MODIFIED by: Paul Krummel, 21 January 1999. Fully documented.
; MODIFIED by: Paul Krummel, 10 February 1999. Added
;  NODATA keyword.
; Modified by: Paul Krummel, 11 May 1999. Fixed Title problem
; Modified by: Paul Krummel, 13 December 1999. Added VAR2, WIND_DIRN2,
;  V2_COLOR, V2_THICK and V2_LINESTYLE keywords.
; Modified by: Paul Krummel, 22 February 2000. Changed way the radial
;  lines and their annotations (degrees) were plotted. It is now
;  much easier to change the number of radial lines plotted.
;  and E_RAD_LINES to control the number and style of the radial lines.
; Modified by: Paul Krummel, 25 July 2000. Fixed bug introduced during
;  23 February 2000 changes above. Now, if N_Rad_Lines is not set at all
;  then it is set to 12, this bombed out before with variable undefined!
; Modified by: Paul Krummel, 11 October 2000. Fixed small bug where by there
;  was a slight offset in the 0 degrees label, due to floating point
;  inaccuracies. In the checking of what alignment to use for each label,
;  now take 'fix' of the angle value to give integer numbers.
;
;-
;
PRO PLOT_WIND_ROSE, VAR, WIND_DIRN, \$
VAR2=Var2, WIND_DIRN2=Wind_dirn2, \$
MAXV=Maxv, \$
TITLE=Title, COLOR=Color, THICK=Thick, \$
LINESTYLE=Linestyle, \$
V2_COLOR=V2_Color, V2_THICK=V2_Thick, \$
V2_LINESTYLE=V2_Linestyle, \$
YTITLE=YTitle, YTHICK=YThick, YMINOR=YMinor, \$
CROSS_HAIRS=Cross_Hairs, E_CROSS=E_Cross, \$
BASELINE=Baseline, E_BASE=E_Base, \$
MAP=Map, E_MAP=E_Map, E_CONT=E_CONT, \$
GRID=Grid, E_GRID=E_Grid, \$
N_CIRC=N_Circ, E_CIRC=E_Circ, \$
COPY_FIRST=Copy_First, NODATA=Nodata

; ++++
; =====>> HELP
on_error,2
if (N_PARAMS(0) LT 2) or keyword_set(help) then begin
doc_library,'PLOT_WIND_ROSE'
if N_PARAMS(0) LT 2 and not keyword_set(help) then \$
message,'Need at least two input parameters, see above for usage.'
return
endif

; ++++
; Find the maximum value of the variable of interest.
; If the NODATA keyword is set then set max_var to 1.0
; and zero the var[] array.
max_var = keyword_set(nodata) ? 1.0 : max(var)
max_var = (n_elements(var2) gt 0) ? max([max_var,var2]) : max_var
if keyword_set(nodata) then var[*]=0.0
;
; If MAXV keyword is set then use its value for max_var.
if n_elements(Maxv) eq 1 then max_var=Maxv
;
; If the COPY_FIRST keyword is set then copy the first
; array element of VAR and WIND_DIRN to the end of these
; arrays.
if keyword_set(Copy_First) then begin
var=[var,var[0]] & wind_dirn=[wind_dirn,wind_dirn[0]]
if n_elements(var2) gt 0 then begin
var2=[var2,var2[0]] & wind_dirn2=[wind_dirn2,wind_dirn2[0]]
endif
endif
;
; Check if the N_Rad_Lines keyword was set, if not default it to 12.
;
; ++++
; Setup the plot coordinates. Make the plot isotropic and
; set the x and y range to +/- the maximum value.
PLOT, var, (450.-wind_dirn)*!DTOR, /nodata, /noerase,\$
/polar, /isotropic,\$
xstyle=5, ystyle=5, \$
yrange=[-max_var,max_var], xrange=[-max_var,max_var]

; ++++
; If the MAP keyword is selected then plot a map behind
; the wind rose centred on the given coordinates. If no
; coordinates are given use lat and lon for Cape Grim,
; Tasmania, Australia.
if keyword_set(MAP) then begin

; Check if there is more than one element in the map
; keyword, if so it should contain the lat and lon
; for the map centre. Set this accordingly.
if n_elements(map) gt 1 then begin
lat=map[0] & lon=map[1]
endif else begin  ; else set the default map centre
; to Cape Grim.
lat=-40.6829441667D & lon=144.689568611D
endelse

; Check how many plots per page are to be made.
; Use different code for just one plot or more
; than one plot and adjust !p.multi accordingly.
if !p.multi[1] gt 1 or !p.multi[2] gt 1 then begin

;+ Code for more than one plot per page +
num_left=!p.multi[0]
if num_left eq 0 then num_left=!p.multi[1]*!p.multi[2]

; Setup the mapping here, forcing the map
; into the plot coordinates set above.
map_set, lat, lon, scale=5.e6, /merc, /noborder, \$
position=[!x.window[0],!y.window[0], \$
!x.window[1],!y.window[1]], \$

; Plot the coastline here, NOTE the default
; is to use hires map.
map_continents, /coasts, /hires, _EXTRA=E_CONT

; Plot a grid on the map if requested.
if KEYWORD_SET(GRID) then map_grid, _EXTRA=E_GRID

; reset !p.multi to stop frame/page advance.
!p.multi[0]=num_left

endif else begin
;+ Code for one plot per page +

; Setup the mapping here, forcing the map
; into the plot coordinates set above.
map_set, lat, lon, /merc, scale=5.e6, /noborder, \$
position=[!x.window[0],!y.window[0], \$
!x.window[1],!y.window[1]], \$
_EXTRA=E_MAP

; Plot the coastline here, NOTE the default
; is to use hires map.
map_continents, /coasts, /hires, _EXTRA=E_CONT

; Plot a grid on the map if requested.
if KEYWORD_SET(GRID) then map_grid, _EXTRA=E_GRID

; reset !p.multi to stop frame/page advance.
if !p.multi[0] le 0 then !p.multi[0]=1

endelse

endif ; End of mapping section.

; ++++
; Setup the plot coordinates again and get the tickmark values.
; Set the title and add !C (new line) to the end of it.
Title = keyword_set(Title) ? Title+'!C  ' : Title
PLOT, var, (450.-wind_dirn)*!DTOR, /nodata, \$
/polar, /isotropic, \$
TITLE=Title, \$
xstyle=5, ystyle=5, \$
yrange=[-max_var,max_var],xrange=[-max_var,max_var],\$
ytick_get=tick_val

; ++++
; Plot the circles here.
; Default circles at 0.5 and 1 of Max_var.
if not KEYWORD_SET(n_circ) then n_circ=[0.5,1.0]
for i=0,n_elements(n_circ)-1 do  \$
PlotS, circle(0, 0, n_circ[i]*Max_var, nPts=361), _EXTRA=E_CIRC

; ++++
; Count the number of tickmarks
ntick_val=n_elements(tick_val)

; Set the negative tickmarks to positive, for use in labelling.
ptick_val=abs(tick_val)

; Set the default thickness and number of minor tick marks
; if they were not defined in the calling routine.
ythick = keyword_set(ythick) ? ythick : 2
yminor = keyword_set(yminor) ? yminor : 5

; ++++
; Plot a y-axis to the left of the rose with the tick mark values
; and optional ytitle. Default thickness and minor tickmarks are
; used unless otherwise specified.
AXIS,-(max_var+0.22*max_var),0, YAXIS=0, YTHICK=YThick, \$
YTICKNAME=format_axis_values(ptick_val), \$
YTICKS=ntick_val-1, YTICKV=tick_val, \$
YTITLE=YTitle, YMINOR=YMinor, \$
YSTYLE=1, YRANGE=[-max_var,max_var]

; ++++
; CROSS HAIRS:
; Place tick marks on the horizontal and vertical radial
; lines through the centre of the plot if the CROSS_HAIRS
; keyword is set.
if KEYWORD_SET(Cross_hairs) then begin
AXIS, 0, 0, YAXIS=0, YTHICK=0.5,\$
YTICKFORMAT='(a1)', \$
YTICKS=ntick_val-1, YTICKV=tick_val, \$
YMINOR=5, YSTYLE=1, \$
YRANGE=[-max_var,max_var], \$
_EXTRA=E_Cross

AXIS, 0, 0, YAXIS=1, YTHICK=0.5,\$
YTICKFORMAT='(a1)', \$
YTICKS=ntick_val-1, YTICKV=tick_val, \$
YMINOR=5, YSTYLE=1, \$
YRANGE=[-max_var,max_var], \$
_EXTRA=E_Cross

AXIS, 0, 0, XAXIS=0, XTHICK=0.5,\$
XTICKFORMAT='(a1)', \$
XTICKS=ntick_val-1, XTICKV=tick_val, \$
XMINOR=5, XSTYLE=1, \$
XRANGE=[-max_var,max_var], \$
_EXTRA=E_Cross

AXIS, 0, 0, XAXIS=1, XTHICK=0.5,\$
XTICKFORMAT='(a1)', \$
XTICKS=ntick_val-1, XTICKV=tick_val, \$
XMINOR=5, XSTYLE=1, \$
XRANGE=[-max_var,max_var], \$
_EXTRA=E_Cross

; If NO radial lines are requested but crosshairs are, then
; label the angles for the four cross hairs. Just set

endif

; ++++
; Plot the radial lines - this is now a more general way of doing this,
; allows for easier changing of the number of radial lines. PBK 22 Feb.
2000.
;
; First check if actually want to plot radial lines.
;
; Check to see if the N_RAD_LINES keyword is set, if so
; default to 12. NO NEED TO DO THIS NOW< IS DONE ABOVE - PBK 25 July 2000.
;
rl_spac=360./float(n_rl)
a=replicate(max_var,n_rl)
; Calculate plot angles, l_theta.
l_theta=!DTOR*(450.-rl_spac*findgen(n_rl))
; Plot the lines
for i=0,n_rl-1 do OPLOT, [0,a[i]], [0,l_theta[i]], /polar, \$
;
; ++++
; Annotate the radial lines - this is now a more general way of doing this,
; allows for easier changing of the number of radial lines. PBK 22 Feb.
2000.
; Setup the radius values for the labels, including an offset.

; Turn the radius and angle into rectangular coords. Use theta from above.
polar_coord=make_array(2,n_rl,/float)
result=cv_coord(from_polar=polar_coord,/to_rect)

; Set the label string here. Now more general.
label=strcompress(fix(rl_spac*indgen(n_rl)),/remove_all)

; Assign the rectangular coords. Adjust y value.
l_x=result[0,*]
l_y=result[1,*]-max_var*0.02

; Set the alignment for xyouts.
l_align=replicate(0.5,n_rl)
al0=where(th_hold gt 0 and th_hold lt 180, cnt_al0)
if cnt_al0 gt 0 then l_align[al0]=0.
al1=where(th_hold lt 360 and th_hold gt 180, cnt_al1)
if cnt_al1 gt 0 then l_align[al1]=1.

; ++++
; Try to set the character size for the labels according
; to the number of plots per page.
nplots=float(total(!p.multi[1:2]))
case 1 of
nplots le 2.: chsize=1.
nplots gt 2. and nplots le 4.: chsize=0.65
nplots gt 4. and nplots le 6.: chsize=0.5
nplots gt 6.: chsize=0.4
endcase

; Plot the labels.
for i=0,n_rl-1 do XYOUTS, l_x[i], l_y[i], label[i], /data, \$
CHARSIZE=chsize, align=l_align[i]
;
; endif for plotting radial lines
endif
; ++++
; Plot two radial lines indicating the baseline sector
; if the BASELINE keyword is set. Can also be used to
; plot any type of indicating line.
if KEYWORD_SET(Baseline) then begin
; If just the baseline keyword is set then plot the
; defaut baseline indicators for Cape Grim, ie 190
; and 280 degrees. If there are two elements or more
; in Baseline keyword, they are angles so plot lines
; along these angles.
if n_elements(Baseline) eq 1 then b_theta=[190.,280.] \$
else b_theta=baseline

; Set the radius to plot out to and convert angles to
b_theta=(450.-b_theta)*!DTOR

; Plot the lines here.
for i=0,n_elements(b_theta)-1 do \$
endif

; ++++
; Plot the actual data here! Turn the winddirn into
; radians and offset it appropriately for plotting
; in compass coordinate system (North-Sount, East-West).
linestyle = keyword_set(linestyle) ? linestyle : 0
thick = keyword_set(thick) ? thick : 4
color = keyword_set(color) ? color : 0

OPLOT, var, (450.-wind_dirn)*!DTOR, /polar, \$
linestyle=linestyle, thick=thick, color=color
;
if n_elements(var2) gt 0 then begin
v2_linestyle = keyword_set(v2_linestyle) ? v2_linestyle : 1
v2_thick = keyword_set(v2_thick) ? v2_thick : 4
v2_color = keyword_set(v2_color) ? v2_color : 0
OPLOT, var2, (450.-wind_dirn2)*!DTOR, /polar, \$
linestyle=v2_linestyle, thick=v2_thick, color=v2_color
endif
;
; ++++

END

```