<jupyter_start><jupyter_text>Importing Data<jupyter_code># -*- coding: utf-8 -*-
"""
Created on  March 2021

@author: <NAME>
"""

import pandas as pd
import glob
import matplotlib.pyplot as plt
import numpy as np
all_files=glob.glob('*.txt')
li=[]
for i in all_files:
    df = pd.read_csv(i,sep=";",skiprows=1)
    j=df.drop(["Deutschland","Unnamed: 19"],axis=1)
    column=list(j.columns[2:])
    f=pd.melt(df,id_vars=["Monat","Jahr"],value_vars=column,var_name="States",value_name=str(j.loc[1][1]))
    li.append(f)

    
from functools import reduce
#df_merged = reduce(lambda  left,right: pd.merge(left,right,how='outer'),li)
data=pd.DataFrame({"YEAR":li[0]["Jahr"],"State":li[1]["States"],"JAN":li[0]["1.0"],
                    "FEB":li[1]["2.0"],"MAR":li[2]["3.0"],"APR":li[3]["4.0"],"MAY":li[4]["5.0"],
                   "JUN":li[5]["6.0"],"JUL":li[6]["7.0"],"AUG":li[7]["8.0"],"SEP":li[8]["9.0"],
                    "OCT":li[9]["10.0"],"NOV":li[10]["11.0"],"DEC":li[11]["12.0"]})<jupyter_output><empty_output><jupyter_text>adding seasonal and annual columns<jupyter_code>data["DEC-FEB"]=data.iloc[:,2:4].sum(axis=1)
data["MAR-MAY"]=data.iloc[:,4:7].sum(axis=1)
data["JUN-AUG"]=data.iloc[:,7:11].sum(axis=1)
data["SEP-NOV"]=data.iloc[:,11:14].sum(axis=1)
data["ANNUAL"]=data.iloc[:,2:14].sum(axis=1)
data
#data= data[(data['YEAR'] > 1959) & (data['YEAR'] <= 2020)].reset_index(drop=True)
data<jupyter_output><empty_output><jupyter_text>line plot of annual rainfall<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import pymannkendall as mk
def STATES(column,division):
        Year = []
        annual = []
        for i in range(len(column)):
            if column[i] == division:
                if 'nan'!= str(data["ANNUAL"][i]) and 'NAN' != str(data["ANNUAL"][i]):
                    Year.append(data["YEAR"][i])
                    annual.append((data["ANNUAL"][i]))
        return Year,annual
colors = ["#808000","#FF0000","#0000FF","#808000","#800080","#008000","#800000","#A52A2A","#FFA500","#000000","#151B54","#FBB917","#806517","#C11B17",
          "#810541","#F6358A","#808000","#FF0000","#0000FF","#808000","#800080","#008000","#800000","#A52A2A","#FFA500","#000000","#151B54","#FBB917",
          "#806517","#C11B17","#810541","#F6358A","#808000","#FF0000","#0000FF","#808000"]
country = data["State"].unique()
name=list(data["State"].unique())
for j in range(len(country)):
    Year, annual = STATES(data["State"],str(country[j]))
    res = mk.original_test(annual,alpha=.05)
    print(res)
    trend_line = np.arange(len(annual)) * res.slope + res.intercept
    plt.style.use('fivethirtyeight')
    plt.rcParams['figure.figsize']=(25,5)
    fix,ax=plt.subplots(1,1)
    ax.plot(Year, annual, color=colors[j] ,alpha = 0.7, marker='o', linestyle='dashed',linewidth=2, markersize=12)
    ax.plot(Year, trend_line)
    plt.title(str(country[j])+" Annual Rainfall\n", size= 25)
    plt.xlabel("Year", size= 25)
    plt.ylabel("Rainfall in\nmm", size= 30)
    plt.tight_layout()
    x=country[j].replace('/', '_')
   # plt.savefig("plots/"+x+".jpg",dpi=150)
    plt.show()
from cycler import cycler
color_c = cycler('color', ['k'])
style_c = cycler('linestyle', ['-', '--', ':', '-.'])
markr_c = cycler('marker', ['', '.', 'o'])
c_cms = color_c * markr_c * style_c
c_csm = color_c * style_c * markr_c


