Re: [ferret_users] event statistics -- length of contiguous "spell"

Hi Russ,

Thank you for suggetions.

I agree that your 1st suggestion makes the logic generally valid

(though with the particular data I am working with it starts from 1st Jan (which was missing but I filled all missing with value=0) and by defintion of active and break spell they can occur only during Jul and Aug so by design of the dataset 1st time step is always 0)

I tried to implement rest of your suggestions but got following error:

yes? def axis/x=1:`id_good[l=@max]`:1 xax
yes? let id_good_2d = if id_good eq i then id_good + 0*i[gx=xax]
yes? shade id_good_2d
*** NOTE: Ambiguous coordinates on X axis: IF ID_GOOD EQ I THEN ID_GOOD + 0*I[GX=XAX]
*** NOTE: Ambiguous coordinates on X axis: IF ID_GOOD EQ I THEN ID_GOOD + 0*I[GX=XAX]
**ERROR: illegal limits: I on grid (G008) is not in the range X=66.3:98.8
Axis extremes are X=0.5:58.5

What could be the problem?

Thanks,

Jagadish

On Thursday, May 15, 2014 2:38 AM, Russ Fiedler <russell.fiedler@xxxxxxxx> wrote:

Hi Jagadish,

I think you also need to consider the case when the first entry of your time series is 1 (i.e. coincides with the start of a spell). In that case your solution (2) would not catch the first event and even numbers would correspond to wet spells rather than dry.

You could either pad the series at the start with a dry day or maybe add one if necessary using the first value in your range of interest.

Something like

def sym first_day 183 ! this is a wet day in your series below
let ID_GOOD = VAR1[L=@EVNT:0.999] + `VAR1[l=(first_day)]`

list/l=($first_day):285 ...

To get the number of days in a spell you could spread the results along a new axis and then count the number of good values.

def axis/x=1:`id_good[l=@max]`:1 xax
let id_good_2d = if id_good eq i then id_good + 0*i[gx=xax] ! Put spell "i" on its own axis.

shade id_good_2d ! we can see each spell wet or dry. This can be masked using mod(var,2) eq 0 or 1

let length_of_spell = id_good_2d[l=@ngd]

plot/title="Spell length" length_of_spell

plot/title="Wet spell length"/sym if mod(id_good,2) eq 1 then length_of_spell

The stats should be easy after this.

Cheers,
Russ

On 15/05/14 05:22, 'jagadish karmacharya' via _OAR PMEL Ferret Users wrote:

Dear steve and Hein,

Thanks for quickly sharing your thoughts.

I have not really gone though all the details you suggestioned but I am glad to report that I have found a solution (in fact 2 different solutions).

First a clearification - What I mean by the spells was: active and break spell has been defined as days when rainfall is above or below a certain threshold provided it last for atleast 3 consequtive days. Previous taking help from the ferret forum I was able implement the criteria such that the result timeseries has value = 1 for active spell (also for break spell but in seperate variable) and missing for other period. That has worked fine as I was then only interested in the total duration of such spell and other composite analysis using those. So my original data and variables looks like following:

yes? sh d
     currently SET data sets:
    1> ../results/12kmRuns/IMD_active_break_JJASonly.90_05.nc (default)
name     title                             I         J         K         L         M         N
ATV_JA   ATV                              ...       ...       ...       1:5844    ...       ...
       (X=66.3E:98.8E, Y=6.3N:38.8N)
BRK_JA   BRK                              ...       ...       ...       1:5844    ...       ...
       (X=66.3E:98.8E, Y=6.3N:38.8N)

yes? list atv_ja[l=@ngd],brk_ja[l=@ngd]
(deleted lines)
         ATV_JA BRK_JA
I / *:     108.0 125.0

Solution 1)
Hein email's yesterday indicating it would be difficult to implement (what I have request for) in ferret prompted me to seek other options which lead to 1st solution in which I avoided the proble mby creating new files for active (and break) period wherein values > 1 during active period and zero for the rest. Next I used:
let ID_atv=atv[l=@evnt:1]

to "tag" the events where successive odd id's corresponds to succesive active spells (& succesive even id's corresponds to succesive gaps between the spells).

Now from that I can easly find duration of each spell. (I then imported those spell lengh in spreadsheet to compute mean, median, max of the spell duration etc, but it should be possible to do it in ferret itself with little more play)

Solution 2)
But quickly going through steve's mail a while ago I though it might be posible to use similar trick:
so what I did is - I filled the missing values in the original data (above) with zero. Then defined the event as before with slight twist:
let id1=var1[l=@evnt:0.999]        ! here any value between 0 and 1 would work
which gave desired result!

