#HARROW
#dat1=read.csv("mm6130257.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
#dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(42.03,length(dat1$Year)),rep(-82.9,length(dat1$Year)),rep(182,length(dat1$Year))))
#RIDGETOWN
dat1=read.csv("mm6137154.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(42.45,length(dat1$Year)),rep(-81.88,length(dat1$Year)),rep(206,length(dat1$Year))))
#VINELAND
dat1=read.csv("mm6139148.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.17,length(dat1$Year)),rep(-79.42,length(dat1$Year)),rep(79,length(dat1$Year))))
#WELLAND
dat1=read.csv("mm6139449.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43,length(dat1$Year)),rep(-79.27,length(dat1$Year)),rep(175,length(dat1$Year))))
#WINDSOR
dat1=read.csv("mm6139527.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(42.27,length(dat1$Year)),rep(-82.97,length(dat1$Year)),rep(190,length(dat1$Year))))
#LONDON
dat1=read.csv("mm6144478.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.03,length(dat1$Year)),rep(-81.15,length(dat1$Year)),rep(278,length(dat1$Year))))
#WOODSTOCK
dat1=read.csv("mm6149625.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.13,length(dat1$Year)),rep(-80.77,length(dat1$Year)),rep(282,length(dat1$Year))))
#BELLEVILLE
dat1=read.csv("mm6150689.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.15,length(dat1$Year)),rep(-77.4,length(dat1$Year)),rep(76,length(dat1$Year))))
#HAMILTON
dat1=read.csv("mm6153193.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.17,length(dat1$Year)),rep(-79.93,length(dat1$Year)),rep(238,length(dat1$Year))))
#ORANGEVILLE
dat1=read.csv("mm6155790.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.92,length(dat1$Year)),rep(-80.08,length(dat1$Year)),rep(412,length(dat1$Year))))
#TORONTO
dat1=read.csv("mm6158355.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.67,length(dat1$Year)),rep(-79.4,length(dat1$Year)),rep(113,length(dat1$Year))))
#HALIBURTON
dat1=read.csv("mm6163171.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.03,length(dat1$Year)),rep(-78.53,length(dat1$Year)),rep(330,length(dat1$Year))))
#PETERBOROUGH
dat1=read.csv("mm6166415.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.23,length(dat1$Year)),rep(-78.37,length(dat1$Year)),rep(191,length(dat1$Year))))
#ATIKOKAN
#dat1=read.csv("mm6020LPQ.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
#dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.8,length(dat1$Year)),rep(-91.58,length(dat1$Year)),rep(442,length(dat1$Year))))
#select the data in 1978-2004
dat3=dat2[dat2[,1]>=1978 & dat2[,1]<=2017,]
#select the data in 1978-2017
dat3=dat2[dat2[,1]>=1978 & dat2[,1]<=2017,]
#marginal model(Sequential Estimation)
#Jan
y=dat3[dat3[,2]!=-9999.9,2]
x1=dat3[dat3[,2]!=-9999.9,1]
x2=dat3[dat3[,2]!=-9999.9,14]
x3=dat3[dat3[,2]!=-9999.9,15]
x4=dat3[dat3[,2]!=-9999.9,16]
model1=lm(y~x1+x2+x3+x4)
#Feb
y=dat3[dat3[,3]!=-9999.9,3]
x1=dat3[dat3[,3]!=-9999.9,1]
x2=dat3[dat3[,3]!=-9999.9,14]
x3=dat3[dat3[,3]!=-9999.9,15]
x4=dat3[dat3[,3]!=-9999.9,16]
model2=lm(y~x2+x3+x4)
#March
y=dat3[dat3[,4]!=-9999.9,4]
x1=dat3[dat3[,4]!=-9999.9,1]
x2=dat3[dat3[,4]!=-9999.9,14]
x3=dat3[dat3[,4]!=-9999.9,15]
x4=dat3[dat3[,4]!=-9999.9,16]
model3=lm(y~x1+x2+x3+x4)
#April
y=dat3[dat3[,5]!=-9999.9,5]
x1=dat3[dat3[,5]!=-9999.9,1]
x2=dat3[dat3[,5]!=-9999.9,14]
x3=dat3[dat3[,5]!=-9999.9,15]
x4=dat3[dat3[,5]!=-9999.9,16]
model4=lm(y~x2+x3+x4)
#May(linear model not good)
quadlon=poly(dat3[,15],2)
y=dat3[dat3[,6]!=-9999.9,6]
x1=dat3[dat3[,6]!=-9999.9,1]
x2=dat3[dat3[,6]!=-9999.9,14]
x3=dat3[dat3[,6]!=-9999.9,15]
x4=dat3[dat3[,6]!=-9999.9,16]
model5=lm(y~x2+poly(x3,2)+x4)
model5=lm(y~x1+x2+x3+x4)
plot(dat3[dat3[,6]!=-9999.9,1],dat3[dat3[,6]!=-9999.9,6])
plot(dat3[dat3[,6]!=-9999.9,14],dat3[dat3[,6]!=-9999.9,6])
plot(dat3[dat3[,6]!=-9999.9,15],dat3[dat3[,6]!=-9999.9,6])
plot(dat3[dat3[,6]!=-9999.9,16],dat3[dat3[,6]!=-9999.9,6])
#June
y=dat3[dat3[,7]!=-9999.9,7]
x1=dat3[dat3[,7]!=-9999.9,1]
x2=dat3[dat3[,7]!=-9999.9,14]
x3=dat3[dat3[,7]!=-9999.9,15]
x4=dat3[dat3[,7]!=-9999.9,16]
model6=lm(y~x1+x2+poly(x3,2)+x4)
model6=lm(y~x1+x2+x3+x4)
#July
y=dat3[dat3[,8]!=-9999.9,8]
x1=dat3[dat3[,8]!=-9999.9,1]
x2=dat3[dat3[,8]!=-9999.9,14]
x3=dat3[dat3[,8]!=-9999.9,15]
x4=dat3[dat3[,8]!=-9999.9,16]
model7=lm(y~x2+poly(x3,2)+x4)
model7=lm(y~x1+x2+x3+x4)
#Aug
y=dat3[dat3[,9]!=-9999.9,9]
x1=dat3[dat3[,9]!=-9999.9,1]
x2=dat3[dat3[,9]!=-9999.9,14]
x3=dat3[dat3[,9]!=-9999.9,15]
x4=dat3[dat3[,9]!=-9999.9,16]
model8=lm(y~x1+x2+poly(x3,2)+x4)
model8=lm(y~x1+x2+x3+x4)
#Sept
y=dat3[dat3[,10]!=-9999.9,10]
x1=dat3[dat3[,10]!=-9999.9,1]
x2=dat3[dat3[,10]!=-9999.9,14]
x3=dat3[dat3[,10]!=-9999.9,15]
x4=dat3[dat3[,10]!=-9999.9,16]
model9=lm(y~x1+x2+x3+x4)
#Oct
y=dat3[dat3[,11]!=-9999.9,11]
x1=dat3[dat3[,11]!=-9999.9,1]
x2=dat3[dat3[,11]!=-9999.9,14]
x3=dat3[dat3[,11]!=-9999.9,15]
x4=dat3[dat3[,11]!=-9999.9,16]
model10=lm(y~x1+x2+x3+x4)
#Nov
y=dat3[dat3[,12]!=-9999.9,12]
x1=dat3[dat3[,12]!=-9999.9,1]
x2=dat3[dat3[,12]!=-9999.9,14]
x3=dat3[dat3[,12]!=-9999.9,15]
x4=dat3[dat3[,12]!=-9999.9,16]
model11=lm(y~x1+x2+x3+x4)
#Dec
y=dat3[dat3[,13]!=-9999.9,13]
x1=dat3[dat3[,13]!=-9999.9,1]
x2=dat3[dat3[,13]!=-9999.9,14]
x3=dat3[dat3[,13]!=-9999.9,15]
x4=dat3[dat3[,13]!=-9999.9,16]
model12=lm(y~x1+x2+x3+x4)
#missing data imputation
j=2
dat3[dat3[,j]==-9999.9,j]=model1$coefficients[1]+model1$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model1$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model1$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model1$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=3
dat3[dat3[,j]==-9999.9,j]=model2$coefficients[1]+model2$coefficients[2]*dat3[dat3[,j]==-9999.9,14]+model2$coefficients[3]*dat3[dat3[,j]==-9999.9,15]+model2$coefficients[4]*dat3[dat3[,j]==-9999.9,16]
j=4
dat3[dat3[,j]==-9999.9,j]=model3$coefficients[1]+model3$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model3$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model3$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model3$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=5
dat3[dat3[,j]==-9999.9,j]=model4$coefficients[1]+model4$coefficients[2]*dat3[dat3[,j]==-9999.9,14]+model4$coefficients[3]*dat3[dat3[,j]==-9999.9,15]+model4$coefficients[4]*dat3[dat3[,j]==-9999.9,16]
j=6
dat3[dat3[,j]==-9999.9,j]=model5$coefficients[1]+model5$coefficients[2]*dat3[dat3[,j]==-9999.9,14]+model5$coefficients[3]*quadlon[dat3[,j]==-9999.9,1]+model5$coefficients[4]*quadlon[dat3[,j]==-9999.9,2]+model5$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=7
newdata=poly(dat3[dat3[,j]==-9999.9,15],2)
dat3[dat3[,j]==-9999.9,j]=model6$coefficients[1]+model6$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model6$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model6$coefficients[4]*quadlon[dat3[,j]==-9999.9,1]+model6$coefficients[5]*quadlon[dat3[,j]==-9999.9,2]+model6$coefficients[6]*dat3[dat3[,j]==-9999.9,16]
j=8
newdata=poly(dat3[dat3[,j]==-9999.9,15],2)
dat3[dat3[,j]==-9999.9,j]=model7$coefficients[1]+model7$coefficients[2]*dat3[dat3[,j]==-9999.9,14]+model7$coefficients[3]*quadlon[dat3[,j]==-9999.9,1]+model7$coefficients[4]*quadlon[dat3[,j]==-9999.9,2]+model7$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=9
newdata=poly(dat3[dat3[,j]==-9999.9,15],2)
dat3[dat3[,j]==-9999.9,j]=model8$coefficients[1]+model8$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model8$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model8$coefficients[4]*quadlon[dat3[,j]==-9999.9,1]+model8$coefficients[5]*quadlon[dat3[,j]==-9999.9,2]+model8$coefficients[6]*dat3[dat3[,j]==-9999.9,16]
j=10
dat3[dat3[,j]==-9999.9,j]=model9$coefficients[1]+model9$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model9$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model9$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model9$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=11
dat3[dat3[,j]==-9999.9,j]=model10$coefficients[1]+model10$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model10$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model10$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model10$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=12
dat3[dat3[,j]==-9999.9,j]=model11$coefficients[1]+model11$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model11$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model11$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model11$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=13
dat3[dat3[,j]==-9999.9,j]=model12$coefficients[1]+model12$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model12$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model12$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model12$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
dat3
length(dat3[is.na(dat3)==TRUE])
19092-79
j=7
newdata=poly(dat3[is.na(dat3[,j])==TRUE,15],2)
dat3[is.na(dat3[,j])==TRUE,j]=model6$coefficients[1]+model6$coefficients[2]*dat3[is.na(dat3[,j])==TRUE,1]+model6$coefficients[3]*dat3[is.na(dat3[,j])==TRUE,14]+model6$coefficients[4]*quadlon[is.na(dat3[,j])==TRUE,1]+model6$coefficients[5]*quadlon[is.na(dat3[,j])==TRUE,2]+model6$coefficients[6]*dat3[is.na(dat3[,j])==TRUE,16]
length(dat3[is.na(dat3)==TRUE])
dat3[is.na(dat3[,j])==TRUE,j]
model6$coefficients[1]+model6$coefficients[2]*dat3[is.na(dat3[,j])==TRUE,1]+model6$coefficients[3]*dat3[is.na(dat3[,j])==TRUE,14]+model6$coefficients[4]*quadlon[is.na(dat3[,j])==TRUE,1]+model6$coefficients[5]*quadlon[is.na(dat3[,j])==TRUE,2]+model6$coefficients[6]*dat3[is.na(dat3[,j])==TRUE,16]
model6
dat2=NULL
#Gore Bay
dat1=read.csv("mm6092920.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.88,length(dat1$Year)),rep(-82.57,length(dat1$Year)),rep(194,length(dat1$Year)))
colnames(dat2)=c("Year","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sept","Oct","Nov","Dec","latitude","longitude","elevation")
#BIG TROUT LAKE
#dat1=read.csv("mm6010735.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
#dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(53.83,length(dat1$Year)),rep(-89.87,length(dat1$Year)),rep(224,length(dat1$Year))))
#LANSDOWNE HOUSE
dat1=read.csv("mm6014353.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(52.23,length(dat1$Year)),rep(-87.88,length(dat1$Year)),rep(255,length(dat1$Year))))
#PICKLE LAKE
dat1=read.csv("mm6016529.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(51.45,length(dat1$Year)),rep(-90.22,length(dat1$Year)),rep(386,length(dat1$Year))))
#RED LAKE
dat1=read.csv("mm6016970.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(51.07,length(dat1$Year)),rep(-93.8,length(dat1$Year)),rep(386,length(dat1$Year))))
#FORT FRANCES
dat1=read.csv("mm6022474.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.65,length(dat1$Year)),rep(-93.43,length(dat1$Year)),rep(342,length(dat1$Year))))
#MINE CENTRE
dat1=read.csv("mm6025205.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.8,length(dat1$Year)),rep(-92.6,length(dat1$Year)),rep(361,length(dat1$Year))))
#DRYDEN
dat1=read.csv("mm6032125.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.78,length(dat1$Year)),rep(-92.83,length(dat1$Year)),rep(413,length(dat1$Year))))
#KENORA
dat1=read.csv("mm6034076.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.78,length(dat1$Year)),rep(-94.37,length(dat1$Year)),rep(227,length(dat1$Year))))
#SIOUX LOOKOUT
#dat1=read.csv("mm6037776.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
#dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(50.12,length(dat1$Year)),rep(-91.9,length(dat1$Year)),rep(383,length(dat1$Year))))
#CAMERON FALLS
dat1=read.csv("mm6041110.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.15,length(dat1$Year)),rep(-88.35,length(dat1$Year)),rep(233,length(dat1$Year))))
#GERALDTON
dat1=read.csv("mm6042716.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.78,length(dat1$Year)),rep(-86.93,length(dat1$Year)),rep(349,length(dat1$Year))))
#THUNDER BAY
dat1=read.csv("mm6048260.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.37,length(dat1$Year)),rep(-89.33,length(dat1$Year)),rep(199,length(dat1$Year))))
#HORNEPAYNE
dat1=read.csv("mm6053575.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.2,length(dat1$Year)),rep(-84.77,length(dat1$Year)),rep(335,length(dat1$Year))))
#SAULT STE MARIE
dat1=read.csv("mm6057591.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(46.48,length(dat1$Year)),rep(-84.52,length(dat1$Year)),rep(192,length(dat1$Year))))
#WAWA
dat1=read.csv("mm6059407.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(47.97,length(dat1$Year)),rep(-84.78,length(dat1$Year)),rep(287,length(dat1$Year))))
#CHAPLEAU
dat1=read.csv("mm6061361.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(47.82,length(dat1$Year)),rep(-83.35,length(dat1$Year)),rep(447,length(dat1$Year))))
#SUDBURY
dat1=read.csv("mm6068153.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(46.62,length(dat1$Year)),rep(-80.8,length(dat1$Year)),rep(348,length(dat1$Year))))
#EARLTON
dat1=read.csv("mm6072230.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(47.7,length(dat1$Year)),rep(-79.85,length(dat1$Year)),rep(243,length(dat1$Year))))
#IROQUOIS FALLS
dat1=read.csv("mm6073810.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.75,length(dat1$Year)),rep(-80.67,length(dat1$Year)),rep(259,length(dat1$Year))))
#KAPUSKASING
dat1=read.csv("mm6073976.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(49.42,length(dat1$Year)),rep(-82.47,length(dat1$Year)),rep(227,length(dat1$Year))))
#MOOSONEE
dat1=read.csv("mm6075425.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(51.27,length(dat1$Year)),rep(-80.65,length(dat1$Year)),rep(10,length(dat1$Year))))
#SMOKY FALLS
dat1=read.csv("mm6077845.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(50.07,length(dat1$Year)),rep(-82.17,length(dat1$Year)),rep(183,length(dat1$Year))))
#TIMMINS
dat1=read.csv("mm6078286.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.57,length(dat1$Year)),rep(-81.38,length(dat1$Year)),rep(295,length(dat1$Year))))
#MADAWASKA
dat1=read.csv("mm6080192.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.5,length(dat1$Year)),rep(-77.98,length(dat1$Year)),rep(316,length(dat1$Year))))
#NORTH BAY
dat1=read.csv("mm6085701.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(46.37,length(dat1$Year)),rep(-79.42,length(dat1$Year)),rep(370,length(dat1$Year))))
#BROCKVILLE
dat1=read.csv("mm6100971.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.6,length(dat1$Year)),rep(-75.67,length(dat1$Year)),rep(96,length(dat1$Year))))
#CORNWALL
dat1=read.csv("mm6101874.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.02,length(dat1$Year)),rep(-74.75,length(dat1$Year)),rep(64,length(dat1$Year))))
#KINGSTON
dat1=read.csv("mm6104142.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.22,length(dat1$Year)),rep(-76.6,length(dat1$Year)),rep(93,length(dat1$Year))))
#MORRISBURG
dat1=read.csv("mm6105460.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.92,length(dat1$Year)),rep(-75.18,length(dat1$Year)),rep(82,length(dat1$Year))))
#OTTAWA
dat1=read.csv("mm6105976.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.38,length(dat1$Year)),rep(-75.72,length(dat1$Year)),rep(79,length(dat1$Year))))
#BEATRICE
#dat1=read.csv("mm6110607.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
#dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.13,length(dat1$Year)),rep(-79.4,length(dat1$Year)),rep(297,length(dat1$Year))))
#OWEN SOUND
dat1=read.csv("mm6116132.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.58,length(dat1$Year)),rep(-80.93,length(dat1$Year)),rep(179,length(dat1$Year))))
#HARROW
#dat1=read.csv("mm6130257.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
#dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(42.03,length(dat1$Year)),rep(-82.9,length(dat1$Year)),rep(182,length(dat1$Year))))
#RIDGETOWN
dat1=read.csv("mm6137154.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(42.45,length(dat1$Year)),rep(-81.88,length(dat1$Year)),rep(206,length(dat1$Year))))
#VINELAND
dat1=read.csv("mm6139148.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.17,length(dat1$Year)),rep(-79.42,length(dat1$Year)),rep(79,length(dat1$Year))))
#WELLAND
dat1=read.csv("mm6139449.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43,length(dat1$Year)),rep(-79.27,length(dat1$Year)),rep(175,length(dat1$Year))))
#WINDSOR
dat1=read.csv("mm6139527.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(42.27,length(dat1$Year)),rep(-82.97,length(dat1$Year)),rep(190,length(dat1$Year))))
#LONDON
dat1=read.csv("mm6144478.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.03,length(dat1$Year)),rep(-81.15,length(dat1$Year)),rep(278,length(dat1$Year))))
#WOODSTOCK
dat1=read.csv("mm6149625.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.13,length(dat1$Year)),rep(-80.77,length(dat1$Year)),rep(282,length(dat1$Year))))
#BELLEVILLE
dat1=read.csv("mm6150689.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.15,length(dat1$Year)),rep(-77.4,length(dat1$Year)),rep(76,length(dat1$Year))))
#HAMILTON
dat1=read.csv("mm6153193.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.17,length(dat1$Year)),rep(-79.93,length(dat1$Year)),rep(238,length(dat1$Year))))
#ORANGEVILLE
dat1=read.csv("mm6155790.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.92,length(dat1$Year)),rep(-80.08,length(dat1$Year)),rep(412,length(dat1$Year))))
#TORONTO
dat1=read.csv("mm6158355.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(43.67,length(dat1$Year)),rep(-79.4,length(dat1$Year)),rep(113,length(dat1$Year))))
#HALIBURTON
dat1=read.csv("mm6163171.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(45.03,length(dat1$Year)),rep(-78.53,length(dat1$Year)),rep(330,length(dat1$Year))))
#PETERBOROUGH
dat1=read.csv("mm6166415.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.23,length(dat1$Year)),rep(-78.37,length(dat1$Year)),rep(191,length(dat1$Year))))
#ATIKOKAN
#dat1=read.csv("mm6020LPQ.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
#dat2=rbind(dat2,cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(48.8,length(dat1$Year)),rep(-91.58,length(dat1$Year)),rep(442,length(dat1$Year))))
#select the data in 1978-2004
dat3=dat2[dat2[,1]>=1978 & dat2[,1]<=2017,]
#select the data in 1978-2017
dat3=dat2[dat2[,1]>=1978 & dat2[,1]<=2017,]
#marginal model(Sequential Estimation)
#Jan
y=dat3[dat3[,2]!=-9999.9,2]
x1=dat3[dat3[,2]!=-9999.9,1]
x2=dat3[dat3[,2]!=-9999.9,14]
x3=dat3[dat3[,2]!=-9999.9,15]
x4=dat3[dat3[,2]!=-9999.9,16]
model1=lm(y~x1+x2+x3+x4)
#Feb
y=dat3[dat3[,3]!=-9999.9,3]
x1=dat3[dat3[,3]!=-9999.9,1]
x2=dat3[dat3[,3]!=-9999.9,14]
x3=dat3[dat3[,3]!=-9999.9,15]
x4=dat3[dat3[,3]!=-9999.9,16]
model2=lm(y~x2+x3+x4)
#March
y=dat3[dat3[,4]!=-9999.9,4]
x1=dat3[dat3[,4]!=-9999.9,1]
x2=dat3[dat3[,4]!=-9999.9,14]
x3=dat3[dat3[,4]!=-9999.9,15]
x4=dat3[dat3[,4]!=-9999.9,16]
model3=lm(y~x1+x2+x3+x4)
#April
y=dat3[dat3[,5]!=-9999.9,5]
x1=dat3[dat3[,5]!=-9999.9,1]
x2=dat3[dat3[,5]!=-9999.9,14]
x3=dat3[dat3[,5]!=-9999.9,15]
x4=dat3[dat3[,5]!=-9999.9,16]
model4=lm(y~x2+x3+x4)
#May(linear model not good)
quadlon=poly(dat3[,15],2)
y=dat3[dat3[,6]!=-9999.9,6]
x1=dat3[dat3[,6]!=-9999.9,1]
x2=dat3[dat3[,6]!=-9999.9,14]
x3=dat3[dat3[,6]!=-9999.9,15]
x4=dat3[dat3[,6]!=-9999.9,16]
model5=lm(y~x2+poly(x3,2)+x4)
#model5=lm(y~x1+x2+x3+x4)
plot(dat3[dat3[,6]!=-9999.9,1],dat3[dat3[,6]!=-9999.9,6])
plot(dat3[dat3[,6]!=-9999.9,14],dat3[dat3[,6]!=-9999.9,6])
plot(dat3[dat3[,6]!=-9999.9,15],dat3[dat3[,6]!=-9999.9,6])
plot(dat3[dat3[,6]!=-9999.9,16],dat3[dat3[,6]!=-9999.9,6])
#June
y=dat3[dat3[,7]!=-9999.9,7]
x1=dat3[dat3[,7]!=-9999.9,1]
x2=dat3[dat3[,7]!=-9999.9,14]
x3=dat3[dat3[,7]!=-9999.9,15]
x4=dat3[dat3[,7]!=-9999.9,16]
model6=lm(y~x1+x2+poly(x3,2)+x4)
#model6=lm(y~x1+x2+x3+x4)
#July
y=dat3[dat3[,8]!=-9999.9,8]
x1=dat3[dat3[,8]!=-9999.9,1]
x2=dat3[dat3[,8]!=-9999.9,14]
x3=dat3[dat3[,8]!=-9999.9,15]
x4=dat3[dat3[,8]!=-9999.9,16]
model7=lm(y~x2+poly(x3,2)+x4)
#model7=lm(y~x1+x2+x3+x4)
#Aug
y=dat3[dat3[,9]!=-9999.9,9]
x1=dat3[dat3[,9]!=-9999.9,1]
x2=dat3[dat3[,9]!=-9999.9,14]
x3=dat3[dat3[,9]!=-9999.9,15]
x4=dat3[dat3[,9]!=-9999.9,16]
model8=lm(y~x1+x2+poly(x3,2)+x4)
#model8=lm(y~x1+x2+x3+x4)
#Sept
y=dat3[dat3[,10]!=-9999.9,10]
x1=dat3[dat3[,10]!=-9999.9,1]
x2=dat3[dat3[,10]!=-9999.9,14]
x3=dat3[dat3[,10]!=-9999.9,15]
x4=dat3[dat3[,10]!=-9999.9,16]
model9=lm(y~x1+x2+x3+x4)
#Oct
y=dat3[dat3[,11]!=-9999.9,11]
x1=dat3[dat3[,11]!=-9999.9,1]
x2=dat3[dat3[,11]!=-9999.9,14]
x3=dat3[dat3[,11]!=-9999.9,15]
x4=dat3[dat3[,11]!=-9999.9,16]
model10=lm(y~x1+x2+x3+x4)
#Nov
y=dat3[dat3[,12]!=-9999.9,12]
x1=dat3[dat3[,12]!=-9999.9,1]
x2=dat3[dat3[,12]!=-9999.9,14]
x3=dat3[dat3[,12]!=-9999.9,15]
x4=dat3[dat3[,12]!=-9999.9,16]
model11=lm(y~x1+x2+x3+x4)
#Dec
y=dat3[dat3[,13]!=-9999.9,13]
x1=dat3[dat3[,13]!=-9999.9,1]
x2=dat3[dat3[,13]!=-9999.9,14]
x3=dat3[dat3[,13]!=-9999.9,15]
x4=dat3[dat3[,13]!=-9999.9,16]
model12=lm(y~x1+x2+x3+x4)
#missing data imputation
j=2
dat3[dat3[,j]==-9999.9,j]=model1$coefficients[1]+model1$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model1$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model1$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model1$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=3
dat3[dat3[,j]==-9999.9,j]=model2$coefficients[1]+model2$coefficients[2]*dat3[dat3[,j]==-9999.9,14]+model2$coefficients[3]*dat3[dat3[,j]==-9999.9,15]+model2$coefficients[4]*dat3[dat3[,j]==-9999.9,16]
j=4
dat3[dat3[,j]==-9999.9,j]=model3$coefficients[1]+model3$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model3$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model3$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model3$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=5
dat3[dat3[,j]==-9999.9,j]=model4$coefficients[1]+model4$coefficients[2]*dat3[dat3[,j]==-9999.9,14]+model4$coefficients[3]*dat3[dat3[,j]==-9999.9,15]+model4$coefficients[4]*dat3[dat3[,j]==-9999.9,16]
j=6
dat3[dat3[,j]==-9999.9,j]=model5$coefficients[1]+model5$coefficients[2]*dat3[dat3[,j]==-9999.9,14]+model5$coefficients[3]*quadlon[dat3[,j]==-9999.9,1]+model5$coefficients[4]*quadlon[dat3[,j]==-9999.9,2]+model5$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=7
newdata=poly(dat3[dat3[,j]==-9999.9,15],2)
dat3[dat3[,j]==-9999.9,j]=model6$coefficients[1]+model6$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model6$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model6$coefficients[4]*quadlon[dat3[,j]==-9999.9,1]+model6$coefficients[5]*quadlon[dat3[,j]==-9999.9,2]+model6$coefficients[6]*dat3[dat3[,j]==-9999.9,16]
j=8
newdata=poly(dat3[dat3[,j]==-9999.9,15],2)
dat3[dat3[,j]==-9999.9,j]=model7$coefficients[1]+model7$coefficients[2]*dat3[dat3[,j]==-9999.9,14]+model7$coefficients[3]*quadlon[dat3[,j]==-9999.9,1]+model7$coefficients[4]*quadlon[dat3[,j]==-9999.9,2]+model7$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=9
newdata=poly(dat3[dat3[,j]==-9999.9,15],2)
dat3[dat3[,j]==-9999.9,j]=model8$coefficients[1]+model8$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model8$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model8$coefficients[4]*quadlon[dat3[,j]==-9999.9,1]+model8$coefficients[5]*quadlon[dat3[,j]==-9999.9,2]+model8$coefficients[6]*dat3[dat3[,j]==-9999.9,16]
j=10
dat3[dat3[,j]==-9999.9,j]=model9$coefficients[1]+model9$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model9$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model9$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model9$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=11
dat3[dat3[,j]==-9999.9,j]=model10$coefficients[1]+model10$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model10$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model10$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model10$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=12
dat3[dat3[,j]==-9999.9,j]=model11$coefficients[1]+model11$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model11$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model11$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model11$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
j=13
dat3[dat3[,j]==-9999.9,j]=model12$coefficients[1]+model12$coefficients[2]*dat3[dat3[,j]==-9999.9,1]+model12$coefficients[3]*dat3[dat3[,j]==-9999.9,14]+model12$coefficients[4]*dat3[dat3[,j]==-9999.9,15]+model12$coefficients[5]*dat3[dat3[,j]==-9999.9,16]
length(dat3[is.na(dat3)==TRUE])
respy=c(dat3[,2],dat3[,3],dat3[,4],dat3[,5],dat3[,6],dat3[,7],dat3[,8],dat3[,9],dat3[,10],dat3[,11],dat3[,12],dat3[,13])
latit=rep(dat3[,14],12)
longi=rep(dat3[,15],12)
eleva=rep(dat3[,16],12)
years=rep(dat3[,1],12)
tindex=c((dat3[,1]-1978)*12+1,(dat3[,1]-1978)*12+2,(dat3[,1]-1978)*12+3,(dat3[,1]-1978)*12+4,(dat3[,1]-1978)*12+5,(dat3[,1]-1978)*12+6,(dat3[,1]-1978)*12+7,(dat3[,1]-1978)*12+8,(dat3[,1]-1978)*12+9,(dat3[,1]-1978)*12+10,(dat3[,1]-1978)*12+11,(dat3[,1]-1978)*12+12)
mon=c(rep(1,length(dat3[,2])),rep(2,length(dat3[,2])),rep(3,length(dat3[,2])),rep(4,length(dat3[,2])),rep(5,length(dat3[,2])),rep(6,length(dat3[,2])),rep(7,length(dat3[,2])),rep(8,length(dat3[,2])),rep(9,length(dat3[,2])),rep(10,length(dat3[,2])),rep(11,length(dat3[,2])),rep(12,length(dat3[,2])))
fital=lm(respy~latit+longi+eleva)
res=respy-fital$coefficients[1]-fital$coefficients[2]*latit-fital$coefficients[3]*longi-fital$coefficients[4]*eleva
respy1=matrix(res,length(dat3[,2]),length(res)/length(dat3[,2]))
mod111 <- astsa::sarima(data.frame(respy1),2,1,0,1,1,1,12)
data.frame(respy1)
mod111 <- astsa::sarima(res,2,1,0,1,1,1,12)
summary(mod111)
mod111
mod111 <- astsa::sarima(res,2,1,0,1,0,1,12)
mod111
mod111 <- astsa::sarima(res,2,1,0,1,0,0,12)
mod111
mod111 <- astsa::sarima(res,2,1,0,0,0,1,12)
mod111
mod111 <- astsa::sarima(res,1,1,0,0,0,1,12)
mod111
mod111 <- astsa::sarima(res,3,1,0,0,0,1,12)
mod111
mod111 <- astsa::sarima(res,3,1,1,0,0,1,12)
mod111
mod111 <- astsa::sarima(res,3,1,2,0,0,1,12)
mod111
mod111 <- astsa::sarima(res,3,2,2,0,0,1,12)
mod111
mod111 <- astsa::sarima(res,4,1,2,0,0,1,12)
mod111
mod111 <- astsa::sarima(res,4,1,2,0,0,1,13)
mod111
mod111 <- astsa::sarima(res,4,1,2,0,0,1,11)
mod111
mod111 <- astsa::sarima(res,4,1,3,0,0,1,12)
mod111
mod111 <- astsa::sarima(res,4,1,3,0,0,1,36)
mod111
mod111 <- astsa::sarima(res,4,1,4,0,0,1,36)
mod111
#KINGSTON
dat1=read.csv("mm6104142.txt",header=FALSE,skip=4,,col.names=c("Year","Jan","Jan Status","Feb","Feb Status","Mar","Mar Status","Apr","Apr Status","May","May Status","Jun","Jun Status","Jul","Jul Status","Aug","Aug Status","Sept","Sept Status","Oct","Oct Status","Nov","Nov Status","Dec","Dec Status","Annual","Annual Status","Spring","Spring Status","Summer","Summer Status","Fall","Fall Status","Winter","Winter Status"))
dat2=cbind(dat1$Year,dat1$Jan,dat1$Feb,dat1$Mar,dat1$Apr,dat1$May,dat1$Jun,dat1$Jul,dat1$Aug,dat1$Sept,dat1$Oct,dat1$Nov,dat1$Dec,rep(44.22,length(dat1$Year)),rep(-76.6,length(dat1$Year)),rep(93,length(dat1$Year)))
dat5=dat2[dat2[,1]>=2009,]
dat5=dat5[-c(9),]
newdata=data.frame(x1=c(2009,2010,2011,2012,2013,2014,2015,2016,2018),x2=dat5[,14],x3=dat5[,15],x4=dat5[,16])
newdata1=data.frame(latit=rep(dat5[,14],12),longi=rep(dat5[,15],12),eleva=rep(dat5[,16],12),years=rep(newdata$x1,12),mon=c(rep(1,length(dat5[,1])),rep(2,length(dat5[,1])),rep(3,length(dat5[,1])),rep(4,length(dat5[,1])),rep(5,length(dat5[,1])),rep(6,length(dat5[,1])),rep(7,length(dat5[,1])),rep(8,length(dat5[,1])),rep(9,length(dat5[,1])),rep(10,length(dat5[,1])),rep(11,length(dat5[,1])),rep(12,length(dat5[,1]))),tindex=c((newdata$x1-1978)*12+1,(newdata$x1-1978)*12+2,(newdata$x1-1978)*12+3,(newdata$x1-1978)*12+4,(newdata$x1-1978)*12+5,(newdata$x1-1978)*12+6,(newdata$x1-1978)*12+7,(newdata$x1-1978)*12+8,(newdata$x1-1978)*12+9,(newdata$x1-1978)*12+10,(newdata$x1-1978)*12+11,(newdata$x1-1978)*12+12))
predict(mod111,newdata1)
mod111
sarima.for(res,10,4,1,3,0,0,1,12)
plot(res)
respy1
mod111 <- astsa::sarima(as.vector(t(respy1)),4,1,3,0,0,1,12)
mod111
plot(as.vector(t(respy1)))
astsa::sarima(as.vector(t(respy1)),3,1,3,0,0,1,12)
astsa::sarima(as.vector(t(respy1)),3,1,3,0,0,1,11)
astsa::sarima(as.vector(t(respy1)),3,1,3,0,0,1,13)
astsa::sarima(as.vector(t(respy1)),4,1,3,1,0,1,12)
astsa::sarima(as.vector(t(respy1)),4,1,3,1,1,1,12)
astsa::sarima(as.vector(t(respy1)),4,1,3,2,1,1,12)
