L\c @seddlZddlZddlmZmZmZddlZddl Z ddl m Z ddl mZddlmZddlmZgggggdZgggggdZgggggd Zggggd Zggggd Zd d e jjd dZedZdddddddded ZdddZdZ dS(iN(tDatasettnum2datetdate2num(tpyplot(tBasemap(tPolygon(tPatchCollectioncCsct|}|jd}|jd}|jd} |jd} |jd} t|dkr|\} } }}||k||k@|| k@|| k@}||}||}| |} | |} | dd|f} n|jd}t|d|j}|jtjj| }t j t j |j }t j |d |d | }t|dkr||d|d !}nt|dkr'|\}}}|||!}tj tj|jjd||jdk}|jdd|f}| |} | |} ||}||}nt|dkr|\}}}|||!}tj tj|jjd||jdk}|jdd|f}| |} | |} ||}||}nt|dkrPxm|D]b}tj |j|k}|jdd|f}| |} | |} ||}||}qWn|| | ||fS( Ntlattlontidtnametprittimetunitstindextcolumnsi(Rt variablestlenRR tclosetnptmatmasked_invalidtpdt to_datetimetdatet DataFrametinverttisnantvaluestsumtshapetiloctkeys(tncfilet_bndst _timebndst_cov1t_cov2t_excludetncdsett_latst_lonst_idst_namestdatatx1tx2ty1ty2t_selt_timest_datest_prt _obsdatestfytlytfracttestt_nn((sfunctions_rainfall.pytreadSAWSonefile sb *       5     5        c' Csd}t|}|jd}|jd}|jd} |jd} |jd} t|dkr|\} } }}||k||k@|| k@|| k@}||}||}| |} | |} | dd|f} n|jd}t|d |j}|jtjj| }t j d d d d }t j |d|d| }t ddH}|j d} tjg| D]}|jjd^q} WdQXtj| dddf}| dddf}| dddf}d}x|D]}| dddf|k}t jg||D]}|j^q=}||}d||dk<|jd}|dkrd}n|dkrt j |d|d|g}n0t j|t j |d|d|ggd}|d}q W|jj|_t j ddd d } tjt| |jdg}!tj|!(|jt j |!d| d|j}|j|t|dkr||d|d!}nt|dkr_|\}"}#}$||"|#!}%tjtj|%jjd|$|%jdk}|j dd|f}| |} | |} ||}||}nt|dkr|\}"}#}$||"|#!}%tjtj|%jjd|$|%jdk}|j dd|f}| |} | |} ||}||}nt|dkrxm|D]b}&tj|j|&k}|j dd|f}| |} | |} ||}||}qWn|| | ||fS(Ns'/work/data/saws/south_africa_2015.pr.nctlatitudet longitudeR R R iR R s 1850-01-01s 2015-12-31tfreqtDRRs%../data/bigsix/sawsupdate/alldata.csvtrit,iitnanttfloatsVILLIERSDORP-SOS ARSt VILLIERSDORPs 2016-01-01s 2017-12-31(!RRRRR RRRRRt date_rangeRtopent readlinestarraytstriptsplittuniqueRtastypetconcatRtzerosRRBtappendR tupdateRRRRR('R"R#R$R%R&R!R'R(R)R*R+R,R-R.R/R0R1R2R3R4R5tinpftxt_unamest_udatest_uprtnt_nametseltseldatesttemptprupdatet appenddatestfillerR6R7R8R9R:((sfunctions_rainfall.pytreadSAWSGs *      7 ,    !0" *  5     5        cCst|}|jd}|jd}|jd} |jd} |jd} t|dkr|\} } }}||k||k@|| k@|| k@}||}||}| |}| |}| dd|f} ntj| }tj| }|jd}t|d|j}|jtjj | }tj ||dk|6d }7d,t||7t?|7d+}/tj@tj|tjA|$d t|d+}8dd|d-d }9|jB|9t|d }:|:|8|/};|$|;}4|$|;}3i| d.6|#d/6|"d06|!d16|$d26|3d36|4d46}<|<S(5scalculates bootstrap confidence interval and significance level for trend, ignoring autocorrelation or accounting for it Parameters ---------- x_var : list independent variable y_var : list dependent variable, same length as x_var q : int, optional, only if trendtype==quantreg quantile for which regression is to be calculated n : int, optional number of bootstrap samples ci : int, optional confidence level. Default is for 95% confidence interval frac : int, optional, only if trendtype==lowess lowess parameter (fraction of time period length used in local regression) it : int, optional, only if trendtype==lowess lowess parameter (numbre of iterations) autocorr : str, optional way of accounting for autocorrelation, possible values: None, "bootstrap" trendtype : str, optional method of trend derivation, possible values: lowess, linreg, quantreg, TheilSen CItype : str, optional method of CI derivation, possible values: "analytical" and "bootstrap". if trendtype is "lowess", CItype will be set to None if CItype is "analytical": autocorrelation will be set to None Results ------- returns library with following elements: slope - slope of the trend CI_high - CI on the slope value CI_low - as above pvalue - trend's significance level trend - trend line, or rather its y values for all x_var trendCI_high - confidence interval for each value of y trendCI_low - as above Remarks ------- the fit function ocassionally crashes on resampled data. The workaround is to use try statement iN(tOLS(tQuantReg(tECDF(tmstatst mannwhitneyuttt kendalltau(tStationaryBootstrapt IIDBootstrapt analyticaltTheilSenRidiitquantregsQuantile regression, CI type: s, autocorrelation adjustment: s RtXtYsY ~ XtqiiRsLinear regression, CI type: thasconsttmissingtdropsTheil-Sen slope, CI type: s raw MK tau:sraw MK pvalue:Rgffffff?tlowesssLowess R8titRDRii_g?g?gY@tslopetCI_hightCI_lowtpvaluettrendt trendCI_lowt trendCI_high(Ctstatsmodels.apitapit#statsmodels.regression.linear_modelRtstatsmodels.tsattsatstatsmodels.formula.apitformulat*statsmodels.regression.quantile_regressionRtstatsmodels.nonparametric.apit nonparametrict0statsmodels.distributions.empirical_distributionRt scipy.statsRRRRtarch.bootstrapRRRRIRRRRtstrRRt column_stackRtfittparamst InterceptRtpvaluestconf_inttpredictRRR`treshapeRtixRPt add_constantRnt theilslopestdictRtmaptgetRMt interpolateRRBt nanpercentilet percentileRt transposetmeanRRtnanmeantppf(=tx_varty_vartn_boottcit trendtypeRR8RRtCItypeRtsmRRtsmfRtnpsmRRRRRRRtn_datatci_lowtci_hightbstxydatatmodeltrestintcptRRRRty_predtresidst bs_trendst bs_slopest bs_intcptsR,tindtgoingdownquietlytx_varOLSttaut interceptR[tpredt temp_interptconfintRRtecdftxtemptsquanomtseetci_adjttcomptcinttoutput((sfunctions_rainfall.pyttrend_CIs,"    '! *     %.      2 . $    ; *020            !"7  7Ri c!CsVddlm}m}ddlm}m}ddlm}d}t|} | |d} t j | | g} t j | (t j | | g} t j | (xt || dD]} | GH| d}| |}|dt j t | }xZt ||D]I}|| d}|| }|||!}t j| }|d kr}|j||d d \}}}}|||\}}ng}xV|j|D]E}|d d }|j||||d d }||d g}qWt j|}||} | d }|d krd|}n|| | ||f<|| | ||fs"  =\=HF! $,