import numpy as np import pandas as pd from netCDF4 import Dataset, num2date, date2num import os, glob from matplotlib import pyplot as plt from mpl_toolkits.basemap import Basemap from matplotlib.patches import Polygon from matplotlib.collections import PatchCollection def readSAWSonefile(ncfile, _bnds=[], _timebnds=[], _cov1=[], _cov2=[], _exclude=[]): # print ncfile, _bnds, _timebnds, _cov1, _cov2, _exclude ncdset=Dataset(ncfile) #print ncdset.variables _lats=ncdset.variables['lat'][:] _lons=ncdset.variables['lon'][:] _ids=ncdset.variables['id'][:] _names=ncdset.variables['name'][:] data=ncdset.variables['pr'][:] if len(_bnds)>0: x1,x2,y1,y2=_bnds _sel=(_latsy1) & (_lons>x1) & (_lons0: _pr=_pr[_timebnds[0]:_timebnds[1]] if len(_cov1)>0: fy,ly,frac=_cov1 test=_pr[fy:ly] _sel=(np.invert(np.isnan(test.values)).sum(0)>frac*test.shape[0]) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] if len(_cov2)>0: fy,ly,frac=_cov2 test=_pr[fy:ly] _sel=(np.invert(np.isnan(test.values)).sum(0)>frac*test.shape[0]) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] if len(_exclude)>0: for _nn in _exclude: _sel=np.invert(_pr.keys()==_nn) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] return _pr,_ids,_names,_lats,_lons def readSAWS(_bnds=[], _timebnds=[], _cov1=[], _cov2=[], _exclude=[]): ncfile="/work/data/saws/south_africa_2015.pr.nc" ncdset=Dataset(ncfile) #print ncdset.variables _lats=ncdset.variables['latitude'][:] _lons=ncdset.variables['longitude'][:] _ids=ncdset.variables['id'][:] _names=ncdset.variables['name'][:] data=ncdset.variables['pr'][:] if len(_bnds)>0: x1,x2,y1,y2=_bnds _sel=(_latsy1) & (_lons>x1) & (_lons0: _pr=_pr[_timebnds[0]:_timebnds[1]] if len(_cov1)>0: fy,ly,frac=_cov1 test=_pr[fy:ly] _sel=(np.invert(np.isnan(test.values)).sum(0)>frac*test.shape[0]) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] if len(_cov2)>0: fy,ly,frac=_cov2 test=_pr[fy:ly] _sel=(np.invert(np.isnan(test.values)).sum(0)>frac*test.shape[0]) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] if len(_exclude)>0: for _nn in _exclude: _sel=np.invert(_pr.keys()==_nn) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] return _pr,_ids,_names,_lats,_lons def readDWS2(_ncfile, _bnds=[], _timebnds=[], _cov1=[], _cov2=[],_exclude=[]): ncdset=Dataset(_ncfile) _lats=ncdset.variables['latitude'][:] _lons=ncdset.variables['longitude'][:] _ids0=ncdset.variables['id'][:] _names0=ncdset.variables['name'][:] data=ncdset.variables['pr'][:] if len(_bnds)>0: x1,x2,y1,y2=_bnds sel=(_latsy1) & (_lons>x1) & (_lons0: _pr=_pr[_timebnds[0]:_timebnds[1]] if len(_cov1)>0: fy,ly,frac=_cov1 test=_pr[fy:ly] _sel=(np.invert(np.isnan(test.values)).sum(0)>frac*test.shape[0]) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] if len(_cov2)>0: fy,ly,frac=_cov2 test=_pr[fy:ly] _sel=(np.invert(np.isnan(test.values)).sum(0)>frac*test.shape[0]) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] if len(_exclude)>0: for _nn in _exclude: _sel=np.invert(_pr.keys()==_nn) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] return _pr,_ids,_names,_lats,_lons def readDWS_old(_bnds=[], _timebnds=[], _cov=[], _exclude=[]): ncfile="/work/data/dws/pr_dws_day_19000101-20171208.nc" ncfile="/work/data/dws/pr_dws_day_19000101-20181231.nc" ncdset=Dataset(ncfile) #print ncdset.variables _lats=ncdset.variables['latitude'][:] _lons=ncdset.variables['longitude'][:] _ids0=ncdset.variables['id'][:] _names0=ncdset.variables['name'][:] data=ncdset.variables['pr'][:] if len(_bnds)>0: x1,x2,y1,y2=_bnds sel=(_latsy1) & (_lons>x1) & (_lons0: _pr=_pr[_timebnds[0]:_timebnds[1]] if len(_cov)>0: fy,ly,frac=_cov test=_pr[fy:ly] _sel=(np.invert(np.isnan(test.values)).sum(0)>.9*test.shape[0]) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] if len(_exclude)>0: for _nn in _exclude: _sel=np.invert(_pr.keys()==_nn) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] return _pr,_ids,_names,_lats,_lons def readCoCT(_bnds=[], _timebnds=[], _cov=[], _exclude=[]): ncfile="/work/data/dws/pr_dws_day_19000101-20171208.nc" ncdset=Dataset(ncfile) #print ncdset.variables _lats=ncdset.variables['latitude'][:] _lons=ncdset.variables['longitude'][:] _ids0=ncdset.variables['id'][:] _names0=ncdset.variables['name'][:] data=ncdset.variables['pr'][:] if len(_bnds)>0: x1,x2,y1,y2=_bnds sel=(_latsy1) & (_lons>x1) & (_lons0: _pr=_pr[_timebnds[0]:_timebnds[1]] if len(_cov)>0: fy,ly,frac=_cov test=_pr[fy:ly] _sel=(np.invert(np.isnan(test.values)).sum(0)>.9*test.shape[0]) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] if len(_exclude)>0: for _nn in _exclude: _sel=np.invert(_pr.keys()==_nn) _pr=_pr.iloc[:,_sel] _ids=_ids[_sel] _names=_names[_sel] _lats=_lats[_sel] _lons=_lons[_sel] return _pr,_ids,_names,_lats,_lons def plotmap(pl, _lats,_lons, ext, color="blue", edgecolor="blue", cmap=plt.cm.RdYlBu_r, msize=4): m=Basemap(projection='merc',resolution='i', lat_0=-30,lon_0=10, llcrnrlat=ext[2],urcrnrlat=ext[3], llcrnrlon=ext[0], urcrnrlon=ext[1], epsg = 4269) m.drawparallels(np.arange(-50.,0.,0.5), labels=[1,0,0,0], dashes=[1,1], linewidth=0.25, color='0.5') m.drawmeridians(np.arange(-10, 100., 0.5), labels=[0,0,1,0], dashes=[1,1], linewidth=0.25, color='0.5') m.drawmapboundary(fill_color='white', color="gray") # m.shadedrelief() m.fillcontinents(color='0.8', lake_color='white', zorder=0) # m.drawcountries(color='0.6', linewidth=1) # m.readshapefile("/work/data/gis/southafrica/bigsix2d", "dams") try: m.readshapefile("/work/data/gis/southafrica/bigsix2d", "dams") except: m.readshapefile("/work/orange/work/data/gis/southafrica/bigsix2d", "dams") # for info, shape in zip(m.dams_info, m.dams): # if info['NAME']=="Nigeria": # x, y = zip(*shape) # m.plot(x, y, marker=None, color='grey', linewidth=1) patches = [] for info, shape in zip(m.dams_info, m.dams): # if info['nombre'] == 'Selva': patches.append( Polygon(np.array(shape), True) ) pl.add_collection(PatchCollection(patches, facecolor= 'deepskyblue', edgecolor='deepskyblue', linewidths=1., zorder=2)) m.scatter(_lons,_lats, color=color, edgecolors=edgecolor, latlon=True, linewidth=3, s=msize, zorder=5, alpha=0.5) # m.arcgisimage(service = 'World_Physical_Map', xpixels = 500, verbose= False) # m.arcgisimage(service = 'World_Shaded_Relief', xpixels = 500, verbose= False) return m def autocorr_test(_xdata, _ydata, verbose=False): from statsmodels.stats.diagnostic import acorr_ljungbox from statsmodels.tsa.stattools import acf #all statst need regularly spaced, continuous time series - just y variable #Durbin-Watson statistics: # calculated correctly with missing data # but no significance level. Apparently critical values for DW are not implemented in any python library #ACF: # crashes on missing data # Ljung-Box: # crashes on missing data too _ydata=np.ma.masked_invalid(_ydata) #autocorrelation in residuals #this is acf function that does not allow nans # print "\nautocorrelation for first three lags:", acf(_ydata)[1:4] #this is from pandas, is nan agnostic pdf=pd.Series(_ydata, index=_xdata, copy=True) if verbose: print "autocorrelation for first three lags:", [pdf.autocorr(i) for i in range(1,4)] #durbin-watson a=_ydata[:-1].astype('float') b=_ydata[1:].astype('float') _stat=np.nansum((b-a)**2)/np.nansum(_ydata**2) if verbose: print "Durbin-Watson statistic (close to 2 if no autocorrelation):", _stat _stat, _pvalue=acorr_ljungbox(_ydata, lags=1, boxpierce=False) if verbose: print "Ljung-Box p-value on lag 1 autocorrelation:", _pvalue print "" return [pdf.autocorr(i) for i in range(1,4)] def trend_CI(x_var, y_var, n_boot=1000, ci=95, trendtype="linreg", q=0.5, frac=0.6, it=3, autocorr=None, CItype="bootstrap", verbose=False): """calculates 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 """ #for linreg import statsmodels.api as sm from statsmodels.regression.linear_model import OLS #for arima import statsmodels.tsa as tsa #for quantreg import statsmodels.formula.api as smf from statsmodels.regression.quantile_regression import QuantReg #for lowess import statsmodels.nonparametric.api as npsm #other from statsmodels.distributions.empirical_distribution import ECDF from scipy.stats import mstats, mannwhitneyu, t, kendalltau from arch.bootstrap import StationaryBootstrap, IIDBootstrap #preparing data if CItype=="analytical" and trendtype=="TheilSen": CItype="bootstrap" x_var=np.array(x_var) y_var=np.ma.masked_invalid(y_var) n_data=len(y_var) ci_low=(100-ci)/2 ci_high=100-ci_low #setting bootstrapping function if autocorr=="bootstrap": bs=StationaryBootstrap(3, np.array(range(len(y_var)))) else: bs=IIDBootstrap(np.array(range(len(y_var)))) if trendtype=="quantreg": if verbose: print "Quantile regression, CI type: "+CItype+", autocorrelation adjustment: "+str(autocorr)+"\n" xydata=pd.DataFrame(np.column_stack([x_var, y_var]), columns=['X', 'Y']) model=smf.quantreg('Y ~ X', xydata) res=model.fit(q=q) intcpt=res.params.Intercept slope=res.params.X pvalue=res.pvalues[1] CI_low=res.conf_int()[0]['X'] CI_high=res.conf_int()[1]['X'] y_pred=res.predict(xydata) #calculating residuals resids=y_var-y_pred #calculate autocorrelation indices autocorr_test(x_var, resids, False) if CItype=="bootstrap": #bootstrapping bs_trends=np.copy(y_pred).reshape(-1,1) bs_slopes=[] bs_intcpts=[] for data in bs.bootstrap(n_boot): ind=data[0][0] model = smf.quantreg('Y ~ X', xydata.ix[ind,:]) try: res = model.fit(q=q) bs_slopes=bs_slopes+[res.params.X] bs_intcpts=bs_intcpts+[res.params.Intercept] bs_trends=np.append(bs_trends,res.predict(xydata).reshape(-1,1), 1) except: goingdownquietly=1 if trendtype=="linreg": if verbose: print "Linear regression, CI type: "+CItype+", autocorrelation adjustment: "+str(autocorr)+"\n" x_varOLS = sm.add_constant(x_var) model = sm.OLS(y_var, x_varOLS, hasconst=True, missing='drop') res = model.fit() intcpt,slope=res.params pvalue=res.pvalues[1] CI_low,CI_high=res.conf_int()[1] y_pred=res.predict(x_varOLS) #calculating residuals resids=y_var-y_pred #calculate autocorrelation indices autocorr_test(x_var, resids) if CItype=="bootstrap": #bootstrapping for confidence intervals bs_slopes=[] bs_intcpts=[] bs_trends=np.copy(y_pred).reshape(-1,1) for data in bs.bootstrap(n_boot): ind=data[0][0] model = sm.OLS(y_var[ind], x_varOLS[ind,:], hasconst=True, missing='drop') try: res = model.fit() bs_slopes=bs_slopes+[res.params[1]] bs_intcpts=bs_intcpts+[res.params[0]] bs_trends=np.append(bs_trends,res.predict(x_varOLS).reshape(-1,1), 1) except: goingdownquietly=1 if trendtype=="TheilSen": if verbose: print "Theil-Sen slope, CI type: "+CItype+", autocorrelation adjustment: "+str(autocorr)+"\n" #significance of MK tau tau,pvalue=kendalltau(x_var, y_var) if verbose: print "raw MK tau:", tau, "raw MK pvalue:", pvalue #TS slope and confidence intervals slope,intercept,CI_low,CI_high=mstats.theilslopes(y_var, x_var, alpha=0.95) #getting slope line's y values y_pred=intercept+slope*x_var #calculating residuals resids=y_var-y_pred #calculate autocorrelation indices autocorr_test(x_var, resids) if CItype=="bootstrap": #bootstrapping for confidence intervals bs_slopes=[] bs_intcpts=[] bs_trends=np.copy(y_pred).reshape(-1,1) for data in bs.bootstrap(n_boot): ind=data[0][0] res=mstats.theilslopes(y_var[ind], x_var[ind], alpha=0.95) bs_slopes=bs_slopes+[res[0]] bs_intcpts=bs_intcpts+[res[1]] bs_trends=np.append(bs_trends, (res[1]+res[0]*x_var).reshape(-1,1), 1) if trendtype=="lowess": if verbose: print "Lowess\n" temp=dict(npsm.lowess(y_var, x_var, frac=frac, it=it, missing="drop")) y_pred=np.array(map(temp.get, x_var)).astype("float").reshape(-1,1) bs_trends=np.copy(y_pred) for data in bs.bootstrap(n_boot): ind=data[0][0] try: temp = dict(npsm.lowess(y_var[ind], x_var[ind], frac=frac, it=it, missing="drop")) temp=np.array(map(temp.get, x_var)).astype("float").reshape(-1,1) pred=pd.DataFrame(temp, index=x_var) temp_interp=pred.interpolate().values bs_trends=np.append(bs_trends, temp_interp, 1) except: goingdownquietly=1 #calculating final values of CI and p-value #skipping when lowess if trendtype=="lowess": CI_low=np.nan CI_high=np.nan slope=np.nan intcpt=np.nan pvalue=np.nan confint=np.nanpercentile(bs_trends, [ci_low,ci_high], 1) trendCI_low=confint[:,0] trendCI_high=confint[:,1] else: if CItype=="bootstrap": #values for slope, intercept and trend can be obtained as medians of bootstrap distributions, but normally analytical parameters are used instead # it the bootstrap bias (difference between analytical values and bootstap median) is strong, it might be better to use bootstrap values. # These three lines would need to be uncommented then # slope=np.median(bs_slopes) # intcpt=np.median(bs_intcpts) # trend=intcpt+slope*x_var #these are from bootstrap too, but needs to be used for this accounts for autocorrelation, which is the point of this script CI_low,CI_high=np.percentile(bs_slopes, [5, 95]) ecdf=ECDF(bs_slopes) pvalue=ecdf(0) #this makes sure we are calculating p-value on the correct side of the distribution. That will be one-sided pvalue if pvalue>0.5: pvalue=1-pvalue if verbose: print bs_trends.shape confint=np.nanpercentile(bs_trends, [ci_low,ci_high], 1).transpose() trendCI_low=confint[:,0] trendCI_high=confint[:,1] # print confint.shape else: #this is for analytical calculation of trend confidence interval #it happens in the same way for each of the trend types, thus it is done here, not under the trendtype subroutines #making sure x are floats xtemp=np.array(x_var)*1.0 #squared anomaly squanom=(xtemp-np.mean(xtemp))**2 temp=((1./len(x_var))+(squanom/sum(squanom)))**0.5 #standard error of estmation see=(np.nansum((np.array(y_var)-np.nanmean(y_pred))**2)/len(x_var))**0.5 #adjusting ci ci_adj=1-((1-ci/100.)/2) #accounting for uncertainty in mean through student's t tcomp=t.ppf(ci_adj, len(x_var)-2) #confidence interval cint=tcomp*see*temp #for trend only trendCI_high=y_pred+cint trendCI_low=y_pred-cint # print cint # print trendtype, "slope:",slope, "pvalue (one sided):", pvalue, "conf interval:", CI_low, CI_high, "autocorrelation adjustment:", autocorr, "\n" output={"slope":slope, "CI_high":CI_high, "CI_low":CI_low, "pvalue":pvalue, "trend": y_pred, "trendCI_low":trendCI_low, "trendCI_high":trendCI_high} return output def trendpyramid(_ts, _how="analytical", _minwindow=10): # _ts is a time series of annual values # minwindow - is the minimum window (number of years) for which trends are calculated from scipy.stats import mstats, kendalltau from arch.bootstrap import StationaryBootstrap, IIDBootstrap from statsmodels.distributions.empirical_distribution import ECDF _nboot=100 _nts=len(_ts) _nlevels=_nts-_minwindow+1 # number of "levels" _slope=np.zeros([_nlevels, _nts]) _slope[:]=np.nan _pval=np.zeros([_nlevels, _nts]) _pval[:]=np.nan for _window in range(_minwindow, _nts+1): print _window _first=_window/2 _last=_nts-_first bs=StationaryBootstrap(3, np.array(range(_window))) for _block in range(_first, _last): _x0=_block-_window/2 _x1=_x0+_window _subts=_ts[_x0:_x1] x_var=np.arange(_window) if _how=="analytical": _trend,intercept,CI_low,CI_high=mstats.theilslopes(_subts, x_var, alpha=0.95) _tau,_pvalue=kendalltau(x_var, _subts) else: bs_slopes=[] for data in bs.bootstrap(_nboot): ind=data[0][0] res=mstats.theilslopes(_subts[ind], x_var[ind], alpha=0.95) bs_slopes=bs_slopes+[res[0]] _trend=np.median(bs_slopes) ecdf=ECDF(bs_slopes) _pvalue=ecdf(0) #this makes sure we are calculating p-value on the correct side of the distribution. That will be one-sided pvalue if _pvalue>0.5: _pvalue=1-_pvalue _slope[_window-_minwindow,_block]=_trend _pval[_window-_minwindow,_block]=_pvalue return [_slope, _pval] def plotpyramid(_slope, _pval, plotparams): fy,ly,vmin,vmax=plotparams ny=ly-fy+1 fig=plt.figure(figsize=(10,10)) pl=plt.subplot(1,1,1) ax=pl.imshow(np.flipud(_slope*10), interpolation="nearest", cmap=plt.cm.bwr_r, vmin=vmin, vmax=vmax) #m.colorbar() pl.contourf(np.flipud(_pval)<0.05, 1, hatches=['', '///'], alpha=0) pl.set_xticks(range(5,ny,10)) pl.set_xticklabels(range(fy+5,ly,10), rotation=75) pl.set_yticks(range(0,ny-10,10)) pl.set_yticklabels(range(ny,10,-10)) pl.set_aspect(0.65) pl.set_xlabel("centered on year") pl.set_ylabel("window width [years]") plt.colorbar(ax, fraction=0.15, shrink=0.4) plt.title("rainfall trend [mm/decade]") plt.show()