1 MENDOZA station vs cmorph (data distribution)

station_vs_cmorph_mza=read.csv(
  file="/home/hadoop/paula/experiments/experiment3-cmorph-correction/station-vs-cmorph-mza-daily.csv",sep=',',header=F)
station_vs_cmorph_mza=cbind(rep("mza"),station_vs_cmorph_mza)
names(station_vs_cmorph_mza)<-c("city","ts","station","cmorph")

station_vs_cmorph_caba=read.csv(
  file="/home/hadoop/paula/experiments/experiment3-cmorph-correction/station-vs-cmorph-caba-daily.csv",sep=',',header=F)
station_vs_cmorph_caba=cbind(rep("caba"),station_vs_cmorph_caba)
names(station_vs_cmorph_caba)<-c("city","ts","station","cmorph")

station_vs_cmorph=rbind(station_vs_cmorph_caba,station_vs_cmorph_mza)
names(station_vs_cmorph)<-c("city","ts","station","cmorph")

month=month(as.Date(as.character(station_vs_cmorph$ts),"%Y%m%d"))
year=year(as.Date(as.character(station_vs_cmorph$ts),"%Y%m%d"))
station_vs_cmorph=cbind(station_vs_cmorph,year,month)
densityplot(~cmorph+station|city,data=station_vs_cmorph,plot.points=F,alpha=0.5,auto.key = T)

histogram(~cmorph+station|city,data=station_vs_cmorph,plot.points=F,alpha=0.5,auto.key = T,breaks=40)

1.1 Scatter Plot and Linear Regression for MZA

mza_data=station_vs_cmorph[which(station_vs_cmorph$city=="mza" ) ,3:4]
lm_mza=lm(mza_data)
print(lm_mza)
## 
## Call:
## lm(formula = mza_data)
## 
## Coefficients:
## (Intercept)       cmorph  
##      1.9203       0.1392
plot(mza_data,main="linear regression for MZA",col='skyblue')
abline(lm_mza,col='orange')

1.2 Scatter Plot and Linear Regression for CABA

caba_data=station_vs_cmorph[which(station_vs_cmorph$city=="caba" ) ,3:4]
lm_caba=lm(caba_data)
print(lm_caba)
## 
## Call:
## lm(formula = caba_data)
## 
## Coefficients:
## (Intercept)       cmorph  
##      5.7510       0.1243
plot(caba_data,main="linear regression for CABA",col='skyblue')
abline(lm_caba,col='orange')

1.3 Time series for MZA

plot(rollapply(mza_data$station,10,mean,by=10),type='l',col='red',ylab='precip',main="Station vs CMORPH cummulative events for MZA (avg 10 data points)",xaxt='n')
lines(rollapply(mza_data$cmorph,10,mean,by=10),col='blue')
legend("topleft",legend=c("station","cmorph"),col=c("red","blue"),cex=0.8,bty="n",lty=1)

1.4 Time series for CABA

plot(rollapply(caba_data$station,30,mean,by=30),type='l',col='red',ylab='precip',main="Station vs CMORPH cummulative events for CABA (avg 30 data points)",xaxt='n')
lines(rollapply(caba_data$cmorph,30,mean,by=30),col='blue')
 legend("topleft",legend=c("station","cmorph"),col=c("red","blue"),cex=0.8,bty="n",lty=1)

2 Anual precipitation volume for year 2010 (MZA and CABA)

vol_precip=station_vs_cmorph %>% group_by(city,year) %>% summarise(station=sum(station),cmorph=sum(cmorph))
barchart(station+cmorph~factor(year)|city,data=vol_precip, auto.key=T,scales=list(x=list(rot=45)))

3 Precipitation days by year for CABA and MZA (using station and cmorph data)

days=station_vs_cmorph %>% group_by(city,year) %>% summarise(n=n())
days_precip=station_vs_cmorph %>% group_by(city,year) %>% filter( station>0 | cmorph >0) %>% summarise(cmorph=sum(cmorph>0),station=sum(station>0))
days_precip=inner_join(days,days_precip,by=c("city","year"))
names(days_precip)<-c("city","year","total_days","cmorph_precip_days","station_precip_days")
barchart(station_precip_days+cmorph_precip_days+total_days~factor(year)|city,data=days_precip,auto.key=T,ylab="Days",scales=list(x=list(rot=45)))

4 Precipitation days by month for CABA and MZA (using station and cmorph data)

