SET MODE VERIFY

! poly_vec_demo.jnl 7/03 *acm*
! Description: demo of filled-polygon vectors
! Uses poly_vectors.jnl and mp_poly_vectors.jnl scripts  
! contributed by E. D. Cokelet, NOAA/PMEL, 2003

! The scripts poly_vectors.jnl and mp_poly_vectors.jnl compute polygon
! coordinates for drawing vectors in Cartesian or curvlilinear coordinates 
! as polygon arrows or sticks, optionally filled with color.  These scripts 
! are similar to scattered_vectors.jnl in that they set up the data, which 
! is then plotted with a separate command, in this case the POLYGON command. 
! A third script, poly_arrow_key.jnl, makes an arrow key at any location on 
! the page.
  
! First lets make a plot in Cartesian coordinates.  We will call the script 
! poly_vectors.jnl to set up the polygons for plotting.

PAUSE
SET WINDOW/siz=0.6
CANCEL MODE logo

USE coads_climatology
SET REGION/X=120:180/Y=20:60/L=1

! poly_vectors.jnl must have the plot parameters (such as axis lengths) 
! already set.  Draw a 2-D plot over the region of interest.  

GO basemap x=120:180 y=20:60 20
LABEL/NOUSER `($ppl$xlen)/2`, `($ppl$yorigin)-1`, 0, 0, 0.16, "UWND, VWND"


! The inputs to poly_vectors.jnl are the components of the vectors to be 
! plotted and the X and Y coordinates of the vectors.  We need to limit 
! the region of the data explicitly, so that the calculations are done 
! only for the vectors we are interested in.

PAUSE
LET x_vec = X[GX=uwnd,X=120:180] + 0*Y[GY=uwnd,Y=20:60]
LET y_vec = 0*X[GX=uwnd,X=120:180] + Y[GY=uwnd,Y=20:60]
LET u_vec = uwnd[X=120:180,Y=20:60,L=1]
LET v_vec = vwnd[X=120:180,Y=20:60,L=1]

GO poly_vectors.jnl x_vec y_vec u_vec v_vec 5 "arrow"

! Now the polygon coordinates have been defined.  Next we issue a POLYGON
! command to draw the arrow polygons. We can use a third argument on the 
! POLYGON command to fill the polygons with colors according to any variable. 
! Here we make a list of the values of sea level pressure which corresponds 
! to the original vector components. The color key is for the fill color. 
! We use the /LINE qualifier to outline the arrows with black.  Call 
! poly_arrow_key.jnl to make the vector key.

PAUSE

LET press = YSEQUENCE(slp[X=120:180,Y=20:60,L=1])
POLYGON/OVER/NOAX/NOLAB/KEY/LINE/PAL=rainbow plt_x_arrow, plt_y_arrow, press

! Make a vector arrow key.

GO poly_arrow_key "winds", "5 m/s"

PAUSE

CAN VAR/ALL
CAN REGION

! For a larger-scale plot we need to sub-sample the vectors before
! sending them as arguments to the poly_vectors script.  This is done
! automatically by the Ferret VECTOR command, but here the user needs
! to do the sub-sampling.

PAUSE

USE coads_climatology
SET REGION/L=1/X=0:360/Y=-80:80

! Make a 2-D plot to set plot parameters.

GO basemap x=0:360 y=-80:80 20
LABEL/NOUSER `($ppl$xlen)/2`, `($ppl$yorigin)-1`, 0, 0, 0.16, "Wind Vectors Colored by SLP"

! Subsample the vector data and set up the vectors to be plotted as polygons
PAUSE

LET x_vec = X[GX=uwnd,X=0:360:10] + 0*Y[GY=uwnd,Y=-80:80:5]
LET y_vec = 0*X[GX=uwnd,X=0:360:10] + Y[GY=uwnd,Y=-80:80:5]
LET u_vec = uwnd[X=0:360:10,Y=-80:80:5,L=5]
LET v_vec = vwnd[X=0:360:10,Y=-80:80:5,L=5]

GO poly_vectors.jnl x_vec y_vec u_vec v_vec 8 "arrow"

PAUSE

! Draw the arrow polygons, colored by SLP, and make a key for the arrows.
! Note the different plot appearance when we leave off the LINE qualifier
! on the POLYGON command.

LET press = YSEQUENCE(slp[X=0:360:10,Y=-80:80:5,L=5])
POLYGON/OVER/NOAX/NOLAB/KEY/PAL=rainbow/LEV=(980,1030,5) plt_x_arrow, plt_y_arrow, press
GO poly_arrow_key "", "8 m/s", red


CAN REGION
CAN VAR/ALL

PAUSE
! Next we will do some polygon-vector plots using map projections. 

! The steps are: Define a subset of the data to plot, define a map projection, 
! and draw some plot in the map projection, using the curvilinear form 
! of the plotting command. Define variables with with U, V, X, and Y for the 
! vectors, call mp_poly_vectors to set up the polygon coordinates and use a 
! curvlinear polygon command to plot them.

USE coads_climatology

SET REGION/X=131:181/Y=21:61/L=5
LET mp_central_meridian = (131+181)/2
LET mp_standard_parallel = (60+20)/2

