| Re: Contour data with two different lon/lat grids onto the same map projection [message #71223] |
Fri, 04 June 2010 10:07  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
mgalloy writes:
> I use main-level programs as a place to start or to place an example at
> the end of a file containing regular routines. Anything more complex
> than a few lines should be in a function/procedure. But the nice part is
> that it is easy to just add the "pro/function my_routine, a, b, etc.."
> line at the top of the main-level program to change it into a routine.
> After coding and getting something working, 99+% of what gets committed
> and kept is in regular routines.
I use main-level programs this way, too. But I see plenty
of people who use them (misuse them?) almost exclusively
and I think that is a big mistake, since it almost always
leads to unanticipated errors once you get more than about
10 lines of code in the darn things.
I certainly didn't mean to pick on you, but I'm almost
certain that if you can't run a program twice and get
the same result, it is programmer error in a main program. ;-)
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 data with two different lon/lat grids onto the same map projection [message #71239 is a reply to message #71223] |
Thu, 03 June 2010 19:22  |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
Tess writes:
> However, if I load the variables from netcdf to idl variables, then
> run this script; it plots fine. If I immediately rerun the same script
> in the same idl session, I get the error:
> "CONTOUR: X,Y, or Z array dimensions are incompatible."
> The first couple of times I ran this script I was adding to it in
> stages, and running it to see what the results were; the error was a
> bit misleading in that way, as I thought it was the result of adding
> the pv grids to the rainfall grids I already had working on the plot.
> Hence my appeal to David and others for "ideas"! If I exit the IDL
> session and reload the variables in a new session to rerun the script,
> it plots fine first time. But only once per session!
>
> Not sure why this happens, but unless it's immediately obvious to
> anyone it's probably not worth pursuing...
Well, I take it by "script" you mean an IDL main
program (that is, not a procedure or function).
There is plenty in just this small code sample
to cause me to well believe such a script might
cause you trouble.
In fact, in spite of Mike Galloy proclaiming here the
other day that he often uses main programs, I rarely
do. And for exactly this reason. They can get you
in a world of strange hurt. You would probably find
out in a minute or two why you can't run the program
more than once if you made a procedure out of it.
I'm mostly confused by statements like this, though:
"if I load the variables from netcdf to idl variables..."
When you read the netCDF variables out of the file they ARE
IDL variables. No need to load anything. Any why is the
graphic being drawn over and over again in the same window
in a loop!?
I'm just not sure what to make of all this. Especially when
I see the limits on the MAP_SET command, which have no logical
connection to the latitude and longitude coordinates discussed
in your previous communications.
Perhaps this is just too small a code fragment to visualize the
rest of it. But, I would say, if it works at least once, I'd be
tempted to let it go as a good show. :-)
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 data with two different lon/lat grids onto the same map projection [message #71240 is a reply to message #71239] |
Thu, 03 June 2010 15:51 |
tessp
Messages: 8
Registered: June 2010
|
Junior Member |
|
|
On Jun 3, 9:46 pm, David Fanning <n...@dfanning.com> wrote:
> pp writes:
>> You need to say more about what you are doing and why it does not
>> work. It is (nearly) irrelevant that the two variables are in
>> different grids. You can just overplot them.
>
> Yes, I agree. The OVERPLOT keyword is basically
> all that is needed here, which is why I was confused
> about what kind of "ideas" you wanted.
>
> You set your map data space, then simply overplot your
> two contours onto it. What grid they are on is irrelevant.
>
> More details as to the problems you have, please. :-)
>
> Cheers,
>
> David
>
This is indeed what I do: as extract from script shows -
****************************
; First use cdf2idl process to read in precip, then pv files.
; Rainfall has total per day, pv has four timesteps per day.
;
; Get the netCDF data into IDL arrays
; Precip:
lon_g=fltarr(360)+lon
lat_g=fltarr(180)+lat
date_g=fltarr(22)+date
pcp=fltarr(360,180,22)+pcp
lonndx_g=long(size(lon, /n_elements))
latndx_g=long(size(lat, /n_elements))
datendx_g=long(size(date, /n_elements))
; PV:
lat=G0_LAT_2
lon=G0_LON_3
time=INITIAL_TIME0_ENCODED
pv=PV_GDS0_THEL*10^6. ; pv now in PVU of 10^-6 K m^2 kg^-1
s^-1
lvl=LV_THEL1
lonndx=long(size(lon, /n_elements))
latndx=long(size(lat, /n_elements))
timendx=long(size(time, /n_elements))
k=6 ; 350K level only
date1=string(time[k])
date=strmid(date1,0,8)
utc=strmid(date1,7,2)
subb=string(date)+' '+string(utc)
counter=0L
daycounter=-0.25
For l=24,timendx-13 Do Begin ; to match rainfall data days
available in gpcp dataset
daycounter=daycounter+0.25 ; 4 PV to each rainfall day
map_set,-25,130,0,LIMIT=[-5,100,-45,160], /isotropic, /mercator,
position=[0.1,0.1,0.8,0.8]
contour,pcp[*,*,
(long(daycounter))],lon_g,lat_g,levels=[((findgen(21)+1)*5.) ,200], $
c_colors=(indgen(21)+1)+4,/isotropic, /cell_fill, $
c_charsize=0.2, /overplot, /noerase
contour, pv[*,*,k,l],lon,lat, levels=((findgen(25)*0.5)), $ ;
positive PV
c_colors=0,/isotropic, /follow,
c_labels=[Replicate(1,25)], $
c_charsize=0.5, /overplot, /noerase, color=0
contour, pv[*,*,k,l],lon,lat, levels=(((findgen(25)-25)*0.5)), $ ;
negative PV
c_colors=0,/isotropic, /follow,
c_labels=[Replicate(1,25)], $
c_charsize=0.5, /overplot, /noerase, color=0
map_continents, color=0, /hires
*********************************
However, if I load the variables from netcdf to idl variables, then
run this script; it plots fine. If I immediately rerun the same script
in the same idl session, I get the error:
"CONTOUR: X,Y, or Z array dimensions are incompatible."
The first couple of times I ran this script I was adding to it in
stages, and running it to see what the results were; the error was a
bit misleading in that way, as I thought it was the result of adding
the pv grids to the rainfall grids I already had working on the plot.
Hence my appeal to David and others for "ideas"! If I exit the IDL
session and reload the variables in a new session to rerun the script,
it plots fine first time. But only once per session!
Not sure why this happens, but unless it's immediately obvious to
anyone it's probably not worth pursuing...
Cheers,
Tess
|
|
|
|
|
|
|
|
|
|
| Re: Contour data with two different lon/lat grids onto the same map projection [message #71255 is a reply to message #71254] |
Wed, 02 June 2010 21:19 |
tessp
Messages: 8
Registered: June 2010
|
Junior Member |
|
|
> Ideas for what? Reducing visual clutter? I'm not at all
> sure what you want ideas about. :-)
>
> Cheers,
>
> David
>
> P.S. Of course, it is possible. The question you should
> be asking is this: Is it a good idea?
>
> --
> 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.")
It's a very good idea. Looking for a link between potential vorticity
and rainfall in certain areas of the globe. Unfortunately, the best
datasets for the two variables come from different projects by
different agencies. If I line contour the vorticity and shade the
rainfall, it should highlight the dynamics of the atmospheric state
without too much visual clutter.... I have done it for a couple of
other situations, fortunately on the same grids, and it's quite
informative! :o)
What I need is some ideas on how to do the contour plots, as my
current method has incompatible x,y,z coordinates due to the grid
mismatches (I think).
Any help from Coyote gratefully received as always. Followed by
payment to Doctors Without Borders, if that's still charity of choice?
Regards
Tess
|
|
|
|
| Re: Contour data with two different lon/lat grids onto the same map projection [message #71256 is a reply to message #71255] |
Wed, 02 June 2010 21:06 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
tessp writes:
> Not sure if this is even possible. I have two sets of meteorological
> data, one for rainfall and one for vorticity. From different sources,
> they are available on different latitude and longitude grids (i.e.
> rainfall on -89.5 to 89.5 lat, 0.5 to 259.5 lon: vorticity on 90 to
> -90 lat, 0 - 358.5 lon), with different grid spacing (1 deg vs 1.5
> deg) just for added fun. I want to plot both sets of data as contours
> over one (global) map projection. Any ideas anyone?
Ideas for what? Reducing visual clutter? I'm not at all
sure what you want ideas about. :-)
Cheers,
David
P.S. Of course, it is possible. The question you should
be asking is this: Is it a good idea?
--
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 data with two different lon/lat grids onto the same map projection [message #71300 is a reply to message #71239] |
Tue, 08 June 2010 18:59 |
tessp
Messages: 8
Registered: June 2010
|
Junior Member |
|
|
> Well, I take it by "script" you mean an IDL main
> program (that is, not a procedure or function).
> There is plenty in just this small code sample
> to cause me to well believe such a script might
> cause you trouble.
>
> In fact, in spite of Mike Galloy proclaiming here the
> other day that he often uses main programs, I rarely
> do. And for exactly this reason. They can get you
> in a world of strange hurt. You would probably find
> out in a minute or two why you can't run the program
> more than once if you made a procedure out of it.
>
> I'm mostly confused by statements like this, though:
>
> "if I load the variables from netcdf to idl variables..."
>
> When you read the netCDF variables out of the file they ARE
> IDL variables. No need to load anything. Any why is the
> graphic being drawn over and over again in the same window
> in a loop!?
>
> I'm just not sure what to make of all this. Especially when
> I see the limits on the MAP_SET command, which have no logical
> connection to the latitude and longitude coordinates discussed
> in your previous communications.
>
> Perhaps this is just too small a code fragment to visualize the
> rest of it. But, I would say, if it works at least once, I'd be
> tempted to let it go as a good show. :-)
>
> Cheers,
>
> David
>
> --
> David Fanning, Ph.D.
David - Yes, script is my word for main program. I write a heap of
different programs to do things, and haven't ever got around to making
them into proceudres; although I agree with Mike that I wouldn't
change the code much if at all, were I to do so....This is mainly due
to time pressures, and the fact that once I have done something and
got the plots I need it's generally on to the next thing on the list.
The code is generally not pretty, but it works - usually. :o)
When I talk about loading the netcdf variables to IDL variables that
is just my clumsy terminology for reading the data out of the netcdf
file.
The map is being drawn over and over again to a postscript file in a
loop because I am plotting each timestep of the data, usually every 6
hours, to get a sequence of maps showing the development of the
atmospheric variables I'm looking at. Stepping through the .ps file
then shows how weather systems develop over the time period of
interest.
The limits on the map_set command are to map just the Australian
region. In my full code I have options to do this, or the whole globe,
or just the southern hemisphere.
Andy from Reading has my full code and data, and is taking a look at
this for me....so he may have some insight into the issues I'm having
when running it more than once.
|
|
|
|
| Re: Contour data with two different lon/lat grids onto the same map projection [message #71365 is a reply to message #71300] |
Thu, 10 June 2010 16:01 |
tessp
Messages: 8
Registered: June 2010
|
Junior Member |
|
|
>
> Andy from Reading has my full code and data, and is taking a look at
> this for me....so he may have some insight into the issues I'm having
> when running it more than once.
Thanks, Andy, for pointing out that I have managed to confuse the two
latitude and longitude arrays. So yes David, you were right! -
"I'm almost certain that if you can't run a program twice and get
the same result, it is programmer error in a main program. ;-)"
Good advice from Andy to do the basic checks by looking at the lat and
lon arrays at various stages of my code, which soon points up the
error. Many thanks for your help Andy, much appreciated.
Thanks David and Mike for your input also.
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
