library("jsonlite")
library("sf")
library("raster")

data=jsonlite::fromJSON(txt = "/home/roux/D/PROJETS/LMI_SENTINELA/SITE_SENTINELLE_GF-AP/DONNEES_HARMO/23-10-2023/consultation_harmonized.json")


harmonized_table <- data.table::as.data.table(data$data)

loc <- read_sf("/home/roux/D/PROJETS/LMI_SENTINELA/SITE_SENTINELLE_GF-AP/REFERENTIELS_CARTO/DONNEES/LOCALITES/Spatial_referential_2019_11_05/Spatial_referential_2019_11_05.shp")
risk <- raster("/home/roux/D/PROJETS/PROGYSAT/LIZMAP/DATA/RASTERS_RISQUE/risque_expo.tif")
hfootprint <- raster("/home/roux/D/PROJETS/PROGYSAT/LIZMAP/DATA/RASTERS_INDICATEURS/eh_ag.tif")
risk_ind <- raster("/home/roux/D/PROJETS/PROGYSAT/test_ind.tif")
names(risk_ind)<- "risque_expo"
values(risk_ind) <- (values(risk_ind)-min(values(risk_ind), na.rm=TRUE))/(max(values(risk_ind), na.rm=TRUE)-min(values(risk_ind), na.rm=TRUE))

cases <- harmonized_table[harmonized_table$active_diagnosis == TRUE & harmonized_table$new_attack == TRUE]
n_cases_per_infect_loc <- tapply(cases$id_consultation, INDEX=cases$infection_place, FUN=length)
n_cases_per_infect_loc <- data.frame(ID_LOC=names(n_cases_per_infect_loc), ncases=n_cases_per_infect_loc)

test=merge(n_cases_per_infect_loc, loc, by="ID_LOC")
test=st_as_sf(test)
plot(test["ncases"])

graphics.off()
rsquared=c()
radjsquared=c()
for (buff in c(250, 500, 1000, 2000)){
  extraction <- extract(risk, test, buffer=buff, fun=median, na.rm=TRUE, sp=TRUE)
  
  step=0.025
  x=seq(0,1-step,step)+step/2
  y=tapply(log(extraction$ncases),INDEX=cut(extraction$risque_expo, breaks=seq(0,1,step)),FUN=min)
  x11(); plot(seq(0,1-step,step)+step/2, y, main=as.character(buff))
  mod=lm(y~0+x)
  smod=summary(mod)
  x11(); plot(extraction$risque_expo, log(extraction$ncases), main=as.character(buff))

  rsquared=c(rsquared,  smod$r.squared) 
  radjsquared=c(radjsquared,  smod$adj.r.squared) 
  
  
  # extraction_all_loc <- extract(risk, loc, buffer=buff, fun=mean, na.rm=TRUE, sp=TRUE)
  # #extraction <- extract(risk_ind, test, buffer=buff, fun=mean, na.rm=TRUE, sp=TRUE)
  # #extraction_all_loc <- extract(risk_ind, loc, buffer=buff, fun=mean, na.rm=TRUE, sp=TRUE)
  # select=extraction_all_loc[extraction_all_loc$ID_LOC%in%setd,]
  # 
  # w=extraction$ncases
  # w[w<=1] <- 0
  # w[w>1] <- 1 
  # 
  # x11()
  # par(mfrow=c(1,2))
  # h<-hist(rep(extraction$risque_expo, w, each=TRUE), freq=FALSE, ylim=c(0,10), breaks=seq(0,1,0.1), main=as.character(buff))
  # hh<-hist(select$risque_expo, breaks=h$breaks, freq=FALSE, ylim=c(0,10))

}
extraction$ID_LOC[c(which(extraction$risque_expo>0.55))]