[ncl-talk] missing values in reg_multlin_stats
Dennis Shea
shea at ucar.edu
Wed Jan 24 22:44:06 MST 2018
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/20180124/67f023a5/attachment.html>
More information about the ncl-talk
mailing list