Re: Map projection

[Pal Jeremy <jpal@xxxxxxxxxxxxxxx>]

  Hi Hyun-Suk!

I have included a script that computes the lat/lon info for a
Lambert conformal conical projection.  It should be the same
Lambert projection used by the MM5.  (I also have one for the
rotated mercator projection for anyone who is interested).

Hopefully this will save you some time.  Let me know if you have
any questions implementing the projection with the 3-argument
shade/fill command.


  - Jeremy


On Sat, 21 Jul 2001, Kang, Hyun-Suk wrote:

> Dear Ferret users,
> 
> I would like to use Ferret for visulalizing MM5 output.
> I use Lambert conformal map projection in MM5. 
> 
> Please let me know if you have any experience using Lambert conformal map projection
> (not equal lat/lon grids but equal distance grids).
> 
> Thanks.
> 
> Hyun-Suk
> ===========================================================
> Kang, Hyun-Suk
> 
> Atmospheric Sciences Program
> School of Earth and Environmental Sciences
> Seoul National Univerisity, Seoul 151-742, KOREA
> 
> Email: hskang@nwp2.snu.ac.kr
> Phone: 82-2-880-6760
> Fax: 82-2-885-4448
> 
> 
let nxx=$1
let nyy=$2
let delx=$3
let clatt=$4
let clonn=$5
let lattr=$6

let xc=nxx*0.5
let yc=nyy*0.5
let my_x = delx*(i[i=1:`nxx`]-xc)
let my_y = delx*(j[j=1:`nyy`]-yc)


let Pi = 3.14159265
let pifac=atan(1.)/45.            ! CONVERT DEGREES TO RADIANS
let psi1=lattr*pifac
let colat=90.-clatt                 ! COLAT OF CENTER
let psi0=colat*pifac
let rad2deg = 180.0 / Pi
let earthr = 6.37e6   ! Earth Radius
let grdfac = 0.716

let rgrdf=1./grdfac
let pi90=2.*atan(1.)              ! PI OVER TWO (90 DEGREES)
let ans=earthr*sin(psi1)*rgrdf
let pf=tan(psi1*0.5)/(-1*ans^rgrdf)
let c1=-1*clonn-90.*rgrdf
let arg1=tan(psi0*0.5)/tan(psi1*0.5)
let c2=ans*(arg1)^grdfac
! let c2=ans*(tan(psi0*0.5)/tan(psi1*0.5))^grdfac

let xlp=atan2(my_y-c2,my_x)
let r=abs(my_x/cos(xlp))
let psi=abs(2.*atan(pf*(r^rgrdf)))        ! CO LATITUDE (RADIANS)
let pfi=pi90-psi

let my_new_lat = pfi/pifac
let my_new_lon = xlp*rgrdf/pifac-c1