grep "module load" /u/arpa/bonafeg/src/carminio/R/*R
module load gdal proj R

new_path <- "/u/arpa/bonafeg/R/x86_64-pc-linux-gnu-library/3.3"
.libPaths(unique(c(new_path, .libPaths())))

source("/u/arpa/bonafeg/src/carminio/R/get_rds.R")
rds_data <- get_rds(from = "2017-01-05 00", to = "2017-01-08 23")

library(ggplot2)
ggplot(rds_data, aes(x=Temp_C, y=Press,
                     col=as.character(Time_launch))) +
  geom_point(size=0.2) +
  scale_color_brewer(name="launch time", palette = "Set2") +
  scale_y_reverse() +
  theme_bw()

ggplot(rds_data, aes(x=WindSpeed, y=H_amsl,
                     col=as.character(Time_launch))) +
  geom_point(size=0.2) +
  scale_color_brewer(name="launch time", palette = "Set2") +
  scale_y_log10() +
  theme_bw()

library(dplyr)
rds_data %>% filter(H_amsl<=3000) %>%
  group_by(Time_launch) %>%
  summarize(Nlev=n()) %>%
  summarize(Nlev_ave=mean(Nlev))

rds_data <- get_rds(from = "2017-01-01 00",
                    to = "2017-01-31 23",
                    nrows=600)

rds_data %>% mutate(H_agl=H_amsl-H_amsl_station) %>%
  group_by(Time_launch) %>%
  summarize(T_sfc=first(Temp_C),
            T_max=max(Temp_C),
            h_T_max=H_agl[which.max(Temp_C)],
            dT=T_max-T_sfc) %>%
  ggplot(aes(x=Time_launch)) +
    geom_point(aes(y=h_T_max),col="navy") +
    scale_y_log10(limits=c(1,500),
                  breaks=c(10,20,50,100,200,500),
                  name="h_Tmax",
                  sec.axis = dup_axis(breaks=c(2,5,10),
                                      name="dT")) +
    geom_point(aes(y=dT), col="sienna") +
    theme_bw()

source("/u/arpa/bonafeg/src/carminio/R/read_bronx.R")
NO2_obs <- get_hourly(stations=c("OSV","CAI"), year=2016, 
                      yearly=T, pollutant="NO2")
NOx_obs <- get_hourly(stations=c("OSV","CAI"), year=2016, 
                      yearly=T, pollutant="NOX")
bind_rows(NO2_obs, NOx_obs) %>%
  transmute(Point=LOC_CODE,
      Pollutant=INQ_CODE,
      Value=VALUE,
      Time=DATEHOUR,
      X=UTM_X_33, Y=UTM_Y_33) -> aq_obs

library(openair)
library(tidyr)
linearRelation(aq_obs %>% mutate(date=Time) %>%
        filter(Point=="OSV") %>%
        spread(key=Pollutant, value=Value),
        x="NOX", y="NO2")

timeVariation(aq_obs %>% mutate(date=Time),
        pollutant="Value", group="Point", type="Pollutant")

source("/u/arpa/bonafeg/src/carminio/R/get_osmer_syn.R")
read_osmer_anag()$code
met_obs <- get_osmer(year=2016, st=c("UDI","TRI","ENE"))

windRose(met_obs,ws="WS.ist.10m",wd="WD.ist.10m",type="name")
windRose(met_obs %>% filter(code=="UDI") %>% 
           mutate(date=time_UTC),
        ws="WS.ist.10m",wd="WD.ist.10m",type="season")

tz(aq_obs$Time)
tz(met_obs$time_UTC)
str(aq_obs)
str(met_obs)
aq_obs %>% mutate(endTime=Time+59*60) %>%
  filter(Pollutant=="NO2") -> aq_obs
met_obs %>% rename(endTime=time_UTC) %>%
  filter(code=="UDI") -> met_obs
aq_met <- left_join(aq_obs, met_obs)
polarPlot(aq_met,
        pollutant="Value",
        x="WS.ist.10m",
        wd="WD.ist.10m",
        type="Point")

source("/u/arpa/bonafeg/src/carminio/R/read_kriging.R")
files <- paste0("O3_",2015:2017,
                "_CivilYear_8hRunningAvgs-DailyMaxes-Nexc120",
                "_DLGS155_02_2011_ukriging_sup.asc")
read_kriging_year(year="2015",file=files[1])->o3_2015
read_kriging_year(year="2016",file=files[2])->o3_2016
read_kriging_year(year="2017",file=files[3])->o3_2017

library(RColorBrewer)
display.brewer.all()
plot(o3_2017, col=brewer.pal(12,"Spectral")[11:2])
source("/u/arpa/bonafeg/src/rutilante/R/draw.marks.R")
draw.marks(file="/u/arpa/bonafeg/src/rutilante/dat/alcune_localita.dat",
           xyz=c(1,2,3),xy.fact=1,cex=0.7)

str(o3_2015)
str(o3_2015@data)
e2015 <- ecdf(o3_2015@data[[1]])
e2016 <- ecdf(o3_2016@data[[1]])
e2017 <- ecdf(o3_2017@data[[1]])
plot(e2015, col="darkred", lwd=2,
  main="Empirical Cumulated Distribution Function",
  xlab="Exceedances (days)")
plot(e2016, col="olivedrab", add=T, lwd=2)
plot(e2017, col="orange", add=T, lwd=2)
legend("bottomright",lwd=2,col=c("darkred","olivedrab","orange"), 
       legend=2015:2017)
abline(v=25, col="red", lty=2)

source("/u/arpa/bonafeg/src/carminio/R/get_fair.R")
met_obs <- get_osmer(year=2016, st="UDI")
met_gap <- extract_model(from="2016-01-01", to="2016-01-31",
                model="GAP", type="", run="y2016",
                v2=NULL, v3=c("U","V"),
                lat=met_obs$lat[1], lon=met_obs$lon[1],
                code="UDI")

str(met_gap)
source("/u/arpa/bonafeg/src/rutilante/R/meteo.functions.R")
met_gap <- spread(met_gap, key="Var", value="Value")
met_gap %>% mutate(WS=ff(U,V), WD=dd(U,V)) -> met_gap

met_Jan <- bind_rows(
        met_obs %>%
        filter(time_UTC>=as.POSIXct("2016-01-01"),
               time_UTC<=as.POSIXct("2016-01-31")) %>%
        transmute(date=time_UTC+60,
                WS=WS.ist.10m, WD=WD.ist.10m,
                Source="obs"),
        met_gap %>%
        transmute(date=Time,
                WS=WS,WD=WD,
                Source="mod"))
windRose(met_Jan, ws="WS",wd="WD",type="Source")

eobs <- ecdf(met_Jan$WS[met_Jan$Source=="obs"])
emod <- ecdf(met_Jan$WS[met_Jan$Source=="mod"])
plot(emod, col="darkred", pch=19,
       main="Empirical Cumulated Distribution Function",
       xlab="wind speed (m/s)")
plot(eobs, col="olivedrab", add=T, pch=19)
legend("bottomright",pch=19,col=c("darkred","olivedrab"), 
       legend=c("mod","obs"))
abline(v=2, lty=2, col="grey")

met_Jan.2 <- spread(met_Jan %>% select(-WD), 
                    key=Source, value=WS)
TaylorDiagram(met_Jan.2, obs = "obs", mod = "mod")