[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