GO mp_stereographic_north  `mp_central_meridian` `mp_standard_parallel`
SET GRID uwnd
GO mp_aspect

! Make a FILL plot using the map projection to set the plot parameters for 
! the mp_poly_vectors script

FILL/NOAXIS/NOKEY/PAL=inverse_grayscale/TITLE="Winds Colored by Pressure" slp, x_page,y_page
GO mp_fland 20 black

PAUSE

! Call mp_poly_vectors to set up the vectors to be plotted as polygons

LET x_vec = X[GX=uwnd,X=131:181] + 0*Y[GY=uwnd,Y=21:61]
LET y_vec = 0*X[GX=uwnd,X=131:181] + Y[GY=uwnd,Y=21:61]
LET u_vec = uwnd[X=131:181,Y=21:61,L=5]
LET v_vec = vwnd[X=131:181,Y=21:61,L=5]

GO mp_poly_vectors x_vec y_vec u_vec v_vec 3 "arrow"


! Draw the polygons, colored by SLP

PAUSE

LET press = YSEQUENCE(slp[X=131:181,Y=21:61,L=5])
POLYGON/OVER/NOAX/NOLAB/KEY/LINE/PAL=rainbow mp_x_arrow, mp_y_arrow, press

GO mp_graticule 131 181 12.5 21 61 10

! Now we illustrate positioning the vector arrow key. By default it is 
! positioned at the lower left.  We can choose to center it horizontally
! or vertically by sending C as the argument for a position.

PAUSE 

GO poly_arrow_key "default key", "3 m/s", red, line
GO poly_arrow_key "left top key", 3, blue, " ", 2, `($ppl$height)-1.5`
GO poly_arrow_key "centered bottom", 3, violet, " ", c, 0.2
GO poly_arrow_key "upper right", 3, green, line, `($ppl$width)-1.6`, `($ppl$height)-1.5`,0.14

CAN DATA/ALL
CAN VAR/ALL
CAN REGION

! Next we'll make a view from space.  We need the mp_mask for this projection, 
! to keep the vectors on other side of the earth from "showing through" the
! map.  Apply it to the polygon data as it is plotted. This time we will 
! use explicit regridding to put the vector components onto a coarser grid.

PAUSE
USE coads_climatology
SET REGION/L=1

! Set up the vectors to be plotted as polygons

DEFINE AXIS/X=0:360:8 xax
DEFINE AXIS/Y=-80:80:4 yax
DEFINE GRID/X=xax/Y=yax gxy

LET x_vec = X[GX=xax] + 0*Y[GY=yax]
LET y_vec = 0*X[GX=xax] + Y[GY=yax]
LET u_vec = uwnd[g=gxy,L=1]
LET v_vec = vwnd[g=gxy,L=1]

! Make a plot with the map projection to set the plot parameters. (FILL plot
! with the palette set to all white, for a blank underlying plot)
PAUSE

GO mp_orthographic 210 45
SET GRID gxy
GO mp_aspect

FILL/NOAXIS/NOKEY/NOLAB/PAL=white/TITLE="Winds Colored by Pressure" u_vec, x_page,y_page
GO mp_fland
GO mp_graticule 0 360 30 -90 90 20

GO mp_poly_vectors x_vec y_vec u_vec v_vec 5 "arrow"

! draw the arrow polygons, colored by SLP, and the vector key.
PAUSE

LET press = YSEQUENCE(slp[g=gxy,L=1])
POLYGON/OVER/NOAX/NOLAB/LINE/KEY/PAL=rainbow mp_mask*mp_x_arrow, mp_mask*mp_y_arrow, mp_mask*press

! Vector key

GO poly_arrow_key " ", "5 m/s", red, " ", 1, `($ppl$yorg)`


CAN DATA/ALL
CAN VAR/ALL
CAN REGION

! Finally, the poly_vector scripts also have the option to plot the vectors 
! as stick vectors, rather than arrows.  All the commands stay the same,
! but with the final argument to poly_vectors.jnl or mp_poly_vectors.jnl
! being "stick" rather than "arrow".  Here is the first plot from this demo
! made with stick vectors:

PAUSE
USE coads_climatology
SET REGION/X=120:180/Y=20:60/L=1

GO basemap X=120:180 Y=20:60 20
LABEL/NOUSER `($ppl$xlen)/2`, `($ppl$yorigin)-1`, 0, 0, 0.16, "Wind Vectors Colored by SLP"

LET x_vec = X[GX=uwnd,X=120:180] + 0*Y[GY=uwnd,Y=20:60]
LET y_vec = 0*X[GX=uwnd,X=120:180] + Y[GY=uwnd,Y=20:60]
LET u_vec = uwnd[X=120:180,Y=20:60,L=1]
LET v_vec = vwnd[X=120:180,Y=20:60,L=1]

GO poly_vectors.jnl x_vec y_vec u_vec v_vec 5 "stick"

LET press = YSEQUENCE(slp[X=120:180,Y=20:60,L=1])
POLYGON/OVER/NOAX/NOLAB/KEY/PAL=rainbow plt_x_arrow, plt_y_arrow, press

GO poly_arrow_key "winds", "5 m/s" red