def mm2inch(*tupl):
    inch = 25.4
    if isinstance(tupl[0], tuple):
        return tuple(i/inch for i in tupl[0])
    else:
        return tuple(i/inch for i in tupl)
def mm2point(mm):
    return mm/(25.4/72)
font = {'family' : 'Times New Roman',
         'weight' : 'normal',
         'size'   : 10}
plt.rc('font', **font)
plt.rcParams['axes.linewidth'] = mm2point(0.0001)
plt.rcParams['ytick.major.width'] = mm2point(0.01)
plt.rcParams['xtick.major.width'] = mm2point(0.01)<jupyter_output><empty_output><jupyter_text>box plot<jupyter_code>%matplotlib
import seaborn as sns
font = {'family' : 'Times New Roman',
         'weight' : 'normal',
         'size'   : 10}
plt.rc('font', **font)
plt.rcParams['figure.figsize']=(16,9)
label=['Brandenburg/\nBerlin','Brandenburg','Baden-\nWuerttemberg',
       'Bayern','Hessen','Mecklenburg-\nVorpommern',
       'Niedersachsen','Niedersachsen/\nHamburg/\nBremen',
       'Nordrhein-\nWestfalen','Rheinland-\nPfalz','Schleswig-Holstein','Saarland',
       'Sachsen','Sachsen-\nAnhalt','Thueringen/\nSachsen-Anhalt','Thueringen']

m1=data.groupby("State",sort=False)["ANNUAL"].mean().round(2)
#vertical_offset = data['ANNUAL'].median() * 0.5
ax = sns.boxplot(x="State", y="ANNUAL", data=data,color="#989898",width=0.6,linewidth=2,showmeans=False,
                meanprops={"marker":"o",
                       "markerfacecolor":"white", 
                       "markeredgecolor":"black",
                      "markersize":"9"})
ax.set_xlabel('State',fontsize=10,alpha=0.8,family="Times New Roman")
ax.set_ylabel('Annual precipitation (in mm)',fontsize=11,alpha=0.8, family="Times New Roman")
ax.set_xticklabels(label,family="Times New Roman")
#plt.title('Annual Rainfall in States in Germany',fontsize=40,)
ax.tick_params(axis='x',rotation=0,grid_alpha=0.8)
ax.tick_params(axis='y',rotation=0,grid_alpha=0.8)
#axis='x', which='major', bottom=False, top=False

plt.xticks(alpha=0.8)
plt.yticks(alpha=0.8)

for xtick in ax.get_xticks():
    ax.text(xtick,m1[xtick]-1,m1[xtick],size='x-small',horizontalalignment='center',color='w', weight='semibold',family="Times New Roman")
plt.tight_layout()
plt.show()
#plt.savefig("plots/boxplot.jpg",dpi=300)<jupyter_output><empty_output><jupyter_text>Trend Table<jupyter_code>#3trend_test table
cols=["State","trend","h","p","z"]
trend_data=pd.DataFrame(columns=cols,index=(range(16)))
for k in range(len(country)):
    Year, annual = STATES(data["State"],str(country[k]))
    res = mk.hamed_rao_modification_test(annual,alpha=.05)
    trend_data["State"][k]=country[k]
    trend_data["trend"][k]=res.trend
    trend_data["h"][k]=res.h
    trend_data["p"][k]=res.p
    trend_data["z"][k]=res.z



#fall trend
def SUMMER(column,division):
        Year = []
        summer = []
        for i in range(len(column)):
            if column[i] == division:
                if 'nan'!= str(data["SEP-NOV"][i]) and 'NAN' != str(data["SEP-NOV"][i]):
                    Year.append(data["YEAR"][i])
                    summer.append((data["SEP-NOV"][i]))
        return Year,summer
cols=["State","trend","h","p","z"]
trend_data_f=pd.DataFrame(columns=cols,index=(range(16)))
for k in range(len(country)):
    Year, summer = SUMMER(data["State"],str(country[k]))
    res = mk.hamed_rao_modification_test(summer,alpha=.05)
    trend_data_f["State"][k]=country[k]
    trend_data_f["trend"][k]=res.trend
    trend_data_f["h"][k]=res.h
    trend_data_f["p"][k]=res.p
    trend_data_f["z"][k]=res.z

