Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
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))]