contour over a map and 180 degrees off or problems with dateline [message #53740] |
Wed, 02 May 2007 10:51  |
teich
Messages: 33 Registered: May 2007
|
Member |
|
|
Hi,
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?
Thanks,
Howie
|
|
|
|
Re: contour over a map and 180 degrees off or problems with dateline [message #53820 is a reply to message #53740] |
Wed, 02 May 2007 12:48  |
teich
Messages: 33 Registered: May 2007
|
Member |
|
|
David, these sort of thing confuse me too!
I wonder if it's possible that contour and map_set could be made more
user friendly one day?
Thanks again,
Howie
On May 2, 3:32 pm, David Fanning <n...@dfanning.com> wrote:
> t...@atmsci.msrc.sunysb.edu writes:
>> If I keep the lon array as is (i.e., -177.5, ..., +177.5) and use your
>> field=Shift(field,36,0) then everything looks in the right location
>> and there is no problem at the dateline. Although this fixed the
>> problem, what if I had an odd number of elements for array lon? Since
>> I compare plots visually, it may be important to get the exact
>> locations. Would you suggest this type of fix (using Shift) for every
>> contour of this type?
>
> Well, I don't know. These kinds of things always confuse me. :-)
>
> It seems to me one of two things has to be wrong. (1) Your
> longitude vector does really describe the physical locations
> of your data (I.e., the left edge of your array is not
> at -177.5 degrees of longitude), or (2) the map projection
> is not centered with respect to the data (I.e., the center of
> the data and the center of the map are at two different places).
>
> I think if you resolve this problem, then shifting, which I would
> only use as a last resort, would probably go away.
>
> 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.")
|
|
|
Re: contour over a map and 180 degrees off or problems with dateline [message #53823 is a reply to message #53740] |
Wed, 02 May 2007 12:32  |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
teich@atmsci.msrc.sunysb.edu writes:
> If I keep the lon array as is (i.e., -177.5, ..., +177.5) and use your
> field=Shift(field,36,0) then everything looks in the right location
> and there is no problem at the dateline. Although this fixed the
> problem, what if I had an odd number of elements for array lon? Since
> I compare plots visually, it may be important to get the exact
> locations. Would you suggest this type of fix (using Shift) for every
> contour of this type?
Well, I don't know. These kinds of things always confuse me. :-)
It seems to me one of two things has to be wrong. (1) Your
longitude vector does really describe the physical locations
of your data (I.e., the left edge of your array is not
at -177.5 degrees of longitude), or (2) the map projection
is not centered with respect to the data (I.e., the center of
the data and the center of the map are at two different places).
I think if you resolve this problem, then shifting, which I would
only use as a last resort, would probably go away.
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.")
|
|
|
Re: contour over a map and 180 degrees off or problems with dateline [message #53827 is a reply to message #53740] |
Wed, 02 May 2007 12:20  |
teich
Messages: 33 Registered: May 2007
|
Member |
|
|
Hi David,
When I say displaced 180 degrees, I mean that features on the map that
should be at 45 East are at 135 West.
If I keep the lon array as is (i.e., -177.5, ..., +177.5) and use your
field=Shift(field,36,0) then everything looks in the right location
and there is no problem at the dateline. Although this fixed the
problem, what if I had an odd number of elements for array lon? Since
I compare plots visually, it may be important to get the exact
locations. Would you suggest this type of fix (using Shift) for every
contour of this type?
I thank you for your help.
Howie
On May 2, 2:52 pm, David Fanning <n...@dfanning.com> wrote:
> t...@atmsci.msrc.sunysb.edu writes:
>> Hi, here are the commands
>
> OK, so what do you mean by "displaced by 180 degrees"?
> Do you mean that the lon vector, which goes from
> -177.500 to 177.5 doesn't really describe where
> the data is located? If you shifted your field
> data in the x direction by 32, would the data look
> correct on your map?
>
> field = Shift(field, 36, 0)
>
> 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.")
|
|
|
Re: contour over a map and 180 degrees off or problems with dateline [message #53834 is a reply to message #53740] |
Wed, 02 May 2007 11:52  |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
teich@atmsci.msrc.sunysb.edu writes:
> Hi, here are the commands
OK, so what do you mean by "displaced by 180 degrees"?
Do you mean that the lon vector, which goes from
-177.500 to 177.5 doesn't really describe where
the data is located? If you shifted your field
data in the x direction by 32, would the data look
correct on your map?
field = Shift(field, 36, 0)
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.")
|
|
|
Re: contour over a map and 180 degrees off or problems with dateline [message #53835 is a reply to message #53740] |
Wed, 02 May 2007 11:41  |
teich
Messages: 33 Registered: May 2007
|
Member |
|
|
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.")
|
|
|
Re: contour over a map and 180 degrees off or problems with dateline [message #53836 is a reply to message #53740] |
Wed, 02 May 2007 11:07  |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
teich@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.")
|
|
|