Re: plane rotation

[Ansley Manke <ansley.b.manke@xxxxxxxx>]

Hi all,
Part II.  The rotated output grid is really a rectangular grid, it just 
does not follow constant lines of longitude and latitude.  So we can put 
the rotated data onto an x-y grid, and not have to deal with it in 
curvilinear coordinates.  Think of the equator.  In the lat-lon grid 
it's a set of points with varying longitude and constant latitude.  In 
the rotated grid both the longitude and latitude coordinates vary.  We 
can use the SAMPLEXY function to get the values of the data field at 
that set of longitudes and latitudes. We can do this for all the points 
on the 2-D grid.  Then create the output, rotated grid in terms of index 
space, I and J (abstract X and Y axes).  In out map projection script, 
we have units of radians, so I'll do the samplexy on the data in radians.

First, the rotate script needs the parameters central_meridian and 
standard_parallel so it will work for data on any part of the globe.

  
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
  \ cancel mode verify  
  ! mp_rotate.jnl --  Sets up variables for a Rotated projection using
  !                  'curvilinear coordinates' code in Ferret v4.50
  !                  
  ! Ansley Manke 4/28/05
  ! based on other map projection scripts
 
  ! Description:  Sets up variables for a Rotated map of the world
  !
  ! Usage:    
  !  go mp_rotate angle [central meridian] [standard_parallel]
  !
  ! arg 1 - angle to rotate from horizontal X axis, in degrees
  ! arg 2 - longitude used for the center of the projection, 
mp_central_meridian
  ! arg 3 - latitude used for the center of the projection, 
mp_standard_parallel
 
  ! Example:
  !  use coads_climatology
  !  go mp_rotate 30
  !  set grid sst
  !  shade/noaxis sst[l=1], x_page, y_page
  !
  let/quiet mp_x = x
  let/quiet mp_central_meridian = $2%(mp_x[i=@max] + mp_x[i=@min])/2%
  let/quiet mp_y = y
  let/quiet mp_standard_parallel = $3%(mp_y[j=@max] + mp_y[j=@min])/2 
%  ! ???
 
  let/quiet Pi = 3.14159265
  let/quiet deg2rad = Pi / 180.0
 
  let/quiet mp_R = $1* deg2rad
 
  let/quiet mp_lambda0 = mp_central_meridian * deg2rad
 
  let/quiet mp_phi0 = mp_standard_parallel * deg2rad
 
  let/quiet mp_lambda = mp_x * deg2rad - mp_lambda0
  let/quiet mp_phi = mp_y * deg2rad - mp_phi0
 
  let/quiet x_page = mp_lambda* cos(mp_R) - mp_phi* sin(mp_R)
  let/quiet y_page = mp_phi* cos(mp_R) + mp_lambda* sin(mp_R)
 
  let/quiet mp_mask = 1
  set mode/last verify
  
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!

 
A sample script, using a subset of the ETOPO data so we can see the 
rotation easily.

  yes? use etopo20
  yes? save/clobber/file=esubset.nc  rose[x=100:150,y=-40:0]
  yes? can data/all

     ! Do the rotated map projection
  yes? use esubset
  yes? set grid rose
  yes? go mp_rotate 30
  yes? shade/pal=land_sea rose, x_page, y_page

     ! The rotation variables are defined in radians.
     ! Define axes in radians, and transfer ROSE onto them

  yes? def axis/x xrad=mp_lambda
  yes? def axis/y yrad=mp_phi
  yes? let roserad = rose[gx=xrad@asn,gy=yrad@asn]

     ! Set up the points to sample ROSERAD at;
     ! all the locations on the rotated grid, as a
     ! simple lists of x, y, and function.

  yes? let xx = xsequence(x_page)
  yes? let yy = xsequence(y_page)
  yes? let rpts = samplexy(roserad, xx, yy)

     ! This is the output grid, in rotated space. (Give
     ! the axes units of INDEX to emphasize this is in
     ! index space).

  yes? let nx = `rose,return=isize`
  yes? let ny = `rose,return=jsize`
  yes? def axis/x=1:`nx`:1/units=index xrot
  yes? def axis/y=1:`ny`:1/units=index yrot

     ! The reshape function puts our sampled points on
     ! the output grid.

  yes? let outvar = x[gx=xrot] + y[gy=yrot]
  yes? shade/pal=land_sea reshape(rpts, outvar)