days=station_vs_cmorph %>% group_by(city,month) %>% summarise(n=n())
days_precip=station_vs_cmorph %>% group_by(city,month) %>% filter( station>0 | cmorph >0) %>% summarise(cmorph=sum(cmorph>0),station=sum(station>0))
days_precip=inner_join(days,days_precip,by=c("city","month"))
names(days_precip)<-c("city","month","total_days","cmorph_precip_days","station_precip_days")
barchart(station_precip_days+cmorph_precip_days+total_days~factor(month)|city,data=days_precip,auto.key=T,ylab="Days",xlab="Months",
         scales=list(x=list(rot=45,at=seq(1,12),label=mymonths)))

5 Building up coordinates CMORPH files for MZA

dir="/home/hadoop/paula/data/2.selected/cmorph_mza_daily/"
cmorph_files=list.files(dir,pattern="[0-9].*")
cmorph_coordenates = do.call(cbind, lapply(paste(dir,cmorph_files,sep=""),
              function(x) {
                        cmorph_coordinate=read.csv(x, stringsAsFactors = FALSE, sep=" ",header=F) 
                        cmorph_coordinate=cmorph_coordinate[,c(2,4)]
                        names(cmorph_coordinate)<-c("ts",basename(x))
                        return(cmorph_coordinate)
              }
              ))
cmorph_coordenates=cmorph_coordenates[,c(1,seq(2,ncol(cmorph_coordenates),2))]
cmorph_coordenates=melt(cmorph_coordenates,id.vars="ts",variable.name="coordinates",value.name="precip")
cmorph_coordenates=cmorph_coordenates %>% mutate(year=year(as.Date(as.character(cmorph_coordenates$ts),"%Y%m%d")),month=month(as.Date(as.character(cmorph_coordenates$ts),"%Y%m%d")))

5.1 Average precipitation volume (and Standard Deviation) for 10x10 CMORPH coordinates

coordinate_stats=cmorph_coordenates %>% group_by(coordinates) %>% summarise(mean=mean(precip),sd=sd(precip))
coordinate_mean=matrix(coordinate_stats$mean,nrow=10,ncol =10)
coordinate_sd=matrix(coordinate_stats$sd,nrow=10,ncol =10)
coordinate_days=cmorph_coordenates %>% group_by(coordinates) %>% summarise(days=sum(precip>0)/sum(precip>=0))


# Loading CMOPRH georeferenciated information for MZA
mza_geo_data=read.csv("/home/harpo/Dropbox/ongoing-work/meteo/scripts/R/cmorph_mza_data_geo.csv",header=T,sep=",")
mza_geo_data$lon=mza_geo_data$lon -360
mza_geo_data=cbind(mza_geo_data,coordinate_stats[,2:3])
mza_geo_data=cbind(mza_geo_data,coordinate_days[,2])

mza_map <- get_map(location = "mendoza", zoom = 7,maptype = "roadmap",color = "bw")

5.1.1 Average

g<-ggmap(mza_map)+geom_tile(data=mza_geo_data,aes(x=lon,y=lat,fill=mean),
                         alpha=0.8,size=0.1,colour='white')+ scale_fill_gradient(low = "white",high = "steelblue")
g

#d3heatmap::d3heatmap(coordinate_mean,symm=T,Colv=NA,Rowv=NA,color="Reds")

5.1.2 Standard Deviation

#d3heatmap::d3heatmap(coordinate_mean,symm=T,Colv=NA,Rowv=NA,color="Reds")

g<-ggmap(mza_map)+geom_tile(data=mza_geo_data,aes(x=lon,y=lat,fill=sd),
                         alpha=0.8,size=0.1,colour='white')+ scale_fill_gradient(low = "white",high = "steelblue")
g

5.2 Percentage of precipitation days for 10x10 CMORPH coordinates

g<-ggmap(mza_map)+geom_tile(data=mza_geo_data,aes(x=lon,y=lat,fill=days),
                         alpha=0.6,colour='white')+ scale_fill_gradient(low = "white",high = "steelblue")
g

#coordinate_days=matrix(coordinate_days$days,nrow=10,ncol =10)
#d3heatmap::d3heatmap(coordinate_days,symm=T,Colv=NA,Rowv=NA,color="Reds")

6 Influence of filtering low precipitation values on the CMOPRH cumulative precipitation

a=read.delim(file="/home/hadoop/paula/data/8.cmorph-bilinear/8.reg-mza-annual-s1-p6.csv", sep=",", header=T)
x=a   %>% group_by(target_year) %>% summarise(precip=sum(precip_3hr_bilinear))
y=a   %>% filter(precip_3hr_bilinear>0.1) %>% group_by(target_year) %>% summarise(precip=sum(precip_3hr_bilinear))
z=rbind(cbind(filter=rep("normal"),x),cbind(filter=rep("filtered"),y))
 barchart(precip~factor(target_year)|filter,data=z)