Dear idl users,
Please prevent me from going crazy:
I want to plot output from an atmospheric model. This output comes as a
regular grid in a ROTATED lat-lon projection. I want to plot these data
as an image, either in a window or as postscript. On top of this image I
want to plot map contours. I am aiming at showing this in the native
projection, so I would not like to bend the image with MAP_IMAGE, if
this is possible, (in contrast to a solution by David Fanning to a
similar problem some time ago). I have been tweaking !d.position,
!w.window, keywords to TV etc. etc. and it is driving me crazy. When I
have found something that works for one area, it does not work when I
change the area, i.e. number of points, aspect ratio, position on the
globe, and the like.
I am using IRIX 6.3 on an SGI O2 and IDL 5.1.
Any suggestions are welcome! My setup is something like the following:
x,y,z are fltarr(nlon,nlat)
x1=x(0,0)+0.5*(x(1,1)-x(0,0))
y1=y(0,0)+0.5*(y(1,1)-y(0,0))
x2=x(0,nlat-1)+0.5*(x(1,nlat-2)-x(0,nlat-1))
y2=y(0,nlat-1)+0.5*(y(1,nlat-2)-y(0,nlat-1))
x3=x(nlon-1,nlat-1)+0.5*(x(nlon-2,nlat-2)-x(nlon-1,nlat-1))
y3=y(nlon-1,nlat-1)+0.5*(y(nlon-2,nlat-2)-y(nlon-1,nlat-1))
x4=x(nlon-1,0)+0.5*(x(nlon-2,1)-x(nlon-1,0))
y4=y(nlon-1,0)+0.5*(y(nlon-2,1)-y(nlon-1,0))
; These are supposedly corner points for the whole image
if ps ne 0 then begin
device,/color,bits=8
device,xsize=18.0,ysize=(18.0*nlat)/nlon
device,xoffset=0.5,yoffset=1.5
endif else begin
; Display is a window:
keyscale=1.00
set_plot,'X'
device,pseudo_color=8
window,0,xsize=740,ysize=round((740L*nlat)/nlon),retain=2,/f ree,colors=256
; More playing around...
px = [0.05,0.85]*!d.x_vsize
py = [0.05,0.85]*!d.y_vsize
sx=px(1)-px(0)+1
sx=round(sx*keyscale)
px(1) = sx - 1 + px(0)
sy=py(1)-py(0)+1
!p.position=[0.05,0.05,0.85,0.85]
endelse
; More tweaking of !p.position etc...
!x.window=[0,1]
!y.window=[0,1]
if polat gt -89.9 then begin ; Rotated lat-lon
map_set,0,polon-90,/cylin,/noborder,/advance, $
limit=[y1,x1,y2,x2,y3,x3,y4,x4],$
CENTRAL_AZIMUTH = -polat-90.0
endif else begin ; Ordinary lat-lon
map_set,/cylin,/noborder,/advance, $
limit=[y4,x1,y2,x3]
endelse
zmod=bytarr(nlon,nlat)
if ps eq 0 then begin
zmod=congrid(zmod,sx,sy)
tv,zmod,px(0)-0.6*(px(1)-px(0))/nlon,py(0)
endif else begin
tv,zmod,-0.005,0.0,xsize=1.005,ysize=1.003,/norm
; Right, more tweaking
endelse
map_grid
map_continents
--
Thanks in advance for any help! /Ole
--
Ole B�ssing Christensen, Ph.D. E-mail: bossing@dmi.dk
Climate Research Div. Phone : (+45) 39 15 74 26
Danish Meteorological Institute Fax: (+45) 39 15 74 60
Lyngbyvej 100, DK-2100 Copenhagen �, Denmark
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