#summer trend
def SUMMER(column,division):
        Year = []
        summer = []
        for i in range(len(column)):
            if column[i] == division:
                if 'nan'!= str(data["JUN-AUG"][i]) and 'NAN' != str(data["JUN-AUG"][i]):
                    Year.append(data["YEAR"][i])
                    summer.append((data["JUN-AUG"][i]))
        return Year,summer
cols=["State","trend","h","p","z"]
trend_data_s=pd.DataFrame(columns=cols,index=(range(16)))
for k in range(len(country)):
    Year, summer = SUMMER(data["State"],str(country[k]))
    res = mk.hamed_rao_modification_test(summer,alpha=.05)
    trend_data_s["State"][k]=country[k]
    trend_data_s["trend"][k]=res.trend
    trend_data_s["h"][k]=res.h
    trend_data_s["p"][k]=res.p
    trend_data_s["z"][k]=res.z

#winter trend
def SUMMER(column,division):
        Year = []
        summer = []
        for i in range(len(column)):
            if column[i] == division:
                if 'nan'!= str(data["DEC-FEB"][i]) and 'NAN' != str(data["DEC-FEB"][i]):
                    Year.append(data["YEAR"][i])
                    summer.append((data["DEC-FEB"][i]))
        return Year,summer
cols=["State","trend","h","p","z"]
trend_data_w=pd.DataFrame(columns=cols,index=(range(16)))
for k in range(len(country)):
    Year, summer = SUMMER(data["State"],str(country[k]))
    res = mk.hamed_rao_modification_test(summer,alpha=.05)
    trend_data_w["State"][k]=country[k]
    trend_data_w["trend"][k]=res.trend
    trend_data_w["h"][k]=res.h
    trend_data_w["p"][k]=res.p
    trend_data_w["z"][k]=res.z

#spring trend
def SUMMER(column,division):
        Year = []
        summer = []
        for i in range(len(column)):
            if column[i] == division:
                if 'nan'!= str(data["MAR-MAY"][i]) and 'NAN' != str(data["MAR-MAY"][i]):
                    Year.append(data["YEAR"][i])
                    summer.append((data["MAY"][i]))
        return Year,summer
cols=["State","trend","h","p","z"]
trend_data_sp=pd.DataFrame(columns=cols,index=(range(16)))
for k in range(len(country)):
    Year, summer = SUMMER(data["State"],str(country[k]))
    res = mk.hamed_rao_modification_test(summer,alpha=.05)
    trend_data_sp["State"][k]=country[k]
    trend_data_sp["trend"][k]=res.trend
    trend_data_sp["h"][k]=res.h
    trend_data_sp["p"][k]=res.p
    trend_data_sp["z"][k]=res.z
    
# table 
trends=pd.DataFrame({"Name":trend_data_s.State,"DEC-FEB":trend_data_w.z,"MAR-MAY":trend_data_sp.z,
                     "JUN-AUG":trend_data_s.z,"SEP-NOV":trend_data_f.z,"ANNUAL":trend_data.z})
#trends.to_excel("plots/Hamadrao_modification_test.xlsx")<jupyter_output><empty_output><jupyter_text>Change point anallysis by Pettitt test<jupyter_code>import pyhomogeneity as hg
cp_data = data.groupby(["State"])["ANNUAL","YEAR"]
cp_df=pd.DataFrame(columns=["states","Change Point","p","h"],index=(range(16)))
states=data["State"].unique()
for q in range(len(country)):
    state=country[q]
    grp=cp_data.get_group(state)
    grp=grp.set_index("YEAR")
    grp.index=pd.to_datetime(grp.index,format="%Y")
    result = hg.pettitt_test(grp,alpha=0.05)
    cp_df["states"][q]=state
    cp_df["Change Point"][q]=result.cp[0:4]
    cp_df["p"][q]=result.p
    cp_df["h"][q]=result.h
