Make an FAQ out of recent thread about rotating a grid

[karlmsmith]

Reported by @AnsleyManke on 28 Apr 2005 20:20 UTC
scripts in ~ansley/ans_ferret/users/ansley/reshape_to_rotated_grid.jnl and
~ansley/FERRET/jnls/go/mp_rotate.jnl

Hi all,
After an off-line conversation with Ned Cokelet of our lab, the question comes
down to defining different kinds of axes. They do not need to be rectangular
latitude-longitude axes. Because the rotated grid has different values of
longitude along a grid line, and likewise for latitude, it will need to be
expressed as a "curvilinear" grid.

So, you can define the new grid using cos and sin relationships. If all you
need is graphical output, you can use the ideas from the map projection scripts
(like mp_mercator.jnl etc) to define the coordinates for the rotated grid, and
then use the curvilinear plotting capabilities of Ferret. You need to create two
new variables, one with the longitudes at each grid point of the new 2-D grid,
and one with the latitudes, also in 2 dimensions.

This would be along these lines, if your variables is called VAR.

let xpts = x[gx=var]
let ypts = y[gy=var]
let angle_deg = 30.

let Pi = 3.14159265
let deg2rad = Pi / 180.0

let angle_rad = angle_deg* deg2rad
let xrad = xpts * deg2rad
let yrad = ypts * deg2rad

let xcoord = xrad* cos(angle_rad) - yrad* sin(angle_rad)
let ycoord = yrad* cos(angle_rad) + xrad* sin(angle_rad)

shade/noaxes var, xcoord, ycoord

If you need to do computations with the data on the new grid, you could use the
RECT_TO_CURV regridding function which was released with Ferret version 5.8, and
regrid var onto the curvilinear grid defined by xcoord and ycoord.

(I'll incorporate the above into a new map projection script, mp_rotate.jnl for
future releases.)

Ansley

Hemerson Tonin wrote:

Hi Ferreters,

Has anybody out there a good idea to define axis (and grid) to be inclined ?

I would like to rotate my grid on XY plane, and my "new X-axis" must have an
angle in relation to Equator.

original (XY plane)
I need (XY plane)

                          / <-- "new X -axis"

;

-----------> Equator and "current X - axis"
---------/---------------------------- Equator

Hemerson Tonin
School of Chemistry, Physics and Earth Sciences
Flinders University
Adelaide - Australia

-------- Original Message --------
Subject: Re: plane rotation

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)

Migrated-From: http://dunkel.pmel.noaa.gov/trac/ferret/ticket/1225