[gradsusr] percentiles

[Arlindo da Silva]

On Sun, Mar 3, 2013 at 7:39 AM, Piotr Djaków <pdjakow at gmail.com> wrote:

> Hello!
>
> I'm trying to create maps of percentiles from E-OBS dataset. I'm succefuly
> run percentiles.py script from pygrads examples, but i'm not very familiar
> with this.
>
> What i want to do, is make percentiles for given day of a month (i.e. for
> Jan 1st). Input is E-OBS database in one big netCDF file with values from
> jan 1st 1950 to 2012. So i want to extract only values for jan 1st
> 1951,1952,1953,....,2009,2010 and calculate percentiles.
>
> In example we have extracting timeseries:
>
> # Extract a timeseries
> # --------------------
> ga('set t 1 41')
> x = ga.exp('tx')
> g = x.grid
>
>
> But how to extract timeseries for January 1st? It's for example t=1,
> t=366...
>
>
Because leap years can be tricky, I'd exract one day at a time and
concatenate them

from numpy import c_
x = ga.exp('tx(time=1jan1950') # first day in sequence
nx, ny = x.shape
X = x.ravel() # must flatten array for c_ to work
for year in range(1951,2013):
   x = ga.exp('tx(time=1jan%d)'%year)
   X = c_[X,x.ravel()] # concatenate arrays
X = X.reshape((nx,ny,-1)) # restore lon/lat dimensions

This should give you an array with shape (nx,ny,nt) where the "t" dimension
has data only for january 1st of each year. Repeat this for other days.

  Arlindo


-- 
Arlindo da Silva
dasilva at alum.mit.edu
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20130304/a9e6cf09/attachment-0003.html