cp_df.to_csv("Pettit test annual.csv")<jupyter_output><empty_output><jupyter_text>Train test splitiing<jupyter_code>#tarin tes split for the year 2010
from sklearn.model_selection import train_test_split

data_2010 = np.asarray(data[['JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL',
       'AUG', 'SEP', 'OCT', 'NOV', 'DEC']].loc[(data['YEAR'] >= 2010) & (data["YEAR"] <=2015)])

#data_2010 = np.asarray(data[['JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL',
       #'AUG', 'SEP', 'OCT', 'NOV', 'DEC']].loc[(data['YEAR'] == 2010)])

#data_2010_brandenburg = np.asarray(data_2010[['JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL',
       #'AUG', 'SEP', 'OCT', 'NOV', 'DEC']].loc[data_2010['State'] == 'Brandenburg/Berlin'])

X_year_2010 = None; y_year_2010 = None
for i in range(data_2010.shape[1]-3):
    if X_year_2010 is None:
        X_year_2010 = data_2010[:, i:i+3]
        y_year_2010 = data_2010[:, i+3]
    else:
        X_year_2010 = np.concatenate((X_year_2010, data_2010[:, i:i+3]), axis=0)
        y_year_2010 = np.concatenate((y_year_2010, data_2010[:, i+3]), axis=0)
        
X_train, X_test, y_train, y_test = train_test_split(X_year_2010,y_year_2010,test_size = 0.20,random_state = 5)
from cycler import cycler
color_c = cycler('color', ['k'])
style_c = cycler('linestyle', ['-', '--', ':', '-.'])
markr_c = cycler('marker', ['', '.', 'o'])
c_cms = color_c * markr_c * style_c
c_csm = color_c * style_c * markr_c
plt.rc('font', **font)
plt.rc('axes', prop_cycle=c_cms)
###### knn model for 2010 data
%matplotlib
from sklearn.neighbors import KNeighborsRegressor
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_absolute_percentage_error
from sklearn.metrics import accuracy_score
neigh = KNeighborsRegressor(n_neighbors=40)
neigh.fit(X_train,y_train)

predicted_knn = neigh.predict(X_test)
#print("Prediction Result: ",predicted1)
mean_abs_per = mean_absolute_percentage_error(predicted_knn,y_test)
#print('\nR-squared test score: {:.3f}'.format(neigh.score(X_test_i,y_test_i))) # R-Squared test score
plt.rcParams['figure.figsize']=(16,4)
plt.plot(y_test,label='Actual')
plt.plot(predicted_knn,label='Predicted by k-NN model')
plt.xlabel('Test Size',alpha=0.8,fontsize=20)
plt.ylabel('Precipitation ( in mm)',alpha=0.8,fontsize=20)
plt.xticks(alpha=0.6,fontsize=20)
plt.yticks(alpha=0.6,fontsize=20)
plt.legend(fontsize=20,loc=1)
plt.text(100,175,"MAPE: "+str(mean_abs_per.round(3)),fontsize=20,horizontalalignment='center')
#plt.text(12,125,"MAPE: "+str(mean_abs_per.round(3)),fontsize=20,horizontalalignment='center')
plt.show()
plt.tight_layout()
#plt.savefig("knn_k20.jpg",dpi=150)

#from sklearn.utils import check_arrays
#mean_abs = mean_absolute_error(predicted_knn,y_test)

print(mean_abs_per)

#plt.savefig("plots/knn_40.jpg",dpi=150)
#SVM  model for 2010 data
from sklearn.svm import SVR
# SVM model
clf = SVR(kernel='rbf', gamma="scale", C=150, epsilon=0.1)
clf.fit(X_train, y_train) 
pred_svr = clf.predict(X_test)
mean_abs_per_svr = mean_absolute_percentage_error(pred_svr,y_test)
#print (mean_absolute_error(y_test, y_pred))
plt.rcParams['figure.figsize']=(16,4)
plt.plot(y_test,label='Actual')
plt.plot(pred_svr,label='Predicted by SVR')
plt.xlabel('Test Size',fontsize=20)
plt.ylabel('Precipitation ( in mm)',fontsize=20)
plt.xticks(alpha=0.6,fontsize=20)
plt.yticks(alpha=0.6,fontsize=20)
plt.text(100,175,"MAPE: "+str(mean_abs_per_svr.round(3)),fontsize=20,horizontalalignment='center')
#plt.text(12,125,"MAPE: "+str(mean_abs_per_svr.round(3)),fontsize=20,horizontalalignment='center')
plt.legend(fontsize=20,loc=1)
plt.tight_layout()
plt.show()
#plt.savefig("plots/SVR_C150.jpg",dpi=150)

