Re: handing missing values in gamma fit

From: Sanjiv Kumar <skumar.water_at_nyahnyahspammersnyahnyah>
Date: Thu Oct 18 2012 - 13:26:44 MDT

Hi Dennis:

Thank you for the reply. Below are some clarification:

On Thu, Oct 18, 2012 at 12:35 PM, Dennis Shea <shea@ucar.edu> wrote:

> I think that there may be grid points where a climatology
> was used. Hence, all Januaries were constant.
>

>>> I think you are right.

>
> =====
> The following does not address the problem but will make
> make all bad points equall _FillValue
>
> Not sure why you did the following:
>
> var = f1->pre(stiR:endiR, :, :)
>

>>> to subset the original data for the selected time period (1950 to
2004). Original data is for 1901 to 2009.

> var = log10(var + 0.01) ; comment this out?
>

>>>>> precipitation data generally do not follow the Gaussian distribution.
Taking a log transformation generally bring the non-nomral distribution to
a normal distribution.

>>>>> adding 0.01 is to avoid zero values, because log(0) is not defined.
original precip data unit is mm/month. Hence the values (0.01) is very
small.

>
> Replace
>
> do nl = 0, nlat-1
> do ml = 0, nlon-1
> if (all(.not.ismissing(var(:, nl, ml))))then
> do k = 0, ds_len-1
> var_spi(k, :, nl, ml) = (/dim_spi_n(var(:,nl ml), dr_len(k),
> False, 0)/)
> end do
> end if
> end do
>

>>>> I did this little tweaking to overcome the missing data issue.
Otherwise the script was stoping without calculating SPI for non-missing
land grid points.

>
> with
>
> wcStrt = systemfunc("date") ; timing only
>
> do k = 0, ds_len-1
> var_spi(k, :, :, :) = (/dim_spi_n(var, dr_len(k), False, 0)/)
> end do
>
> This is *MUCH* faster than manually looping. 8 seconds on CGD systems/
>

>>>> I agree, this is much faster. I implemented this where i was sure
there is not missing grid points (e.g. many CMIP5 models), it took around
10 minutes instead of 10+ hours.

>
> Add
>
> nNaN = num(isnan_ieee(var_spi))
> print(nNaN)
> if (nNaN.gt.0) then
> var_spi = where(isnan_ieee(var_spi), var_spi@_FillValue, var_spi)
> end if
>
> nBig = num(var_spi.gt.1e10)
> print(nBig)
> if (nBig.gt.0) then
> var_spi = where(var_spi.gt.1e30, var_spi@_FillValue, var_spi)
> end if
>
> wallClockElapseTime(wcStrt,"**LOOP",0)
>
> >>> this piece of code is helpful. I think this can take care of my
problem at least for now.

>
>
>
>
>
> On 10/18/2012 11:52 AM, Mary Haley wrote:
>
>> This will be fixed, but unfortunately not in the next version (6.1.0) of
>> NCL.
>>
>> If you really need a fix, we can get you a new binary soon.
>>
>> --Mary
>>
>> On Oct 17, 2012, at 10:28 AM, Sanjiv Kumar wrote:
>>
>> Hi Mary:
>>>
>>> thank you for the reply. You may also want to fix another problem with
>>> 'dim_spi_n'
>>>
>>> i get following error message:
>>>
>>> >>>>>>>>>>>>>>>
>>>
>>> Error in anvnrmd(). Prob. not in [0,1.0]
>>>
>>> >>>>>>>>>>
>>>
>>> I the output file i get some non-wanted values e.g.
>>> "9.99999993381581e+36" and "NaN".
>>>
>>> I am guessing this issue is related to the above error message.
>>>
>>> Let me know if you have any question.
>>>
>>> -Sanjiv Kumar
>>>
>>>
>>>
>>> On Wed, Oct 17, 2012 at 8:22 AM, Mary Haley <haley@ucar.edu
>>> <mailto:haley@ucar.edu>> wrote:
>>>
>>> This will be fixed in V6.1.0, which we are trying to get out the
>>> door this week or next.
>>>
>>> --Mary
>>>
>>> On Oct 15, 2012, at 2:50 PM, Sanjiv Kumar wrote:
>>>
>>> > Hi,
>>> >
>>> > I am trying to use 'dim_spi_n' in NCL 6.1.0. IT works great if i
>>> do not have missing data. But when data array contains missing
>>> data i get following error message:
>>> >
>>> > >>>>>>
>>> > Error in gamfitd - empty data array
>>> > >>>>>
>>> >
>>> > For example, i have a global land only precipitation so all
>>> ocean data point are empty.
>>> > Since, the program encounter first missing ocean data point, it
>>> gives the above error and stops. It does not calculate the SPI
>>> values for the land points.
>>> >
>>> > Is there some way to work around this issue.
>>> >
>>> > Thank you in advance.
>>> >
>>> > -sanjiv
>>> >
>>> > ______________________________**_________________
>>> > ncl-talk mailing list
>>> > List instructions, subscriber options, unsubscribe:
>>> > http://mailman.ucar.edu/**mailman/listinfo/ncl-talk<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<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<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
Received on Thu Oct 18 13:27:01 2012

This archive was generated by hypermail 2.1.8 : Tue Oct 23 2012 - 11:10:04 MDT