# [ncl-talk] how can I avoid dividing by zero?

Sri Nandini snandini at marum.de
Wed Aug 22 06:40:33 MDT 2018

```

Hello

Thank you for your quick reply.
Im not sure how to limit numbers reasonably above this.
Attached is the map i plotted.

And below my entire original script. would be grateful for any suggestions.

f= addfile("SLP_PHIS_Z3_PI.nc", "r") ;
T41    = f->Z3(:,{850},:,:)

printVarSummary(T41)     ; [Time|12]x[ilev | 26] x [lat | 96] x [lon | 144]

T41 at _FillValue = -9.96921e+36
;print(T41)

aveX    = dim_avg_n_Wrap(T41,0)
printVarSummary(aveX)                      ; (lat,lon)

aveX at _FillValue =9.96921e+36
;print(where(aveX.ne.0, aveX, aveX at _FillValue))

aveX = 1. / where(aveX.ne.0, aveX, aveX at _FillValue)
; print(aveX)
;==============================================================

f1= addfile("totalwinds_PI.nc", "r") ;
T411    = f1->U(:,{850},:,:)

printVarSummary(T411)       ; [Time|12]x[ilev | 26] x [lat | 96] x [lon | 144]
T411 at _FillValue = -9.96921e+36

aveX1    = dim_avg_n_Wrap(T411,0)
printVarSummary(aveX1)                                ; (lat,lon)

aveX1 at _FillValue = 9.96921e+36

aveX1 = 1. / where(aveX1.ne.0, aveX1, aveX1 at _FillValue)
;print(aveX1)

;==============================================================
f2= addfile("th_PI.nc", "r") ;
T412    = f2->TH(:,{900},:,:)

printVarSummary(T412)       ; [Time|12]x[ilev | 27] x [lat | 96] x [lon | 144]
T412 at _FillValue = -9.96921e+36

aveX2    = dim_avg_n_Wrap(T412,0)
printVarSummary(aveX2)                                ; (lat,lon)

aveX2 at _FillValue = 9.96921e+36

aveX2 = 1. / where(aveX2.ne.0, aveX2, aveX2 at _FillValue)
;print(aveX2)

;==============================================================
;    Read latitudes
;    The 'eady_growth_rate' function requires that 'lat' and 'th' agree
;    Use 'conform' the propogate the lat values
;==============================================================

xlat = f->lat                               ;  [lat | 96]
printVarSummary(xlat)

XLAT = conform ( aveX2 , xlat , 0 )
printVarSummary(XLAT)

copy_VarMeta(aveX2,XLAT)
printVarSummary(XLAT)

XLAT = 1. / where(XLAT.ne.0, XLAT, XLAT at _FillValue)
;print(XLAT)

;==============================================================

egr = eady_growth_rate(aveX2, aveX1, aveX, XLAT, 0,  1)
printVarSummary(egr)

printMinMax(egr,0)

;print(egr)
print("-----------------------------------------")

;==============================================================

wks = gsn_open_wks("pdf","Eady_850")              ; send graphics to PNG file

;---Set some basic plot options

res               = True
res at gsnMaximize   = True       ; maximize plot in frame
;res at gsnAddCyclic  = False

res at cnFillOn      = True
res at cnLinesOn     = False
;res at cnFillMode    = "RasterFill"                 ; slow here
res at cnFillMode    = "CellFill"                   ; faster

res at tiMainString         = "Baroclinicity at 850hPa"
res at gsnCenterString      = " "
res at gsnLeftString        = " "
res at gsnRightString       = " "

;res at cnFillPalette        = "precip2_17lev"    ;BrRiEgr

;res at mpMaxLatF                   = 80           ; choose a different subregion
;res at mpMinLatF                   = 30
;res at mpMaxLonF                   = 70
;res at mpMinLonF                   = -40
res at mpFillOn      = False

res at lbBoxEndCapStyle  = "TriangleBothEnds"
res at lbTitleString       = "Maximum Eady growth rate per day"
res at lbTitlePosition     = "bottom"              ; title position
res at lbTitleFontHeightF  = .012                ; make title smaller
;res at pmLabelBarOrthogonalPosF =-0.01              ;-- move the labelbar upward
;res at lbLabelFontHeightF = 0.01

;==============================================================
;Eady growth rate (1/day)

egr        = egr*86400
egr at units  = "annual 1/day"

res at cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
res at cnMinLevelValF       =  -2               ; set min contour level
res at cnMaxLevelValF       =  2               ; set max contour level
res at cnLevelSpacingF      =  0.2              ; set contour spacing

res at gsnRightString  =  egr at units
contour = gsn_csm_contour_map(wks, egr(:,:),res)

and this is the output: there is no error but upon printing i can see the gridpoints which are nan and some which have fillvalues assigned.
Variable: T41
Type: float
Total Size: 663552 bytes
165888 values
Number of Dimensions: 3
Dimensions and sizes:    [time | 12] x [lat | 96] x [lon | 144]
Coordinates:
time: [136358.5..136692.5]
lat: [ -90..  90]
lon: [   0..357.5]
Number Of Attributes: 5
lev :    867.1607600000013
mdims :    1
units :    m
long_name :    Geopotential Height (above sea level)
cell_methods :    time: mean time: mean

Variable: aveX
Type: float
Total Size: 55296 bytes
13824 values
Number of Dimensions: 2
Dimensions and sizes:    [lat | 96] x [lon | 144]
Coordinates:
lat: [ -90..  90]
lon: [   0..357.5]
Number Of Attributes: 7
_FillValue :    -9.96921e+36
lev :    867.1607600000013
mdims :    1
units :    m
long_name :    Geopotential Height (above sea level)
cell_methods :    time: mean time: mean
average_op_ncl :    dim_avg_n over dimension(s): time

Variable: T411
Type: float
Total Size: 663552 bytes
165888 values
Number of Dimensions: 3
Dimensions and sizes:    [time | 12] x [lat | 96] x [lon | 144]
Coordinates:
time: [136358.5..136692.5]
lat: [ -90..  90]
lon: [   0..357.5]
Number Of Attributes: 5
lev :    867.1607600000013
mdims :    1
units :    m/s
long_name :    Zonal wind
cell_methods :    time: mean time: mean

Variable: aveX1
Type: float
Total Size: 55296 bytes
13824 values
Number of Dimensions: 2
Dimensions and sizes:    [lat | 96] x [lon | 144]
Coordinates:
lat: [ -90..  90]
lon: [   0..357.5]
Number Of Attributes: 7
_FillValue :    -9.96921e+36
lev :    867.1607600000013
mdims :    1
units :    m/s
long_name :    Zonal wind
cell_methods :    time: mean time: mean
average_op_ncl :    dim_avg_n over dimension(s): time

Variable: T412
Type: float
Total Size: 663552 bytes
165888 values
Number of Dimensions: 3
Dimensions and sizes:    [time | 12] x [lat | 96] x [lon | 144]
Coordinates:
time: [136358.5..136692.5]
lat: [ -90..  90]
lon: [   0..357.5]
Number Of Attributes: 5
ilev :    903.3002900000016
mdims :    2
units :    K
long_name :    Potential Temperature
cell_methods :    time: mean time: mean

Variable: aveX2
Type: float
Total Size: 55296 bytes
13824 values
Number of Dimensions: 2
Dimensions and sizes:    [lat | 96] x [lon | 144]
Coordinates:
lat: [ -90..  90]
lon: [   0..357.5]
Number Of Attributes: 7
_FillValue :    -9.96921e+36
ilev :    903.3002900000016
mdims :    2
units :    K
long_name :    Potential Temperature
cell_methods :    time: mean time: mean
average_op_ncl :    dim_avg_n over dimension(s): time

Variable: xlat
Type: double
Total Size: 768 bytes
96 values
Number of Dimensions: 1
Dimensions and sizes:    [lat | 96]
Coordinates:
lat: [ -90..  90]
Number Of Attributes: 2
long_name :    latitude
units :    degrees_north

Variable: XLAT
Type: double
Total Size: 110592 bytes
13824 values
Number of Dimensions: 2
Dimensions and sizes:    [96] x [144]
Coordinates:

Variable: XLAT
Type: double
Total Size: 110592 bytes
13824 values
Number of Dimensions: 2
Dimensions and sizes:    [lat | 96] x [lon | 144]
Coordinates:
lat: [ -90..  90]
lon: [   0..357.5]
Number Of Attributes: 8
ilev :    903.3002900000016
mdims :    2
units :    K
long_name :    Potential Temperature
cell_methods :    time: mean time: mean
average_op_ncl :    dim_avg_n over dimension(s): time
_FillValue_original :    9.96921e+36
_FillValue :    9.969209968386869e+36

Variable: egr
Type: double
Total Size: 110592 bytes
13824 values
Number of Dimensions: 2
Dimensions and sizes:    [lat | 96] x [lon | 144]
Coordinates:
lat: [ -90..  90]
lon: [   0..357.5]
Number Of Attributes: 3
_FillValue :    9.969209968386869e+36
long_name :    maximum eady growth rate
units :
(0)    maximum eady growth rate : min=nan   max=nan
(0)    -----------------------------------------

On Aug 22, 2018 2:33:35 PM, Barry Lynn wrote:
> Hi:
> I am not sure of all the steps, but small numbers are not NaNs.  Perhaps you are dividing by close to zero, when you should limit to a number reasonably above?
>
> Is that possible to do?
>
>
> On Wed, Aug 22, 2018 at 2:40 PM Sri Nandini <> snandini at marum.de> > wrote:

> > Yes.
> >
> > 1) I use geopotential height and the values were originally like this: 1080.752
> >
> >    T41    = f->Z3(:,{850},:,:)
> >
> >    printVarSummary(T41)                        ; [Time|12]x [lat | 96] x [lon | 144]
> >    T41 at _FillValue = -9.96921e+36
> >    print(T41)

> > 2) Then i performed: dim_avg_n_Wrap and values were like: 0.0009252814
> >
> >     aveX    = dim_avg_n_Wrap(T41,0)
> >    printVarSummary(aveX)                      ; (lat,lon)
> >    aveX at _FillValue =9.96921e+36
> >    print(where(aveX.ne.0, aveX, aveX at _FillValue))
> >  values like this:  1080.752
> >
> >    aveX = 1. / where(aveX.ne.0, aveX, aveX at _FillValue)
> >    print(aveX)
> > 0.0009252814
> >
> > 3) I use this in the eady growth rate equation and that gives me nan values for few gridpoints.
> > > >    egr = eady_growth_rate(aveX2, aveX1, aveX, XLAT, 0,  1)
> >    printVarSummary(egr)
> >
> >    printMinMax(egr,0)
> >
> > Variable: egr
> > Type: double
> > Total Size: 110592 bytes
> >             13824 values
> > Number of Dimensions: 2
> > Dimensions and sizes:    [lat | 96] x [lon | 144]
> > Coordinates:
> >             lat: [ -90..  90]
> >             lon: [   0..357.5]
> > Number Of Attributes: 3
> >   _FillValue :    9.969209968386869e+36
> >   long_name :    maximum eady growth rate
> >   units :
> > (0)    maximum eady growth rate : min=nan   max=nan

> > 4) All the inputs into the eady growth rate equation have data points but some are almost zeros, as such i assume it gets masked out to nan?

> > On Aug 22, 2018 12:15:37 PM, Barry Lynn wrote:

> > > > > Hi:
> > > Are the points that actually have data producing NaNs from the calculations?
> > >
> > > Barry
> > >
> > >
> > > On Wed, Aug 22, 2018 at 12:56 PM Sri Nandini <> > > snandini at marum.de> > > > wrote:

> > > > Thank you
> > > >
> > > > I used it and have seen that i have data points except some gridpoints which point to nan and the fillvalue both.
> > > > > > > > I think then there is no data over these points.

> > > > On Aug 22, 2018 10:46:56 AM, Barry Lynn wrote:

> > > > > > > > > Hello Sri:
> > > > > I would try it this way.
> > > > >
> > > > > Set > > > > > y> > > > >  to a _FillValue value where it is equal to 0, and then can do the divide:
> > > > > > > > > >    y> > > > > @_FillValue = 1.E36 ; or something like this.
> > > > > > > > > >    print(where(y.ne.0,y,> > > > > @_FillValue)) ;  I am assuming this works.
> > > > > > > > > > If nothing prints then maybe you have no values.
> > > > > > > > > > In any case, try the next two lines.
> > > > > > > > > >    y    = > > > > > where> > > > > (y.ne.0,y,y at _FillValue)
> > > > >   yinv = 1. / y> > > > > However, just to check, I would

> > > > >
> > > > > On Wed, Aug 22, 2018 at 11:39 AM Sri Nandini <> > > > > snandini at marum.de> > > > > > wrote:
> > > > > > Hello
> > > > > >
> > > > > >
how can I remove the zero values to avoid dividing by zero the error I got is as the following:
> > > > > >
> > > > > >
fatal:divide: Division by 0, Can't continue
> > > > > >
fatal:Div: operator failed, can't continue
> > > > > >
> > > > > >
I have tried using   yinv = 1. / where(y.ne.0, y, y at _FillValue) but now it gives me all nan values
> > > > > >
> > > > > >
Heres the bit of my code
> > > > > >
> > > > > >
T41    = f->Z3(:,{850},:,:)
> > > > > >
> > > > > >
printVarSummary(T41)
> > > > > >
T41 at _FillValue = -9.96921e+36
> > > > > >
> > > > > >
aveX    = dim_avg_n_Wrap(T41,0)
> > > > > >
printVarSummary(aveX)                      ; (lat,lon)
> > > > > >
T41 at _FillValue = -9.96921e+36
> > > > > >
> > > > > >
aveX = 1. / where(aveX.ne.0, aveX, aveX at _FillValue)
> > > > > >
> > > > > >
aveX1 and aveX2 below are of same type as aveX.
> > > > > >
XLAT = conform ( aveX2 , xlat , 0 )
> > > > > >
printVarSummary(XLAT)
> > > > > >
> > > > > >
egr = eady_growth_rate(aveX2, aveX1, aveX, XLAT, 0,  1)
> > > > > >
> > > > > >
Output
> > > > > >
===========================
> > > > > >
Variable: egr
> > > > > >
Type: double
> > > > > >
Total Size: 110592 bytes
> > > > > >
13824 values
> > > > > >
Number of Dimensions: 2
> > > > > >
Dimensions and sizes:   [lat | 96] x [lon | 144]
> > > > > >
Coordinates:
> > > > > >
lat: [ -90..  90]
> > > > > >
lon: [   0..357.5]
> > > > > >
Number Of Attributes: 3
> > > > > >
_FillValue :  -9.969209968386869e+36
> > > > > >
long_name :   maximum eady growth rate
> > > > > >
units :
> > > > > >
(0)     maximum eady growth rate : min=nan   max=nan
> > > > > >
> > > > > >
> > > > > >
Any advice is much appreciated
> > > > > >
_______________________________________________
> > > > > >
ncl-talk mailing list
> > > > > > ncl-talk at ucar.edu
> > > > > >
List instructions, subscriber options, unsubscribe:
> > > > > > http://mailman.ucar.edu/mailman/listinfo/ncl-talk

> > > > >

> > > > > --
> > > > > Barry H. Lynn, Ph.D> > > > > Senior Associate Scientist, Lecturer,
> > > > > > > > > > The Institute of the Earth Science,
> > > > > The Hebrew University of Jerusalem,
> > > > > Givat Ram, Jerusalem 91904, Israel
> > > > > Tel: 972 547 231 170
> > > > > Fax: (972)-25662581
> > > > >

> > > > > > > > > > C.E.O, Weather It Is, LTD
> > > > > Weather and Climate Focus
> > > > > http://weather-it-is.com
> > > > > Jerusalem, Israel
> > > > > Local: 02 930 9525
> > > > > Cell: 054 7 231 170
> > > > > Int-IS: x972 2 930 9525
> > > > >

> > > > >
> > > > >
> > > > >
> > > > > > > > > >

> > > >

> > >

> > > --
> > > Barry H. Lynn, Ph.D> > > Senior Associate Scientist, Lecturer,
> > > > > > The Institute of the Earth Science,
> > > The Hebrew University of Jerusalem,
> > > Givat Ram, Jerusalem 91904, Israel
> > > Tel: 972 547 231 170
> > > Fax: (972)-25662581
> > >

> > > > > > C.E.O, Weather It Is, LTD
> > > Weather and Climate Focus
> > > http://weather-it-is.com
> > > Jerusalem, Israel
> > > Local: 02 930 9525
> > > Cell: 054 7 231 170
> > > Int-IS: x972 2 930 9525
> > >

> > >
> > >
> > >
> > > > > >

> >

>

> --
> Barry H. Lynn, Ph.D> Senior Associate Scientist, Lecturer,
> > The Institute of the Earth Science,
> The Hebrew University of Jerusalem,
> Givat Ram, Jerusalem 91904, Israel
> Tel: 972 547 231 170
> Fax: (972)-25662581
>

> > C.E.O, Weather It Is, LTD
> Weather and Climate Focus
> http://weather-it-is.com
> Jerusalem, Israel
> Local: 02 930 9525
> Cell: 054 7 231 170
> Int-IS: x972 2 930 9525
>

>
>
>
> >

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180822/eb052562/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Eady_850.pdf
Type: application/pdf
Size: 164846 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180822/eb052562/attachment-0001.pdf>
```

More information about the ncl-talk mailing list