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