i.e.
let VAR1 = IF ATV_JA THEN ATV_JA ELSE 0
let ID_GOOD = VAR1[L=@EVNT:0.999]
let ID_BAD = ATV_JA[L=@EVNT:1]        ! just for illustration, this doesn't work, using value between 0 & 1 also doesn't give desired result

list/l=181:245 atv_ja,id_bad,var1,id_good
(deleted lines)
                     ATV_JA ID_BAD VAR1 ID_GOOD
30-JUN-1990 00 / 181:   ....   0.00 0.000   0.000
01-JUL-1990 00 / 182:   ....   0.00 0.000   0.000
02-JUL-1990 00 / 183: 1.000   1.00 1.000   1.000
03-JUL-1990 00 / 184: 1.000   2.00 1.000   1.000
04-JUL-1990 00 / 185: 1.000   3.00 1.000   1.000
05-JUL-1990 00 / 186: 1.000   4.00 1.000   1.000
06-JUL-1990 00 / 187:   ....   4.00 0.000   2.000
07-JUL-1990 00 / 188:   ....   4.00 0.000   2.000
(deleted lines)
19-AUG-1990 00 / 231:   ....   4.00 0.000   2.000
20-AUG-1990 00 / 232:   ....   4.00 0.000   2.000
21-AUG-1990 00 / 233:   ....   4.00 0.000   2.000
22-AUG-1990 00 / 234: 1.000   5.00 1.000   3.000
23-AUG-1990 00 / 235: 1.000   6.00 1.000   3.000
24-AUG-1990 00 / 236: 1.000   7.00 1.000   3.000
25-AUG-1990 00 / 237: 1.000   8.00 1.000   3.000
26-AUG-1990 00 / 238:   ....   8.00 0.000   4.000
27-AUG-1990 00 / 239:   ....   8.00 0.000   4.000
28-AUG-1990 00 / 240:   ....   8.00 0.000   4.000
29-AUG-1990 00 / 241: 1.000   9.00 1.000   5.000
30-AUG-1990 00 / 242: 1.000 10.00 1.000   5.000
31-AUG-1990 00 / 243: 1.000 11.00 1.000   5.000
01-SEP-1990 00 / 244:   .... 11.00 0.000   6.000
02-SEP-1990 00 / 245:   .... 11.00 0.000   6.000

I would like to thank your both for sharing your ideas that help me come up with these solutions.

Best regards,

Jagadish

On Wednesday, May 14, 2014 6:39 PM, Steve Hankin <steven.c.hankin@xxxxxxxx> wrote:

Hi Jagadish,

Here is an example of using the "closest index above" and "closest index below" filters in Ferret to compute the length of contiguous "spells"

