[Thread Prev][Thread Next][Index]
Re: Wind rose plot
On Fri, 29 Oct 2004, Ansley Manke wrote:
> Hi Francisco,
> ... if you mean a script that depicts wind data in a circular plot with
> direction and magnitude it should be possible using the POLYGON command.
!-------------------------------------------------------------------
! A compass rose for wind direction only
! step 0 : suppose that you are provided with the frequencies of
! occurance of the 16 wind directions N,NNE,NE,ENE,E, ... NW,NNW
let freq={7,5,2,6,6,3,6,6,7,10,11,9,6,3,6,7}
! step 1 : put them on a suitable axis (units of degrees)
def axis/x=0:337.5:22.5 x16 ; def grid/x=x16 g16
let r16=freq[g=g16,gx=@asn]
! step 2 : extend the grid to complete the circle. maybe there is a
! way to avoid this using modulo but I couldn't see it
def axis/x=0:360:22.5 x17 ; def grid/x=x17 g17
let r17=missing(r16[g=g17],r16[i=1])
! step 3 : a 22.5 degree increment is too coarse - move to 1 degree
def axis/x=0:360:1 x360 ; def grid/x=x360 g360
let r360=r17[g=g360,gx=@nrst] ! @nrst is nearest neighbor regridding
! step 4 : represent the frequency distribution in polar coords
let d2r=atan(1.0)/45 ! degree to radian factor
let x360=r360*sin(d2r*x[g=g360])
let y360=r360*cos(d2r*x[g=g360])
! step 5 : set axis lengths for square figure and shade polygon
ppl axlen,6,6
polygon/nolab/pal=red/hlim=-10:10/vlim=-10:10 x360,y360
plot/o/vs/nolab/sym=4 0,0
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! wind rose when frequencies are provided for speed and direction bins
! Continue with a similar example (16 directions) but with say 8
! speed classes centered at 1,2,... m/s. With the data on a grid the
! 3-argument "shade" becomes available, and modulo becomes easy
def axis/x=0:337.5:22.5/modulo x16 ! 16 directions N,NNE,... as before
! use edges and from_data in defining the "speed" axis
let spd=ysequence({0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5})
def axis/y/name=y8/edges/from_data spd
def grid/x=x16/y=y8 gxy
! step 0 : make a fake dataset
let xy=x[g=gxy]+y[g=gxy]
let v=randu(xy)
! step 1 : introduce the polar coords of the grid
let d2r=atan(1.0)/45
let xv=y[g=xy]*sin(d2r*x[g=xy])
let yv=y[g=xy]*cos(d2r*x[g=xy])
! step 2 : set axis lengths for square figure and shade
ppl axlen,6,6
shade/nolab/hlim=-10:10/vlim=-10:10 v,xv,yv
! this simple-minded approach has a few problems
! a) there is no provision for the zero speed class
! b) the azimuth resolution is coarse
! c) only half of the first and last speed bins are drawn
! the first two are easy to solve (and the first might even be
! absent if the input data included the zero class). Completing
! the outermost ring, seems like it should be easier than my
! ad-hoc solution below
! step 3 : improvements
def axis/x=0:360:1/modulo x360 ! for spreading the azimuth
let vzero=0.5 ! the zero speed class, when not included in the grid
def axis/y=0:9:1 y0 ! a y-axis that includes zero speed
def grid/x=x360/y=y0 gxy0
let v0=missing(v[g=gxy0,gx=@nrst],vzero)
let xv0=y[g=gxy0]*sin(d2r*x[g=gxy0])
let yv0=y[g=gxy0]*cos(d2r*x[g=gxy0])
! shade/nolab/hlim=-10:10/vlim=-10:10 v0,xv0,yv0 ! outer ring incomplete
let v0f=if(y[g=gxy0] le 8)then v0
shade/nolab/hlim=-10:10/vlim=-10:10 v0f,xv0,yv0 ! outer ring complete
!--------------------------------------------------------------------
Hope these help,
Mick
[Thread Prev][Thread Next][Index]
Dept of Commerce /
NOAA /
OAR /
PMEL /
TMAP
Contact Us | Privacy Policy | Disclaimer | Accessibility Statement