scatter plot of directions for ensemble data at different time steps

Ufficio Meteo cmt.meteo at ARPA.VENETO.IT
Wed Mar 4 03:25:28 EST 2009


Hello all.

I'm trying to plot ensemble ECMWF data for wind.
It is not a much difficult job to plot the time series of standard deviation and mean for intensities.
The problem comes with directions. Sorry, it could be just a mistake of mine..... but I cannot identify it!
I'm working with version 1.8.

I start with 51 grib files (one for each ensemble member), where latitude, longitude, level and time vary.
I want to have a diagram for each grid point and level, so I let only time to vary (in 40 steps) when I run the plotting job.
An example of the resulting diagram is attached as a png file.
I can have a scatter plot for the (u,v) couples, but it seems that I must not set time, even if such a setting covers the whole time range. Otherwise, the plot is empty. That is, I must avoid any 'set t .....' command in the part of the script ruling the scatter plot, anyway I obtain too many points in the plot (maybe all members x all times = 2040 points in the u-v plane!!!!!).
How can I have one plot for each time step, instead of a plot for the whole time range?

This is one of the ctl files ("p2.ctl", the grib file is attached):

*************************************************************************************************
dset ^p2.grib

index ^p2.grib.idx

undef 9.999E+20

title p2.grib

* produced by grib2ctl v0.9.12.5p41

dtype grib 255

options yrev

ydef 4 linear 45.000000 0.5

xdef 6 linear 10.500000 0.500000

tdef 40 linear 18Z23feb2009 6hr

* z has 3 levels, for prs

zdef 3 levels

1000 925 850

vars 12

no10Usfc 0 165,1,0 ** surface 10 metre u wind component m s**-1

no10Vsfc 0 166,1,0 ** surface 10 metre v wind component m s**-1

no2Tsfc 0 167,1,0 ** surface 2 metre temperature K

SFsfc 0 144,1,0 ** surface Snow fall m(of water equivalent)

Tprs 3 130,100,0 ** (profile) Temperature K

TPsfc 0 228,1,0 ** surface Total precipitation m

U850mb 0 131,100,850 ** 850 mb U-velocity m s**-1

V850mb 0 132,100,850 ** 850 mb V-velocity m s**-1

Z850mb 0 129,100,850 ** 850 mb Geopotential (at the surface = orography) m**2 s**-2

var121sfc 0 121,1,0 ** surface undefined

var122sfc 0 122,1,0 ** surface undefined

var123sfc 0 123,1,0 ** surface undefined

ENDVARS

*************************************************************************+

Then, here is the grads job to plot the data:

********************************************************************************************************************* 
function rose(args)

'reinit'

llon=subwrd(args,1)

llat=subwrd(args,2)

llev=subwrd(args,3)

*Scelgo il titolo con la localita' corrispondente alle coordinate

'open c.ctl'

'set lat 'llat

'set lon 'llon

nomepunto=''

coordinate=llat' 'llon

if (coordinate='46.00 12.00') ; nomepunto='FELTRE'; endif

if (coordinate='46.50 12.00') ; nomepunto='LIVINALLONGO'; endif

if (coordinate='46.50 12.25') ; nomepunto='MISURINA'; endif

if (coordinate='46.50 12.50') ; nomepunto='AURONZO'; endif

if (coordinate='46.25 12.00') ; nomepunto='AGORDO'; endif

if (coordinate='46.25 12.25') ; nomepunto='LONGARONE'; endif

if (coordinate='46.25 12.50') ; nomepunto='CIMOLAIS'; endif

if (coordinate='46.00 11.50') ; nomepunto='MARCESINA'; endif

if (coordinate='46.00 11.75') ; nomepunto='ARSIE'; endif

if (coordinate='46.00 12.00') ; nomepunto='FELTRE'; endif

if (coordinate='46.00 12.25') ; nomepunto='VITTORIO VENETO'; endif

if (coordinate='45.75 11.25') ; nomepunto='VALLI DEL PASUBIO'; endif

if (coordinate='45.75 11.50') ; nomepunto='CARRE VI'; endif

if (coordinate='45.75 11.75') ; nomepunto='BASSANO'; endif

if (coordinate='45.75 12.00') ; nomepunto='MONTEBELLUNA'; endif

if (coordinate='45.75 12.25') ; nomepunto='VILLORBA'; endif

if (coordinate='45.75 12.50') ; nomepunto='ODERZO'; endif

if (coordinate='45.75 12.75') ; nomepunto='PORTOGRUARO'; endif

if (coordinate='45.75 13.00') ; nomepunto='SAN MICHELE AL TAGLIAMENTO'; endif

if (coordinate='45.50 10.75') ; nomepunto='LAZISE'; endif

if (coordinate='45.50 11.00') ; nomepunto='VR TORRICELLE'; endif

if (coordinate='45.50 11.25') ; nomepunto='MONTECCHIA DI CROSARA'; endif

if (coordinate='45.50 11.50') ; nomepunto='ALTAVILLA VIC.NA'; endif

if (coordinate='45.50 11.75') ; nomepunto='GRISIGNANO DI ZOCCO'; endif

if (coordinate='45.50 12.00') ; nomepunto='PIANIGA'; endif

if (coordinate='45.50 12.25') ; nomepunto='MESTRE'; endif

if (coordinate='45.50 12.50') ; nomepunto='CAVALLINO'; endif

