INTERNAL/REVIEW: An example of creating a 3D polar plot

Topic:

The following routine shows an example of how to plot a 3-D dataset given in spherical coordinates using IDL's object graphics system. The example uses XOBJVIEW so that it is easy to rotate and zoom.Discussion:
The example routine "polar3d.pro" takes three arguments. These are Phi, Theta and Rho (or longitute, latitude and radius) where Phi and Theta are assumed to be in radians. The routine uses the QHULL (new in IDL 5.5) to find the connectivity between the specified points.

The three arguments must be arrays of equal sizes. If no arguments are specified, then a default doughnut shaped dataset is used.

This routine works for visualizing data that has only one distinct radius value for each lon/lat pair. An example of such a dataset would be an antenna radiation pattern.

The example also demonstrates how to expand the functionality of XOBJVIEW.
Solution:

;+
; NAME:
; POLAR3D
;
; PURPOSE:
; Create an interactive 3-D polar surface plot.
;
;
; CATEGORY:
; Visual
;
; CALLING SYNTAX:
; POLAR3D, Phi, Theta, Rho
;
; INPUT ARGUMENTS:
; Phi = Longitude.  Angle from the positive X-axis to the vector from the
;       origin to the projection of each point onto the X-Y plane.
;
; Theta = Latitude. Angle between the X-Y plane and the vector from the
;         origin to each point.
;
; Rho = Radius.  Distance from the origin to each point.
;
; OUTPUTS:
; No outputs.
;
; SIDE EFFECTS:
; None.
;
; RESTRICTIONS:
; Works for datasets that have a unique Rho for each unique Phi/Theta pair.
;
; PROCEDURE:
; Phi, Theta and Rho must all have the same number of elements describing
; the points.  All the (Phi, Theta) pairs are first projected onto the unit sphere.
; QHULL is then used to obtain the convex hull of the points on the unit sphere.
; The connectivity list for the convex hull is finally used to connect the original
; points.
;
; MODIFICATION HISTORY:
; -
pro polar3D_buttons, event
  ;
  ; Event handler for the buttons added to XOBJVIEW
  widget_control, event.id, GET_UVALUE=state
  widget_control, event.id, GET_VALUE=name
  if (strlowcase(name) eq 'gouraud shading') then begin
    state->getproperty, SHADING=shading
    state->setproperty, SHADING=1-shading  ; Toggle SHADING property
    xobjview, REFRESH=event.top
  endif else begin
    state->getproperty, HIDE=hide
    state->setproperty, HIDE=1-hide   ; Toggle HIDE property
    xobjview, REFRESH=event.top
  endelse
end
 
pro polar3D_loadct, event
  ;
  ; Event handler for colortable pulldown list
  widget_control, event.id, GET_UVALUE=obj
  obj[0]->loadct, event.index
  obj[0]->getproperty, RED=r, GREEN=g,BLUE=b
  rgb=transpose([[r],[g],[b]])
  obj[1]->getproperty, DATA=data
  rho=sqrt(data[0,*]^2+data[1,*]^2+data[2,*]^2)
  obj[1]->setproperty, VERT_COLORS=rgb[*,bytscl(rho)]
  xobjview, REFRESH=event.top
end
 
pro polar3D_cleanup, base
  ; clean up the objects
  widget_control, base, GET_UVALUE=obj
  obj_destroy, obj
end
 
pro polar3D, Phi, Theta, Rho
  if (n_elements(Phi) eq 0) then begin
    N=1000  ; make up some fake data if no data is given  (1000 points)
    theta=acos(randomu(1,N)*2-1)-!pi/2
    phi=2*!pi*randomu(2,N)-!pi
    rho=cos(theta)
  endif
  N=n_elements(phi)
  ; find the points on the unit sphere used for triangulation
  xyz=cv_coord(FROM_SPHERE=transpose([[phi],[theta],[dblarr(N)+1]]),/TO_RECT)
  qhull,xyz,tr  ; find the convex hull (which is the same as the spherical triangulation)
  ; find the connectivity array
  ntr=n_elements(tr)
  conn=reform([replicate(3,1,ntr/3),temporary(tr)],ntr*4/3,/OVERWRITE)
  ; convert actual data to rectangular
  xyz=cv_coord(FROM_SPHERE=transpose([[phi],[theta],[rho]]),/TO_RECT)
  ; create objects
  pal=obj_new('IDLgrPalette')
  pal->loadct,13
  pal->getproperty, red=r, green=g,blue=b
  rgb=transpose([[r],[g],[b]])
  blob=obj_new('IDLgrpolygon',xyz,POLYGONS=conn,VERT_COLORS=rgb[*,bytscl(rho)],/SHADING,DEPTH_OFFSET=1)
  grid=obj_new('IDLgrpolygon',xyz,POLYGONS=conn,STYLE=1,DEPTH_OFFSET=0)
  ax=objarr(3)
  cr=dblarr(2,3)
  for i=0,2 do begin
    ; create axes
    ax[i]=obj_new('IDLgrAxis',i,RANGE=[min(xyz[i,*]),max(xyz[i,*])])
    ax[i]->getproperty, CRANGE=c
    cr[*,i]=c
  endfor
  for i=0,2 do begin
    ; set character sizes and location for axes intersection
    ax[i]->setproperty, LOCATION=[cr[0,0],cr[0,1],cr[0,2]]
    ax[i]->getproperty, TICKTEXT=tt
    r0=cr[*,0]
    r1=cr[*,i>1]
    r0=r0[1]-r0[0]
    r1=r1[1]-r1[0]
    tt->setproperty, CHAR_DIM=[r0/20,r1/20]
  endfor
  ; visulize
  dev=strlowcase(!d.name)  ; find the current screen size
  if (dev eq 'x') or (dev eq 'win') or (dev eq 'mac') then device, GET_SCREEN_SIZE=gss else gss=[800,600]
  xobjview,[blob,grid,ax], XSIZE=long(gss[0]*0.8), YSIZE=long(gss[1]*0.8), TLB=tlb, TITLE='Polar3D'
  ;
  ;  Add more buttons and a pulldown list to the XOBJVIEW window
  rowbase=widget_base(tlb,/ROW)
  base=widget_base(rowbase, /ROW, /NONEXCLUSIVE, EVENT_PRO='polar3D_buttons',$
    KILL_NOTIFY='polar3D_cleanup', UVALUE=[blob,grid,ax,pal])
  b1=widget_button(base, VALUE='Show Triangles', UVALUE=grid)
  b2=widget_button(base, VALUE='Show X-Axis', UVALUE=ax[0])
  b3=widget_button(base, VALUE='Show Y-Axis', UVALUE=ax[1])
  b4=widget_button(base, VALUE='Show Z-Axis', UVALUE=ax[2])
  b5=widget_button(base, VALUE='Gouraud shading', UVALUE=blob)
  loadct, GET_NAMES=names  ; fill in the colortable names in the dropdown list
  ct=widget_droplist(rowbase, VALUE=names, EVENT_PRO='polar3D_loadct',UVALUE=[pal,blob])
  widget_control, b1, /SET_BUTTON
  widget_control, b2, /SET_BUTTON
  widget_control, b3, /SET_BUTTON
  widget_control, b4, /SET_BUTTON
  widget_control, b5, /SET_BUTTON
  widget_control, ct, SET_DROPLIST_SELECT=13
  widget_control, tlb, /UPDATE
  widget_control, tlb, GET_UVALUE=uval
end