Re: Bug in where function?

I did this Monday evening and, inadvertently, did not post
==========================================================
If 'a' and 'b' are one dimensional

         c = a ; same dimensionality as 'a'
         c = 0.0 ; array syntax (init to 0.0)
         ii = ind(a.gt.0.0)
         c(ii) = b(ii)/a(ii)

----
A slight variation on Kyle's suggestion is the following.
Note: This approach has a bug which has been fixed in
NCL v6.2.0 so it can not be used thru 6.1.2
    c = 1.0/where(a .gt. 0, b/a, a@_FillValue)
This assumes 'a' has an _FillValue
'c' will have _FillValue in all instances where a.le.0
  To get 0s, you would have to do the following on the next lines
          c@_FillValue = 0.0    ; all _FillValue set to 0.0
          delete(c@_FillValue)
A possible unfortunate side effect would be cases where 'a'
had_FillValue prior to the calculation. Then it would
be impossible to distinguish between the original
suite of a@_FillValue and the above values.
---
quick and dirty 'solution'
undef("my_zero_where")
  function my_zero_where(a,b)  ; a can be any dimensionality
                               ; rank[a] = rank[b]
  local dima, dimb, ranka, rankb, c, a1d, b1d, c1d, ii, jj, iir, jjr, n\
      , dimiir, dimjjr, niir, njjr
  begin
     dima  = dimsizes(a)
     ranka = dimsizes(dima)
     dimb  = dimsizes(b)
     rankb = dimsizes(dimb)
     if (ranka.ne.rankb) then
         print("my_zero_where: a,b do not conform")
         exit
     end if
     if (all(a.gt.0.0)) then
                            ; nothing special need be done
         c = b/a            ; array syntax ... 'a' any dimensionality
     else
         c     = (/ a /)           ; create array; same a@_FillValue
         if (ranka.eq.1) then
             ii    = ind(a.gt.0.0)
             c(ii) = b(ii)/a(ii)
             jj    = ind(a.le.0.0)
             c(jj) = 0.0
         else
             a1d   = ndtooned(a)
             b1d   = ndtooned(b)
             c1d   = (/ a1d /)
             ii      = ind(a1d.gt.0.0)
             c1d(ii) = b1d(ii)/a1d(ii)
             jj      = ind(a1d.le.0.0)
             c1d(jj) = 0.0
                   ; 'resolve' [remap] 1d indices to to nD
             iir   = ind_resolve(ii,dima)
             dimiir= dimsizes(iir)
             niir  = dimiir(0)       ; number of elements a1d > 0.0
             if (ranka.ne.3) then
                 print("my_zero_where: the nD part supports 3D only")
                 print("               The current array is "+ranka+"D")
                 exit
             end if
             do n=0,niir-1
                c(iir(n,0),iir(n,1),iir(n,2)) = c1d(ii(n))
             end do
             jjr   = ind_resolve(jj,dima)
             dimjjr= dimsizes(jjr)   ; number of elements a1d <= 0.0
             njjr  = dimjjr(0)       ; number of elements a1d > 0.0
             do n=0,njjr-1
                c(jjr(n,0),jjr(n,1),jjr(n,2)) = c1d(jj(n))
             end do
         end if
     end if
     return(c)
  end
;++++++++++++++++++++ TEST DRIVER ++++++++++++++++++++
;                         MAIN
  N   = 100
  A1D = random_normal(0,5,N)
  B1D = random_uniform( -5, 5, N)
  A1D@_FillValue = 1e20
  IX  = ispan(0,N-1,5)
  A1D(IX) = A1D@_FillValue
  I0  = ispan(3,N-1,5)
  A1D(I0) = 0.0
  C1D =  my_zero_where(A1D,B1D)
  print(C1D)
On 1/20/14, 7:06 AM, Robert Schuster wrote:
> Thank you! Working with missing values is a good idea!
> It is good to know, that the expressions for both cases are evaluated for all values in the arrays and not only for those where the condition is true.
>
> Robert
>
> Am 20.01.2014 um 14:41 schrieb Kyle Griffin <ksgriffin2@wisc.edu>:
>
>> Unfortunately, Robert, that is exactly how the where function is defined to work. I suggest turning the zeros into missing values (at least temporarily) for the calculation and then subbing the values from the original function back in with another where statement, but others may have better suggestions to work around this issue.
>>
>>
>> Kyle
>>
>> ----------------------------------------
>> Kyle S. Griffin
>> Department of Atmospheric and Oceanic Sciences
>> University of Wisconsin - Madison
>> Room 1421
>> 1225 W Dayton St, Madison, WI 53706
>> Email: ksgriffin2@wisc.edu
>>
>>
>> On Mon, Jan 20, 2014 at 3:05 AM, Robert Schuster <rxschuster@gmail.com> wrote:
>> Hello,
>>
>> I have two arrays that I want to divide by each other, one of them contains zeros:
>>
>> a = (/0,1/)
>> b = (/1,1/)
>>
>> I tried to use the where function to avoid a divide by zero error:
>>
>> c = where(a .gt. 0, b/a, 0)
>>
>> But interestingly the divide operation is performed on the whole arrays and not only where the condition a .gt. 0 is true:
>>
>> fatal:divide: Division by 0, Can't continue
>> fatal:Div: operator failed, can't continue
>>
>> Is that the way the where function is supposed to work?
>>
>> Robert
>> _______________________________________________
>> 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
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk