[ncl-talk] missing values in reg_multlin_stats
Dennis Shea
shea at ucar.edu
Fri Jan 26 12:33:52 MST 2018
Hi Anahita,
THX for ftp'ing the stability file.
I've attached two simple scripts.
(a) reg_simple.ncl_1:
Create area weighted (30N-90N) time series of 'CLDLOW' and
'theta_ANN'. These are input to the 'reg_multlin_stats' function.
(b) reg_simple.ncl_:
Same as (a) but masks out all non-land locations.
So, input to the reg_multlin_stats' function are just two time series with
no missing data.
For fun: I have added a very basic 2-dimensional PDF to each.
===
Re: Your offline question about your use of the ttest within
'SDE_significance.ncl'
That is a rather large script so I skimmed throgh it. I think you use of
'ttest' is correct.
===
A word of caution:
The 4-dimensional model data are on hybrid levels: (time,lev,lat,lon)
The 'omega.nc' variable is on isobaric levels: (time,lev,lat,lon)
Just to be explicitly clear about the vertical coordinate, even though the
dimension names are the same the are very different variables.
Hybrid levels are a function of surface pressure at each grid point. At a
particular location:
*p(k) = a(k)*p0 + b(k)*psfc*
The 'lev' coordinate variable for the model variables is a 'nominal' level.
It is the pressure level **if** the surface pressure was 1000 hPa. It is
just a placeholder. Think if psfc was say 732 hPa.
===
Good Luck
D
On Wed, Jan 24, 2018 at 10:44 PM, Dennis Shea <shea at ucar.edu> wrote:
> Hello,
>
> You did ftp most of the files. However, the 'stability_E_total.nc' file is
> not present.
>
> ---
> Furher, offline you sent two scripts. One was a "simple script"; the other
> was (let's say) not so simple'.
>
>
> ====
> [1]
> "my code is very simple:"
>
> I edited your "simple code" to make it as simple as possible
> ; **edited
> code**
> ana=addfile("clouds_E_total.nc","r")
> CLDLOW=ana->CLDLOW ; (180,96,144)
> printVarSummary(CLDLOW)
>
> ana=addfile("stability_E_total.nc","r") ; missing file ... not on ftp
> site
> theta_ANN=ana->theta_ANN
> printVarSummary(theta_ANN) ; (???, ?96?,?144?) <=== missing
> file
>
> CLDLOW_1d=ndtooned(CLDLOW)
> theta_ANN_1d=ndtooned(theta_ANN) ; is this the same shape & size as
> CLDLOW????
>
> print("nmsg_CLDLOW="+num(ismissing(CLDLOW_1d))) ; =0
> print("nmsg_theta_ANN="+num(ismissing(theta_ANN_1d)))
>
> ; the following inputs 180*96*180 values for each variable.
> ; Is that what you want to do???
>
> opt = True
> opt at print_anova = True ; optional
> opt at print_data = True
> b = reg_multlin_stats(CLDLOW_1d,theta_ANN_1d,opt) ; partial regression
> coef
>
>
> Not sure what to say. Does variable theta_ANN have 180 time steps?
> The _ANN (?Annual?) suffix would imply to me that these are annual values
> rather than monthly values.
> If this is the case you would have to make the CLDLOW variables also be
> Annual values.
> The dimensionalities must match.
>
> You said that the stability has missing values. The dimensionalites MUST
> match . If so,
> you can use the following two step approacj
>
> (a)
> CLDLOW at _FillValue = theta_ANN at _FillValue
> CLDLOW = where(ismissing(theta_ANN), CLDLOW at _FillValue, CLDLOW)
>
> (b)
> The way you have the independent and dependent variables as 1D would
> weight grid points near the poles the same as values at low latitudes. I do
> *not* think this is appropriate. High latitude grid points represent much
> smaller areas.
>
> After step (a) you should compute weighted areal averages
>
> http://www.ncl.ucar.edu/Document/Functions/Contributed/wgt_areaave_Wrap.
> shtml
>
> The files you sent do not have any gaussian weights [ gwt ] so you could
> the cos(rad*lat)
> or
> http://www.ncl.ucar.edu/Document/Functions/Contributed/latGauWgt.shtml
>
> for the wgty.
>
> The 'wgt_areaave' function allows missing values.
>
> This will yield one global time series which is what I think you want.
>
> =====================
>
> Your 2nd offline script is quite a bit more complicated. Also, it looks
> like you only want latitudes north of 30N. Correct? I'd suggest the
> following approach. It will make for much smaller arrays.
>
> diri = "./"
> diri = "/project/cas/shea/netCDF_Misc/"
>
> latStrt = 30.0
> latLast = 90.0
>
> ana=addfile(diri+"clouds_E_total.nc","r") ; *hybrid* levels
> CLDLOW=ana->CLDLOW(:,{latStrt:latLast},:)
> printVarSummary(CLDLOW) ; [time | 180] x [lat | 32]
> x [lon | 144]
>
> print(CLDLOW&lat)
>
> ---
> Your 2nd script is treating variable RH's hybrid levels to be mixed with
> omega's isobaric levels. In my opinion, this is not appropriate.
>
> This script also mask out ocean areas.
>
> You can use an approach like (b)
>
> lsm = .... ; land-sea mask
>
> Then, use an approach where CLDLOW , theta_ANN are masked via (say)
>
> http://www.ncl.ucar.edu/Document/Functions/Built-in/mask.shtml
>
> Then do the global averges; then the regression.
>
> ----
> Please when you send code to ncl-talk, it shoul be clean and structured.
> Please do not include numerous commented lines. It is tedious to read code
> that has been heavily commented.
>
> Good luck
> D
>
>
>
>
>
> On Wed, Jan 10, 2018 at 7:55 PM, Dennis Shea <shea at ucar.edu> wrote:
>
>> The "reg_multlin_stats" documentation
>> http://www.ncl.ucar.edu/Document/Functions/Contributed/reg_
>> multlin_stats.shtml
>> does state that but the code (located in contributed.ncl) has:
>>
>> undef ("reg_multlin_stats")
>> function reg_multlin_stats(Y[*]:numeric, XP:numeric, opt) ; should be
>> float or double
>> ;
>> ; Nomenclature:
>> ; Y - dependent variable (size NY)
>> ; missing values (_FillValue) are not allowed. <=============
>> ; XP - one [*] or more [*][*] independent variables.
>> ; Missing values (_FillValue) are not allowed. <=============
>>
>> [snip]
>>
>> if (isatt(Y,"_FillValue") .and. any(ismissing(Y))) then
>> print("reg_multlin_stats: Y has missing values. Not allowed!")
>> exit
>> end if
>>
>> if (isatt(XP,"_FillValue") .and. any(ismissing(XP))) then
>> print("reg_multlin_stats: XP has missing values. Not allowed!")
>> exit
>> end if
>> [snip]
>> ==========
>>
>> Still, there might be a way to use the function. My thought is to
>> manually loop ove ?grid points? skipping any points that have missing
>> values.
>>
>> How are you using it? Please be clear and concise.
>>
>> I am gone through 22 Jan 2018 so maybe someone else can help.
>>
>>
>>
>>
>>
>>
>>
>>
>>
>> On Wed, Jan 10, 2018 at 4:00 PM, Anahita Amiri Farahani via ncl-talk <
>> ncl-talk at ucar.edu> wrote:
>>
>>> Dear all,
>>>
>>> I wanted to use reg_multlin_stats to calculate linear regression between
>>> low cloud and stability, the low cloud and stability are my dependent and
>>> independent variables. Stability is calculated based on the potential
>>> temperature at the surface and 700 hPa. So we have missing values for some
>>> areas. When I read the description of this function, it says While
>>> missing values are allowed, it is recommended that users not input any
>>> missing independent values. It just confuses the results. Missing values
>>> should be indicated by the _FillValue
>>> <https://www.ncl.ucar.edu/Document/Language/fillval.shtml> attribute.
>>>
>>>
>>> But when I use this function I got this error: reg_multlin_stats: XP
>>> has missing values. Not allowed!
>>>
>>> and it popped out from NCL, Any suggestions?
>>>
>>> Thanks,
>>> Ana
>>>
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> ncl-talk at ucar.edu
>>> 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/20180126/913ad77c/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: reg_simple.ncl_1
Type: application/octet-stream
Size: 2624 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180126/913ad77c/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: reg_simple.ncl_2
Type: application/octet-stream
Size: 2994 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180126/913ad77c/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pdf_cld_theta_1.png
Type: image/png
Size: 29895 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180126/913ad77c/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pdf_cld_theta_2.png
Type: image/png
Size: 31255 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180126/913ad77c/attachment-0001.png>
More information about the ncl-talk
mailing list