if (coordinate='45.50 12.75') ; nomepunto='LARGO DI JESOLO'; endif

if (coordinate='45.25 11.00') ; nomepunto='ISOLA DELLA SCALA'; endif

if (coordinate='45.25 11.25') ; nomepunto='ROVERCHIARA'; endif

if (coordinate='45.25 11.50') ; nomepunto='MONTAGNANA'; endif

if (coordinate='45.25 11.75') ; nomepunto='MONSELICE'; endif

if (coordinate='45.25 12.00') ; nomepunto='BOVOLENTA'; endif

if (coordinate='45.25 12.25') ; nomepunto='LAGUNA SUD'; endif

if (coordinate='45.00 11.25') ; nomepunto='CARBONARA DI PO'; endif

if (coordinate='45.00 11.50') ; nomepunto='BAGNOLO DI PO'; endif

if (coordinate='45.00 11.75') ; nomepunto='ARQUA POLESINE'; endif

if (coordinate='45.00 12.00') ; nomepunto='ADRIA'; endif

if (coordinate='45.00 12.25') ; nomepunto='PORTO VIRO'; endif

if (coordinate='45.00 12.50') ; nomepunto='PILA'; endif

if (coordinate='45.75 10.75') ; nomepunto='GARDA NORD'; endif

if (coordinate='45.75 11.00') ; nomepunto='ALA (TN)'; endif

*Epsigramma in basso, con la statistica del modulo del vento

*'set parea 1 8 0.5 5'

'set parea 1 11 0.5 3.75'

*Definisco le variabili da plottare nel diagramma in basso:

*media e varianza del modulo del vento

'set lev 'llev

'set t 1 40'

if llev=1000

'define campoc=sqrt(no10Usfc*no10Usfc+no10Vsfc*no10Vsfc)'

'define campocu=no10Usfc'

'define campocv=no10Vsfc'

else

'define campoc=sqrt(U850mb*U850mb+V850mb*V850mb)'

'define campocu=U850mb'

'define campocv=v850mb'

endif

i=1

'define ps=0.0'

'define pv=0.0'

while(i<51)

'open p'i'.ctl'

'set dfile 2'

'set lat 'llat

'set lon 'llon

'set lev 'llev

'set t 1 40'

if llev=1000

'define campop'i'=sqrt(no10Usfc*no10Usfc+no10Vsfc*no10Vsfc)'

'define campop'i'u=no10Usfc'

'define campop'i'v=no10Vsfc'

else

'define campop'i'=sqrt(U850mb*U850mb+V850mb*V850mb)'

'define campop'i'u=U850mb'

'define campop'i'v=V850mb'

endif

'define ps=ps+campop'i

'close 2'

i=i+1

endwhile

'define ps=ps+'campoc

'define pm=ps/51.0'

*Plot delle medie nel diagramma in basso

'set grads off'

'set gxout line'

i=1

while(i<51)

'set vrange 0 20'

'set ccolor 1'

'set cmark 3'

'set cthick 1'

'set csmooth on'

'set missconn on'

'd campop'i

i=i+1

endwhile

i=1

while(i<51)

'define pvd=pm-campop'i

'define pv=pv+pvd*pvd'

i=i+1

endwhile

'define pvd=pm-'campoc

'define pv=pv+pvd*pvd'

'define pstd=sqrt(pv/50.0)'

'set ccolor 2'

'set cmark 3'

'set cthick 3'

'set csmooth on'

'set missconn on'

'set vrange 0 20'

'd campoc'

'set ccolor 4'

'set cmark 1'

'set cthick 3'

'set cstyle 0'

'set missconn off'

'd pm'

*Barre di deviazione standard nel diagramma in basso

'set gxout bar'

'set bargap 50'

'set baropts outline'

'set barbase 0'

'set ccolor 4'

'set cthick 10'

'define pp1=pm-pstd/2.0'

'define pp2=pm+pstd/2.0'

'set vrange 0 20'

'd pp1;pp2'

*Diagrammi in alto: scatter plots con le punte dei vettori vento

*'set parea 1 8 5.5 10.5'

'set parea 1 11 4.25 7.5'

*Plot nel diagramma in alto

'set grads off'

'set grid off'

'set xlab off'

'set ylab off'

'set gxout scatter'

'set vrange -20 20'

'set vrange2 -20 20'

i=1

while(i<51)

'set ccolor 1'

'set digsize 0.05'

'set cmark 3'

'd campop'i'u;campop'i'v'

i=i+1

endwhile

'set ccolor 2'

'set digsize 0.05'

'set cmark 3'

'd campocu;campocv'

'draw title 'nomepunto

'enable print rose.gmf'

'print'

'disable print'

'quit'

*******************************************************************************************

Antonino Claudio Bonan

Meteorologo

Centro Meteorologico di Teolo - ARPAV

via Marconi 55

35037 Teolo (PD) - Italy


I'm protected by SpamBrave
http://www.spambrave.com/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20090304/ba267893/attachment.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: rose_sfc_11.50_45.00.png
Type: image/png
Size: 92785 bytes
Desc: not available
Url : http://gradsusr.org/pipermail/gradsusr/attachments/20090304/ba267893/attachment.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: p2.grib
Type: application/octet-stream
Size: 105600 bytes
Desc: not available
Url : http://gradsusr.org/pipermail/gradsusr/attachments/20090304/ba267893/attachment.obj 


More information about the gradsusr mailing list