Re: field for hovmoller diagram

From: Saji N Hameed <saji_at_nyahnyahspammersnyahnyah>
Date: Tue, 6 Oct 2009 18:30:04 +0900

Hi Andrea!

(i don't know how to say "hello" in Italian )

I think gsn_csm_hov does not work as you expected - i.e., create a hovmoeller
plot from your data. It is rather a specialized function that sets some
special resources needed to plot a nice hovmoeller diagram. Which means,
you have to provide the data in the form that it needs (a 2D data set).

Please check the attached function make_and_plot_hovmoller_diagram(wks,var,res,opt). Based on the information provided through the 4th argument 'opt', you
can tell it to cut the data along a specified dimension of your array.
It is currently written to handle on 3D data sets. You need to extent
it if you have a 4D or higher dimesion data set (unlikely). In the
case you have one, you can reduce it to 3D before passing it on. The third
argument should contain the options for plotting your data.

An example of using the function is provided at the bottom.

It will need a number of helper functions, which are found on the top
of the function definition itself.

You can adapt it for your needs...

Best wishes,

saji

---
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
function cr()
begin
  return(inttochar(10))
end
procedure pp(aString)
begin
  print(""+aString)
end
procedure terminate_program(aString)
begin
  print("-----------%< -------")
  pp(aString)
  print("Program Terminated")
  print("-----------%< -------")
  exit
end
function rank_of(var)
local dims
begin
  dims=dimsizes(var)
  return(dimsizes(dims))
end
function reorder_dims(var,newdims)
begin
  rank=rank_of(var)
  if rank.eq.5
    return(var($newdims(0)$|:,$newdims(1)$|:,$newdims(2)$|:,\
                            $newdims(3)$|:,$newdims(4)$|:))
  end if
  if rank.eq.4
    return(var($newdims(0)$|:,$newdims(1)$|:,$newdims(2)$|:,$newdims(3)$|:))
  end if
  if rank.eq.3
    return(var($newdims(0)$|:,$newdims(1)$|:,$newdims(2)$|:))
  end if
  if rank.eq.2
    return(var($newdims(0)$|:,$newdims(1)$|:))
  end if
end
function ensure_last_dim(var,last_dim)
; Will reorder the array to make the
; last dimension to be "dim"
begin
  rank=rank_of(var)
  dims=getvardims(var)
  dindex=ind(dims.eq.last_dim)
  if ismissing(dindex)
    terminate_program("Missing the dimension of "+last_dim)
  end if
  newdims=dims
  newdims(dindex)=dims(rank-1)
  newdims(rank-1)=dims(dindex)
  return(reorder_dims(var,newdims))
end
function cut_from_last_dim(var,slice)
; cuts a slice from the last dimesion
; slice=(/slice1,slice2/)
begin
  slice1=slice(0)
  slice2=slice(1)
  dims=getvardims(var)
  rank=rank_of(var)
  if rank.eq.5
    return(var($dims(0)$|:,$dims(1)$|:,$dims(2)$|:,\
                            $dims(3)$|:,$dims(4)$|slice1:slice2))
  end if
  if rank.eq.4
    return(var($dims(0)$|:,$dims(1)$|:,$dims(2)$|:,$dims(3)$|slice1:slice2))
  end if
  if rank.eq.3
    return(var($dims(0)$|:,$dims(1)$|:,$dims(2)$|slice1:slice2))
  end if
  if rank.eq.2
    return(var($dims(0)$|:,$dims(1)$|slice1:slice2))
  end if
end
function make_and_plot_hovmoller_diagram(wks,var,res,opt)
begin
  rank=rank_of(var)
  if rank .gt. 3
    terminate_program("Extend me to handle >= 4 dimension arrays")
  end if
  if rank.eq.2
    return(gsn_csm_contour(wks,var,res))
  end if  
  if .not. isatt(opt,"avg_over")
    msg1="You need to specify which dimension to average"
    msg2="Example: opt_at_avg_over="+"x"
    terminate_program(msg1+cr+msg2)
  end if
  dim=opt_at_avg_over
  if .not. isatt(opt,dim)
    msg1="You need to specify the slice of "+dim+" to average"
    msg2="Example: opt_at_x=(/0,10/)"
    terminate_program(msg1+cr+msg2)
  end if
  slice=opt@$dim$
  new_var=cut_from_last_dim(ensure_last_dim(var,dim),slice)
printVarSummary(dim_avg_Wrap(new_var))
  return(gsn_csm_contour(wks,dim_avg_Wrap(new_var),res))
end
; Bogus 3D Data -- stolen from http://www.ncl.ucar.edu/Applications/Scripts/conwomap_3.ncl
; -------------------------------------------------
  mx  = 31
  ny  = 21 
  nt  = 9 
	
  u   = random_normal(0.0, 3.0, (/ny,mx,nt/))  
  u!0 = "y"
  u!1 = "x"
  u!2 = "t"
  u&x = ispan(0,mx-1,1)
  u&y = ispan(0,ny-1,1)
  u&t = ispan(0,nt-1,1)
  u@_FillValue = -999.
wks=gsn_open_wks("x11","hov")
res=True
res_at_tiMainString = "Default Hovmueller"       ; title
opt=True
opt_at_avg_over="y"
opt_at_y=(/0,10/)
plot = make_and_plot_hovmoller_diagram(wks,u,res,opt)
* Andrea Borrelli <andrea.blizzard_at_gmail.com> [2009-10-04 11:35:34 +0200]:
> Hello to all,
> 
> I'd plot an Hovmoller diagram on Tropical Pacific of the field wind stress U
> (is an OPA8.2 output), but I don't understand how "gsn_csm_hov" works.
> 
> Wind stress U is in this form:
> 
> Variable: stress_u2
> Type: float
> Total Size: 535760 bytes
>             133940 values
> Number of Dimensions: 3
> Dimensions and sizes:   [time | 5] x [y | 148] x [x | 181]
> Coordinates:
>             time: [   0..   4]
> Number Of Attributes: 10
>   long_name :   Wind Stress along i
>   units :       N/m2
>   coordinates : nav_lon nav_lat
>   valid_min :   1.000000020040877e+20
>   valid_max :   -1.000000020040877e+20
>   short_name :  sozotaux
>   online_operation :    ave(only(x))
>   interval_operation :  5400
>   interval_write :      2592000
>   _FillValue :  1e+20
> 
> This field have three dimensions but "gsn_csm_hov" function it requires only
> two.
> What could I do?
> 
> Thanks a lot for your help
> Andrea
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-- 
Saji N. Hameed
APEC Climate Center          				
1463 U-dong, Haeundae-gu,                               +82 51 745 3951
BUSAN 612-020, KOREA                    		saji_at_apcc21.net
Fax: +82-51-745-3999
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Oct 06 2009 - 03:30:04 MDT

This archive was generated by hypermail 2.2.0 : Fri Oct 09 2009 - 08:29:22 MDT