print(mean_abs_per_svr)
from sklearn.ensemble import RandomForestRegressor
regressor = RandomForestRegressor(n_estimators=500, max_depth=None, max_features="auto", min_samples_leaf=1, min_samples_split=2, bootstrap=True)
regressor.fit(X_train, y_train)

predicted_rfr = regressor.predict(X_test)
mean_abs_per_rfr = mean_absolute_percentage_error(predicted_rfr,y_test)
plt.rcParams['figure.figsize']=(16,4)
plt.plot(y_test,label='Actual')
plt.plot(predicted_rfr,label='Predicted by RFR')
plt.xlabel('Test Size',size=20)
plt.ylabel('Precipitation ( in mm)',size=20)
plt.xticks(alpha=0.6,fontsize=20)
plt.yticks(alpha=0.6,fontsize=20)
plt.legend(fontsize=20,loc=1)
plt.text(100,175,"MAPE: "+str(mean_abs_per_rfr.round(3)),fontsize=20,horizontalalignment='center')
#plt.text(12,125,"MAPE: "+str(mean_abs_per_rfr.round(3)),fontsize=20,horizontalalignment='center')
plt.tight_layout()
plt.show()
plt.savefig("RFR.jpg",dpi=150)
mean_abs_per_rfr = mean_absolute_percentage_error(predicted_rfr,y_test)
print(mean_abs_per_rfr)
plt.savefig("plots/RFR_N500.jpg",dpi=150)
from cycler import cycler
color_c = cycler('color', ['k'])
style_c = cycler('linestyle', ['-', '--', ':', '-.'])
markr_c = cycler('marker', ['', '.', 'o'])
c_cms = color_c * markr_c * style_c
c_csm = color_c * style_c * markr_c

plt.rc('axes', prop_cycle=c_cms)
#fig.patch.set_facecolor('white')
fig,(ax1,ax2,ax3,ax4) = plt.subplots(4,1)
ax1.plot(y_test,label="Actual")
ax2.plot(predicted_knn,label="k-NN")
ax3.plot(pred_svr,label="SVR")
plt.plot(predicted_rfr,label="RFR")
plt.legend()
plt.show()
from sklearn import linear_model

# linear model
reg = linear_model.ElasticNet(alpha=0.1)
reg.fit(X_train, y_train)
pred_lin = reg.predict(X_test)

mean_abs_per_lin = mean_absolute_percentage_error(pred_lin,y_test)
plt.rcParams['figure.figsize']=(16,4)
plt.plot(y_test,label='Actual')
plt.plot(pred_lin,label='Predicted by Linear model')
plt.xlabel('Test Size',size=20)
plt.ylabel('Precipitation ( in mm)',size=20)
plt.xticks(alpha=0.6,fontsize=20)
plt.yticks(alpha=0.6,fontsize=20)
plt.legend(fontsize=20)
plt.text(13,138,"MAPE: "+str(mean_abs_per_lin.round(3)),fontsize=20,horizontalalignment='center')
plt.tight_layout()
plt.show()
plt.savefig("plots/linear.jpg",dpi=150)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2038670921719542, 'fraction_numerical': 0.04298854890182091, 'mean_word_length': 5.3219936708860756, 'pattern_counts': {'":': 24, '<': 32, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '557480', 'license': 'MIT', 'n_tokens_mistral': 6966, 'n_tokens_neox': 6417, 'n_words': 827, 'path': 'trend_analysis_germany.ipynb', 'repo_name': 'talukmdm/Trend-analysis-and-prediction-of-precipitation', 'size': '16000'}