Hi, here are the commands. If you want to see the file, I can put
those on my public_html or email directly fo you. Note, there are
some routines in here that will be foreign to anyone unfamiliar with
my idl library (winset, colorbar_w), but they do not relate to my
question. The lat and lon coordinates are called and the subroutine
is after the main program.
Thanks,
Howie
;Sample to read binary, unformatted, sequential access, starting
;with an 80-byte title, followed by a REAL*4 array(72,46) of data;
field=fltarr(72,46)
title=strarr(80)
OPENR, 10, 'CO_IIASA_DOM', /F77_UNFORMATTED
READU, 10, title, field
CLOSE, 10
help,field
lat=fltarr(46)
lon=fltarr(72)
;Now call module getcoord to get GISS 4x5 lats and lons
getcoord,lat,lon
;May need to change lon to 0 to 360?
;lon=lon+180.
winset ;set printing device
loadct,39
ncolors=10
print,'number of levels used is: ',ncolors+1
;levels=[0,10,20,30,40,50,60,80,100,120,150,200,500,1000,200 0,50000]
colors=60+findgen(ncolors)*(254-60)/(ncolors-1)
;levels=min(field)+findgen(16)*(max(field)-min(field)/ncolor s)
levels=min(field)+findgen(ncolors+1)*((200000.-0.)/ncolors)
map_set,0,0,title='CO_IIASA_DOM!C',$
position=[0.15,0.35,0.85,0.65]
contour,field,lon,lat,c_color=colors,$
levels=levels,/cell_fill,/closed,/overplot
map_continents
lats = [-90, -60, -30, 0, 30, 60,90]
lons = [-180,-135,-90,-45,0,45,90,135,180]
; Create string equivalents of latitudes:
latnames = strtrim(lats, 2)
lonnames = strtrim(lons, 2)
; Label the equator:
;latnames(where(lats eq 0)) = '0'
; Draw the grid:
;MAP_GRID,/BOX_AXES, LABEL=2, LATS=lats, LATNAMES=latnames, LATLAB=15,
$
;LONLAB=-2.5, LONDEL=20, LONS=-15, ORIENTATION=-30
MAP_GRID,/BOX_AXES, LABEL=1, LATS=lats, LATNAMES=latnames, $
LONS=lons, LONNAMES=lonnames,LATLAB=7,LONLAB=20
!p.position=[0.2,0.3,0.8,0.7]
colorbar_w,c_color=colors,levels=levels,format='(f10.1)'
END
PRO getcoord,lat,lon
im=72
jm=46
lon(0)=-177.5
for i=1,im-1 do begin
lon(i)= lon(i-1)+5.0
endfor
lat(0)= -89.0
lat(1)= -86.0
for j=2,jm-2 do begin
lat(j)= lat(j-1)+4.0
endfor
lat(45)= 89.0
END
On May 2, 2:07 pm, David Fanning <n...@dfanning.com> wrote:
> t...@atmsci.msrc.sunysb.edu writes:
>> I want to contour data on a map projection and have something like
>> this:
>
>> field=fltarr(72,46)
>> lon=fltarr(72)
>> lat=fltarr(46)
>
>> My lon arrary is indexed from -177.5 to +177.5. My lat array is
>> indexed from -89.0, -86.0, ..., +86.0, +89.0.
>
>> After mapset (using default map and centered at 0,0), I do
>> contour,field,lon,lat. What I get is displaced by 180 degrees. If I
>> change my lon to go from 0.25 to 357.5, then the map looks correct,
>> but the contours do not connect at the dateline.
>
>> Can anyone help with this?
>
> Can we see the commands you are using to display the map
> projection and the contour?
>
> Cheers,
>
> David
> --
> David Fanning, Ph.D.
> Fanning Software Consulting, Inc.
> Coyote's Guide to IDL Programming:http://www.dfanning.com/
> Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