define axis/t=1-jan-1900:1-jan-2000:30/units=days t30let spells = IF RANDU(L[gt=t30]) GT 0.5 THEN 1let inverseSpells = IF MISSING(spells,2) EQ 2 THEN 1let closestAbove = inverseSpells[l=@cia]let closestBelow = inverseSpells[l=@cib]let rawSize = closestAbove+closestBelowlet contigSize = IF rawSize NE 0 THEN rawSize-1 list/l=1:50 spells, inverseSpells, closestAbove, closestBelow, contigSize TIME: 17-DEC-1899 00:00 to 26-JAN-1904 00:00 Column 1: SPELLS is IF RANDU(L[GT=T30]) GT 0.5 THEN 1 Column 2: INVERSESPELLS is IF MISSING(SPELLS,2) EQ 2 THEN 1 Column 3: CLOSESTABOVE is INVERSESPELLS[L=@CIA] Column 4: CLOSESTBELOW is INVERSESPELLS[L=@CIB] Column 5: CONTIGSIZE is IF RAWSIZE NE 0 THEN RAWSIZE-1 SPELLS INVERSE CLOSEST CLOSEST CONTIGSIZE01-JAN-1900 / 1: .... 1.000 0.000 0.000 ....31-JAN-1900 / 2: 1.000 .... 1.000 1.000 1.00002-MAR-1900 / 3: .... 1.000 0.000 0.000 ....01-APR-1900 / 4: 1.000 .... 1.000 1.000 1.00001-MAY-1900 / 5: .... 1.000 0.000 0.000 ....31-MAY-1900 / 6: 1.000 .... 1.000 1.000 1.00030-JUN-1900 / 7: .... 1.000 0.000 0.000 ....30-JUL-1900 / 8: 1.000 .... 1.000 1.000 1.00029-AUG-1900 / 9: .... 1.000 0.000 0.000 ....28-SEP-1900 / 10: .... 1.000 0.000 0.000 ....28-OCT-1900 / 11: .... 1.000 0.000 0.000 ....27-NOV-1900 / 12: .... 1.000 0.000 0.000 ....27-DEC-1900 / 13: 1.000 .... 1.000 1.000 1.00026-JAN-1901 / 14: .... 1.000 0.000 0.000 ....25-FEB-1901 / 15: 1.000 .... 2.000 1.000 2.00027-MAR-1901 / 16: 1.000 .... 1.000 2.000 2.00026-APR-1901 / 17: .... 1.000 0.000 0.000 ....26-MAY-1901 / 18: .... 1.000 0.000 0.000 ....25-JUN-1901 / 19: 1.000 .... 2.000 1.000 2.00025-JUL-1901 / 20: 1.000 .... 1.000 2.000 2.00024-AUG-1901 / 21: .... 1.000 0.000 0.000 ....23-SEP-1901 / 22: .... 1.000 0.000 0.000 ....23-OCT-1901 / 23: .... 1.000 0.000 0.000 ....22-NOV-1901 / 24: .... 1.000 0.000 0.000 ....22-DEC-1901 / 25: .... 1.000 0.000 0.000 ....21-JAN-1902 / 26: 1.000 .... 2.000 1.000 2.00020-FEB-1902 / 27: 1.000 .... 1.000 2.000 2.00022-MAR-1902 / 28: .... 1.000 0.000 0.000 ....21-APR-1902 / 29: 1.000 .... 1.000 1.000 1.00021-MAY-1902 / 30: .... 1.000 0.000 0.000 ....20-JUN-1902 / 31: 1.000 .... 3.000 1.000 3.00020-JUL-1902 / 32: 1.000 .... 2.000 2.000 3.00019-AUG-1902 / 33: 1.000 .... 1.000 3.000 3.00018-SEP-1902 / 34: .... 1.000 0.000 0.000 ....18-OCT-1902 / 35: .... 1.000 0.000 0.000 ....17-NOV-1902 / 36: 1.000 .... 2.000 1.000 2.00017-DEC-1902 / 37: 1.000 .... 1.000 2.000 2.00016-JAN-1903 / 38: .... 1.000 0.000 0.000 ....15-FEB-1903 / 39: .... 1.000 0.000 0.000 ....

On 5/14/2014 8:56 AM, Steve Hankin wrote:

On 5/13/2014 5:19 PM, 'jagadish karmacharya' via _OAR PMEL Ferret Users wrote:

Dear ferret users,

I have a 1 dimensional time series with intermittent values equal to 1 and missing values in between (1's represent spells of certain events). I want to calculate number of spells of of 1's and other stats such as mean and medium duration of spells of 1's.

Hi Jagadish,

Following up on the ideas that Hein has begun ...

Your description is not sufficiently specific on the definition of a "spell" to give a precise answer on how (or whether) the tools available in Ferret can do what you are after.

if the definition of a spell is that it a threshold on the *density* of 1's found per unit length of the time series, then the @SBX (running average) filter will give you something very close to what you are after

e.g. yes? let tdensity = mytimeseries[l=@sbx:11]
       yes? PLOT tdensity    ! have a look ...
       yes? let myspells = IF tdensity gt threshold THEN 1

if the definition of a spell is that it must be a strictly contiguous sequence of 1's, it will take more imagination to see how to apply the available Ferret tools. Here's one approach to identifying the contiguous spells using existing tools:

use the tdensity approach above varying the box width to identify the longest spell.
yes? let tdensity = mytimeseries[l=@sbx:11]
yes? let elevenOrLonger = IF tdensity EQ 1 THEN 1
yes? LIST elevenOrLonger[L=@SUM] ! if answer is zero, there are no spells this long
define a new time series in which you mask out the blocks of contiguous 1s just found
Iterate the procedure to find number and location of spells of shorter and shorter duration. (Note only odd-number length spells can be found this way.)

Or ... as Hein says, you can write raw code to perform any custom calculation you want by using either Python in pyFerret or defining your own external function in vanilla Ferret. ==> Then contribute the idea back to the Ferret user's list! A filter that, at each point, measured the length of contiguous block of non-missing values surrounding it, would be a nice, generic Ferret tool. (Should we add such a filter into Ferret? It would be a fairly easy addition: e.g. LET contiguous_spell_length = mytimeseries[l=@cntg])

I cant help wondering if a clever use of the nearest index above and nearest index below filters might be used to build a contiguous length filter using a few LET definitions ...

- Steve

I guess @evnt transformation would be helpful for this, in particular defining ID number as outlined here http://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2013/msg00277.html would provide an solution. I made an attempt accordingly but haven't been able to solve the problem.

I would appreciate any lead on this.

Thanks,

Jagadish