Re: [ferret_users] Robinson map projection

[Marco Steinacher <steinacher@xxxxxxxxxxxxxxxx>]

Hi Patrick,

On 31.01.2013 09:10, Patrick Brockmann wrote:
> Following http://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2011/msg00129.html
> I am wondering if anyone has succeeded to code a Robinson projection with ferret ?
> It is not so easy because this projection is achieved through the use of tabular coordinates rather than mathematical formulae.
> Any news on this topic are welcome.

Savric et al. (2011) [1] proposed a polynomial approximation of the
Robinson projection. The attached Ferret script by Laurent Bopp is a
modification of the mp_winkel_i script and makes use of this
approximation for "Robinson-like" projections. Maybe that's good enough
for your purpose.

Cheers,
Marco

[1]
http://cartography.oregonstate.edu/pdf/2011_Savric_etal_PolynomialNaturalEarth.pdf

-- 
***************************************
Dr. Marco Steinacher

Climate and Environmental Physics
Physics Institute, University of Bern
Sidlerstrasse 5, CH-3012 Bern

Phone ++41 (0)31 631 34 02
Fax   ++41 (0)31 631 87 42
steinacher@xxxxxxxxxxxxxxxx
http://www.climate.unibe.ch/~steinach
OpenPGP Key ID: 0x62937F7F
***************************************

\cancel mode verify	
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
! mp_winkel_i.jnl --  Sets up variables for a Robison I projection using
!                  'curvilinear coordinates' code in Ferret v4.50
!		    		
!L. Bopp - Robinson projection based on a polynomial Savric et al. 2011
! 2013
! X = R ∙ λ ∙ (A0 + A2 ∙ φ2 + A4 ∙ φ4) (Eq. 2)
!Y = R ∙ (A1 ∙ φ + A3 ∙ φ3 + A5 ∙ φ5)
!where:
! and Y are projected coordinates;
!φ and λ are the latitude and longitude;
!R is the radius of the generating globe;
!A0 = 0.8507;
!A1 = 0.9642;
!A2 = -0.1450;
!A3 = -0.0013;
!A4 = -0.0104; and
!A5 = -0.0129. 5/97

! Description:  Sets up variables for a Winkel I map of the world
!
! Usage:                   arg1               arg2
!  go mp_winkel_1 [central meridian] [standard parallel]
!
! arg 1 - longitude used for the center of the projection
! arg 2 - latitude used for the center of the projection

! Example:
!  use coads_climatology 
!  go mp_robinson_like
!  set grid sst
!  shade sst[l=1], x_page, y_page
!
! Note 1: If you intend to plot an overlay which is a subregion
!         of the original plot you will need to specify the two
!         optional arguments as in:
!
!  use coads_climatology 
!  go mp_robinson_like
!  set grid sst
!  shade sst[l=1], x_page, y_page
!  go mp_robinson_like `mp_central_meridian` `mp_standard_parallel`
!  set region/x=40e:110e/y=60s:20s
!  shade/over slp[l=1], x_page, y_page
!
! Note 2: The x-range for the subregion to be overlayed must
!         be: ( `mp_cent_mer` - 180 < x < `mp_cent_me`r + 180 )

! This journal file comes from equations in the following book:
!
!    An Album of Map Projections
!    U.S. Geological Survey
!    Professional Paper 1453
!
! Equation numbers refer to the equation number in this book.

if `$2%0% gt 90` then
    query/ignore $3"<The standard parallel must be between -90 and 90"
elif `$2%0% lt (-90)` then
    query/ignore $3"<The standard parallel must be between -90 and 90"
endif

let/quiet mp_x = x
let/quiet mp_central_meridian = $1%(mp_x[i=@max] + mp_x[i=@min])/2%
let/quiet mp_y = y
let/quiet mp_standard_parallel = $2%(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
let/quiet mp_lambda = mp_x * deg2rad
let/quiet mp_lambda0 = mp_central_meridian * deg2rad
let/quiet mp_phi = mp_y * deg2rad


! X = R ∙ λ ∙ (A0 + A2 ∙ φ2 + A4 ∙ φ4) (Eq. 2)
!Y = R ∙ (A1 ∙ φ + A3 ∙ φ3 + A5 ∙ φ5)
!where:
! and Y are projected coordinates;
!φ and λ are the latitude and longitude;
!R is the radius of the generating globe;
let/quiet A0 = 0.8507
let/quiet A1 = 0.9642
let/quiet A2 = -0.1450
let/quiet A3 = -0.0013
let/quiet A4 = -0.0104
let/quiet A5 = -0.0129

let/quiet x_page = mp_R * (mp_lambda-mp_lambda0) * (A0 + A2*mp_phi*mp_phi + A4*mp_phi*mp_phi*mp_phi*mp_phi)
let/quiet y_page = mp_R * (A1*mp_phi + A3*mp_phi*mp_phi*mp_phi + A5*mp_phi*mp_phi*mp_phi*mp_phi*mp_phi)

!let x_page = mp_R * (mp_lambda-mp_lambda0) 
!let/quiet y_page = mp_R * mp_phi

let/quiet mp_mask = 1

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
set mode/last verify