# [ncl-talk] trend_manken data size limit?

Mon Apr 6 16:54:04 MDT 2015

```Hi Dennis, Mary et al.
thanks for the run down.

The data are TRMM precip rate (0.25 deg, every 3 hours) from 1998 to
2010.  The single grid point
is over southwest Mexico City (Rio Magdalena Basin).  It basically only
rains in the afternoon/evening during June, July, August, September, hence
all of the 0.0 values.

I have used regline and this is what I will use for my results (a
significant decrease in precip rate over the 13 year period, but,
surprisingly, a significant increase in flowrate (Rio Magdalena Basin).

saludos,
Dave

On Mon, Apr 6, 2015 at 5:42 PM, Dennis Shea <shea at ucar.edu> wrote:

> You sent your NCL script and test data offline.
>
> The time series that it fails on has 29737 values... ~84% (24878) are 0.0
> The is a really unusual distribution.
>
> NCL's  Mann-Kendall code is a fortran subroutine which does two
> [1] Compute the Mann-Kendall (MK) trend statistic; [2] Compute the
> Theil-Sen (TS) robust estimate of  linear trend.
>
> [1] I am sure the MK statistic and probability estimate were computed.
> [2] The TS estimate requires that the 'slopes' calulated within the do
> loop  be sorted into ascending order. The code fails in the fcode used for
> sorting. Specifically, at the following location of DPSORTU
>
> [SNIP]
>
>              65 I = I + 1
>                 IF (I.EQ.J) GO TO 55
>  1341 ----> T = A(I+1)
>                 IF (A(I).LE.T) GO TO 65
>                 K = I
> [SNIP]
>
> The subroutine used is a 'famous' one translated from ALGOL (original
> author: ?Wilkinson?) to fortran (Singleton). It is part of the ACM codes
> and is available via netlib.
>
> C Sorts A into ascending order, from A(1) to A(N).
> C Reference: Richard C. Singleton, Algorithm 347, SORT. <=======
> C Comm. ACM 3, 321 (March 1969).
> C Algorithm is Copyright 1969 Association of Computing Machinery,
>
>
> I speculate that the preponderance of 0.0 values has 'somehow' causes this
> code to fail. Why it required > 25,000 values to fail .. I do not know.
> Again, I speculate, if there had been more 0.0 in the 1st 25,000 it may
> have failed with fewer points.
>
> There is no chance that I am going to change this sorting code. It
> requires more knowledge about sorting than I have. NOTE: the code used by
> NCL's 'dim_pqsort' is based on this code so I doubt this would be an
>
> ----
> OK .... MK is 'distribution free' BUT I think some knowledge of how MK
> works should be used. Basically, it checks to see if successive values go
> up (+1), down (-1) or are tied (0). The *amount *of the change is not used
> for MK. With series like 0,0,0,#,0,0,0,0,0,#,#,0,0,0,0,0,0,#,0,0,#,0,0,...
> etc. It will often go from 0.0 to some positive number (+1) and from some
> positive number to 0.0 (-1).  I would speculate that the +1 and -1 counts
> almost nullify each other. All cases where the series has 0,0 would be a
> tie. This will be *by far* the largest number. The algorithm uses a
> correction for ties but I am not sure if was developed for cases where the
> number of ties is >> greater that the sum of the +1 and -1.
>
> Also, some knowledge of hoe Theil-Sen estimates are made. In the process
> of calculating  the sums and ties  'slopea' are calculated. The median of
> the slopes is the TS estimate. Since the majority of slopes are 0, the TS
> median value will be 0.0.
>
> Just using regline ...
>
>      rc    = regline(year_fraction,precip_ave)          ; slope
>      print(rc)
>
>
> Variable: rc     <==== ALL
> Type: float
> Total Size: 4 bytes
>             1 values
> Number of Dimensions: 1
> Dimensions and sizes:    [1]
> Coordinates:
> Number Of Attributes: 8
>   units :    mmhr^-1/3hrs
>   _FillValue :    -999.99
>   yintercept :    2.898926
>   yave :    0.112996
>   xave :    2006.162
>   nptxy :    29737
>   rstd :    0.000660348
>   tval :    -2.102962
> (0)    -0.001388687
>
>
>
>   ii = ind(precip_ave.gt.0.0)  ; observations of when it rains
>   RC = regline(year_fraction(ii),precip_ave(ii))
>   print(RC)
>
> Variable: RC  <=== Conditioned on rain only times.
> Type: float
> Total Size: 4 bytes
>             1 values
> Number of Dimensions: 1
> Dimensions and sizes:    [1]
> Coordinates:
> Number Of Attributes: 7
>   _FillValue :    -999.99
>   yintercept :    -6.783739
>   yave :    0.6915334
>   xave :    2005.787
>   nptxy :    4859
>   rstd :    0.003920668
>   tval :    0.9505657
> (0)    0.003726852
>
>
>
>
>
>
> On Sun, Apr 5, 2015 at 9:04 AM, Mary Haley <haley at ucar.edu> wrote:
>
>> Dave,
>>
>> Could you write the precip_ave variable to a NetCDF file and send us the
>> file, so we can try trend_manken on it here?
>>
>> You can do something like:
>>
>> fout->pcp = precip_ave
>>
>> and then if the file is too large (say over 25 megabytes), then you can
>> use our ftp to upload it. See
>>
>> http://www.ncl.ucar.edu/report_bug.shtml#HowToFTP
>>
>> Thanks,
>>
>> --Mary
>>
>>
>> On Sat, Apr 4, 2015 at 12:56 PM, David Adams <dave.k.adams at gmail.com>
>> wrote:
>>
>>> Hi all,
>>> I am using the trend_manken
>>> function to look at the significance of  precipitation intensity trends
>>> (using TRMM data) for 15 years (every 3 hours data).  When I use 10,000
>>> data points it takes about 3 seconds to run, 20,000 data points about 10
>>> seconds, 25,000 about 15 seconds.  When I use the entire time series 29,737
>>> data points, it never finishes.  I have waited more than 20 minutes and it
>>> keeps spinning.  I switched from double to floating point, but that makes
>>> no difference.
>>>
>>> ;---------------------------------
>>> ; Mann Kendall test
>>> ;---------------------------------
>>> pa = trend_manken(precip_ave, False, 0)
>>> print(pa)
>>>
>>>
>>>
>>> Any ideas?
>>>
>>> saludos,
>>> Dave
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>>
>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150406/9142d2ee/attachment.html
```