#!/usr/bin/env Rscript
#
# Copyright 2020 Andrea Trentini (http://atrent.it)
#
#    This program is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program.  If not, see <http://www.gnu.org/licenses/>.*

#options(error=traceback) 
options(readr.num_columns = 0) # toglie messaggio in read_csv

source("./include.r")

suppressMessages(library(corrgram))

DIRGRAFICIC <- "GraficiCOVID"

covid=read_csv("COVID/pcm-dcp/dati-province/dpc-covid19-ita-province.csv")

# prendere i dati del COVID per provincia (generando anche il differenziale giorno per giorno)
# mediare i dati aria per inquinante e per provincia
# confrontare (plot, correlazione, qq) shiftando di 14 giorni (incubazione)

genGrafici = function(provincia,inquinante){
	covid_provincia=covid[covid$denominazione_provincia==provincia,]
	covid_provincia=covid_provincia[order(covid_provincia$data),]
	last_covid=max(covid_provincia$data)
	quanticovid=length(covid_provincia$data)
	#covid_lodi$data=anytime(covid_lodi$data)
	#covid_lodi$data=dmy_hms(covid_lodi$data)
	summary(covid_provincia) # non produce nulla!!!
	print(covid_provincia)

	ptrn=paste("^",provincia,".*",inquinante,".*\\.csv$",sep="")
	CSVs <- dir(pattern = ptrn,include.dirs=TRUE,recursive=TRUE)
	if(length(CSVs)==0){
		print("CSV mancanti")
		return()
	}
	print(CSVs)

	# carico i CSV, aggrego con media, prendo solo i dati post COVID
	aria=lapply(CSVs,read_csv) # genera n tibbles
	aria <- do.call(rbind, aria) # collassa in un solo tibble
	aria$Data <- dmy_hms(aria$Data)
	aria=aria[aria$Data>=COVIDPRE & aria$Data<=as.Date(last_covid)-14,] # inutile prendere dati prima, prendo fino a 14 giorni fa (assume che i dati COVID siano aggiornati a oggi)
	aria=aria[order(aria$Data),]
	if(length(aria$Valore)==0){ 
		print("Dati aria vuoto")
		return()
	}
	#print(aria)

	#aria <- aggregate(aria, by = list(as.Date(aria$Data, format = "%Y/%m/%d")), mean)
	aria <- aggregate(aria, by = list(as.Date(aria$Data)), mean)

	aria$Data <- NULL
	colnames(aria)[1] <- "Data"
	if(length(aria$Valore) < quanticovid){ 
		print("Dati aria con buchi")
		return()
	}
	if(length(aria$Valore) > quanticovid){ 
		print("Troppi dati aria?!? (dopo aggregate)")
		return()
	}
	print(aria)
	
	# calcolo incrementale
	covid_provincia$nuovi_casi <- c(0,diff(covid_provincia$totale_casi))

	jpeg(paste(DIRGRAFICIC,"/",provincia,"-",inquinante,".jpg",sep=""), width = 20, height = 15, units = "cm", res = 200)
	par(mfrow = c(2, 1))

	#plot(covid_provincia$totale_casi~covid_provincia$data,type="l",main="totale")
	plot(covid_provincia$nuovi_casi ~ covid_provincia$data,main="Nuovi casi giorno per giorno",type="l",xlab="data",ylab="nuovi casi")
	regressione <- lm(covid_provincia$nuovi_casi ~ covid_provincia$data)
	#cat("$$$ coefficiente trend line: ", regressione$coefficients[2], "\n")
	abline(regressione, col = "orange", lwd = 2)

	#qqplot(aria$Valore,covid_provincia$nuovi_casi,main="qqplot")
	#abline(0, 1, col = 'red')

	plot(aria$Valore~aria$Data,type="l",main=paste(inquinante,"(anticipato 14gg rispetto a 'nuovi casi')"),xlab="data",ylab="livello inquinante")
	#plot(covid_provincia$nuovi_casi,aria$Valore) #scatter
	regressione <- lm(aria$Valore ~ aria$Data)
	#cat("$$$ coefficiente trend line: ", regressione$coefficients[2], "\n")
	abline(regressione, col = "orange", lwd = 2)
	
	correlazione = cor(covid_provincia$nuovi_casi,aria$Valore)
	print(paste("correlazione:",correlazione))

	#legend("topright", paste("indice correlazione:", correlazione), text.col = "red")
	mtext(paste("indice correlazione:", correlazione),3)

	title(paste("[ COVID vs.",inquinante, "@", provincia,"]"), outer = TRUE ,line = -1)
	dev.off()
}

#########################################
# prime prove
#genGrafici("Lodi","PM10SM2005")
#genGrafici("Brescia","BiossidodiZolfo")
#genGrafici("Milano","PM10SM2005")
#stop()
#########################################

for(provincia in unique(covid[covid$denominazione_regione == "Lombardia",]$denominazione_provincia)){
	print(provincia)
	for(inquinante in c("ParticellesospesePM25","ParticolatoTotaleSospeso","PM10","PM10SM2005")){
		print(inquinante)
		genGrafici(provincia,inquinante)
	}
}