Re: Another Map Projection Problem with Map_Proj_Init [message #78617 is a reply to message #78416] |
Mon, 05 December 2011 16:43  |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
David Fanning writes:
> I have discovered another map projection problem this morning.
> It seems Map_Proj_Init is calculating incorrect UV ranges
> (they are wrong by 10 degrees) when the center latitude is
> not 0.0. Here is a short test program:
>
> FOR j=0,10 DO BEGIN
> center_lon = 0 > j*36.0d < 360.0
> mapStruct = Map_Proj_Init(117, CENTER_LONGITUDE=center_lon)
> print, center_lon, mapStruct.uv_box[[0,2]]
>
> ll = Map_Proj_Inverse(mapStruct.uv_box[[0,2]], $
> mapStruct.uv_box[[1,3]], MAP_STRUCTURE=mapStruct)
> print, center_lon, ll[0,0], ll[0,1]
> print, ""
> ENDFOR
> END
>
> Here are the results. The numbers should be symmetric,
> but you can see they are not.
>
> 0.000000 -20015077. 20015077.
> 0.000000 180.00000 -180.00000
>
> 36.0000 -19570298. 19347908.
> 36.0000 -140.00000 -150.00000
>
> 72.0000 -19125518. 19792688.
> 72.0000 -100.00000 -110.00000
>
> 108.000 -19792688. 19125518.
> 108.000 -70.000000 -80.000000
>
> 144.000 -19347908. 19570298.
> 144.000 -30.000000 -40.000000
>
> 180.000 -20015077. 18903129.
> 180.000 0.00000000 -10.000000
>
> 216.000 -19570298. 19347908.
> 216.000 40.000000 30.000000
>
> 252.000 -19125518. 19792688.
> 252.000 80.000000 70.000000
>
> 288.000 -19792688. 19125518.
> 288.000 110.00000 100.00000
>
> 324.000 -19347908. 19570298.
> 324.000 150.00000 140.00000
>
> 360.000 -20015077. 20015066.
> 360.000 180.00000 179.99990
>
> I couldn't figure out why my grid lines were always getting
> screwed up when I calculated them using the UV range! :-(
It turns out, according to the folks in technical support,
that to get correct UV_BOX values from Map_Proj_Init you
MUST use the LIMIT keyword, even (well, *especially*) if
you don't plan to limit the map projection in any way.
Another way to say this is that if you plan to use
Map_Proj_Init and set the CENTER_LONGITUDE to anything
other than 0, you MUST also use the LIMIT keyword to
produce correct results. This is, as you might expect,
not documented. :-)
I've updated my cgMap object to handle the problem
correctly.
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
|