################################################################################
# Course SMACH - January 2014
# project - Script to read weather data for climate change projection
# version $Date:: 2014-01-13 0#$ by $Author: brun $(v0: 2013-10-9)
################################################################################
library(maps)
library(reshape)
library(ggplot2)
#setwd("D:/CourseSmach/data/weather/DRIAS")
setwd("D:/data/modelia_formation/Volterra2016/data/weather/DRIAS")
fichier ="tasmintasmaxtasprwsmax_metro_IPSL_WRF_histo_QT_REF_19710101-20051231.txt"
weather_drias1 = read.table(file = fichier,skip=54,sep=";",stringsAsFactor=F)
fichier ="tasmintasmaxtasprwsmax_metro_IPSL_WRF_rcp4.5_QT_RCP4.5_20060101-21001231.txt"
weather_drias2 = read.table(file = fichier,skip=54,sep=";",stringsAsFactor=F)
weather_drias=rbind(cbind(period="histo",weather_drias1),cbind(period="rcp45",weather_drias2))
names(weather_drias)= c("period","Latitude","Longitude","Date","tasmin","tasmax","tas","pr","wsmax")
rm(fichier,weather_drias1,weather_drias2)

head(weather_drias)
attr(weather_drias, "nom_variable")= c("Date : Date du jour sous la forme 'AAAAMMJJ'",
                                       "Latitude, Longitude : position du point de grille SAFRAN en degres decimaux (WGS84)",
                                       "tasmin : Temperature minimale journaliere a 2 m [K]",
                                       "tasmax : Temperature maximale journaliere a 2 m [K]",
                                       "tas : Temperature moyenne a 2 m [K]",
                                       "pr : Precipitations totales [mm jr-1]",
                                       "wsmax : Vent maximal sans rafales [m/s]")

weather_drias$date=as.Date(paste(weather_drias$Date),"%Y%m%d")
weather_drias$year=as.numeric(format(weather_drias$date,"%Y"))
weather_drias$Tmin=weather_drias$tasmin-273.15
weather_drias$Tmax=weather_drias$tasmax-273.15
weather_drias$Tmean=(weather_drias$Tmin+weather_drias$Tmax)/2


weather_drias$site=NA
weather_drias[weather_drias$Latitude==43.5677&weather_drias$Longitude==1.49670,"site"]="Toulouse"
weather_drias[weather_drias$Latitude==45.7731&weather_drias$Longitude==4.85737,"site"]="Lyon"
weather_drias[weather_drias$Latitude==49.3266&weather_drias$Longitude==2.17227,"site"]="Paris"

table_site = unique(weather_drias[,c("site","Latitude","Longitude")])
map("france")
points(table_site$Longitude,table_site$Latitude,col="red",cex=2,pch=16)
################################################################################
# to complete....
sweather=subset(weather_drias,site=="Toulouse")
synth=cast(weather_drias,period+Latitude+Longitude+year~.,mean, value ="Tmean")
names(synth)[names(synth)=="(all)"]="Tmean_annual"
plot(synth$year,synth$Tmean_annual,col=synth$period)
# with ggplot
c <- ggplot(synth, aes(year, Tmean_annual,colour=factor(period)))
c + stat_smooth(method=lm, aes(fill = factor(period)))+ geom_point()
# end of file