RaddictsParis/M7_outilsCartos

## M7_outilsCartos
library("ggmap")
library("stringr")
library("plotGoogleMaps")
library("pgirmess")
library("png")
library("proto")
library("gridExtra")
library("spdep")
library("RColorBrewer")
library("rMaps")
library("leafletR")
library("rCharts")
library("rgdal")

# Fonctions auxiliaires ###################

geom_custom <- function (mapping = NULL, data = NULL, stat = "identity", position = "identity", show_guide = FALSE, ...) {
  GeomCustom$new(mapping = mapping, data = data, stat = stat,
                 position = position, show_guide = show_guide, ...)
}
GeomCustom <- proto(ggplot2:::Geom, {

  objname <- "custom"

  new <- function(., data = NULL, mapping = NULL, grob =NULL, ...) {
    .super$new(., data = data, mapping = mapping, grob=grob, inherit.aes = FALSE, ...)
  }

  draw <- function(., data, scales, coordinates, grob, just=c("centre", "centre"), ...) {
    with(coord_transform(coordinates, data, scales), {
      width <- grobWidth(grob)
      height <- grobHeight(grob)
      grob <- editGrob(grob, vp=viewport(x = x, y = y,
                                         width=width, height=height, just = just), ...)
      ggname(.$my_name(), grob)
    })
  }

  default_aes <- function(.)
    aes(x=0.5, y=0.5)
  default_stat <- function(.) StatIdentity

})
PNG.to.rasterGrob<-function(Width=0.6,Height=0.6,Unit="cm",path){
  NAME<-toupper(gsub("^.*/(.*).png$", "\\1", path))
  assign(paste('i',NAME,sep=''),readPNG(path),inherits = TRUE)
  assign(paste('g',NAME,sep=''),rasterGrob(eval(parse(text=paste('i',NAME,sep=''))),width = unit(Width, Unit), height = unit(Height, Unit), interpolate=TRUE),inherits = TRUE)
  print(paste("L'objet du type raster s'appel:",paste('g',NAME,sep=''),sep=' '))
}

# Début script ###################

##PROJECTIONS
LIIE<-"+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs"
WGS84<-"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

##GADM ADMINISTRATIVE BOUNDARIES
load(url("http://biogeo.ucdavis.edu/data/gadm2/R/FRA_adm5.RData"))
cible<-geocode("25 Rue du Petit Musc 75004 Paris",output="latlon")
scible<-SpatialPointsDataFrame(coords=cible, data=cible, coords.nrs = numeric(0), proj4string = CRS(WGS84), match.ID = TRUE)
id_entite_geo<-over(scible,gadm)$ID_5
entite_geo<-subset(gadm,ID_5%in%id_entite_geo)

##OSM LAYERS
path.osm<-"http://download.geofabrik.de/europe/france/ile-de-france-latest.shp.zip"
path.out<-"YOUR_PATH"
download.file(path.osm, "parisosm.zip")
dir.create(path.out)
unzip("parisosm.zip")
unzip("parisosm.zip",exdir=path.out)
files<-list.files(path.out,pattern=".shp")
ifiles<-substr(files,1,str_locate(files,".shp")[,1]-1)
lapply(1:length(ifiles),function(i){assign(toupper(ifiles[i]),readOGR(dsn=path.out,layer=ifiles[i]),inherits=TRUE)})

##LA GRILLE
NB_SIDE_BLOCS<-10
start_point<-spTransform(scible,CRS(LIIE))
SIDE<-NB_SIDE_BLOCS/2
coord<-matrix(c(coordinates(start_point)[1],coordinates(start_point)[2]+100*SIDE,coordinates(start_point)[1],coordinates(start_point)[2]-100*SIDE),nr=2,byrow=TRUE)
TEMP1<-pave(coord,round(SIDE*1.25,0),SIDE,output="spdf")
TEMP2<-pave(coord,round(SIDE*1.25,0),SIDE,output="spdf",ydown=FALSE)
TEMP1<-spChFIDs(TEMP1,as.character(paste("S1",rownames(TEMP1@data),sep="")))
TEMP2<-spChFIDs(TEMP2,as.character(paste("S2",rownames(TEMP2@data),sep="")))
BLOCS<-spRbind(TEMP1,TEMP2)
co<-coordinates(BLOCS)
BLOCS@data$ID2<-1:length(BLOCS)
BLOCS<-spChFIDs(BLOCS,as.character(BLOCS@data$ID2))
proj4string(BLOCS)<-CRS(LIIE)
BLOCS<-spTransform(BLOCS,CRS(proj4string(LANDUSE)))
UBLOCS<-unionSpatialPolygons(BLOCS,rep(1,length(BLOCS)))

##OVER
oLANDUSE<-over(LANDUSE,UBLOCS)
oLANDUSE<-LANDUSE[which(!is.na(oLANDUSE)),]
oNATURAL<-over(NATURAL,UBLOCS)
oNATURAL<-NATURAL[which(!is.na(oNATURAL)),]
oPLACES<-over(PLACES,UBLOCS)
oPLACES<-PLACES[which(!is.na(oPLACES)),]
oPOINTS<-over(POINTS,UBLOCS)
oPOINTS<-POINTS[which(!is.na(oPOINTS)),]
oPOINTS2<-over(POINTS,UBLOCS_BIS)
oPOINTS2<-POINTS[which(!is.na(oPOINTS2)),]
oRAILWAYS<-over(RAILWAYS,UBLOCS)
oRAILWAYS<-RAILWAYS[which(!is.na(oRAILWAYS)),]
oROADS<-over(ROADS,UBLOCS)
oROADS<-ROADS[which(!is.na(oROADS))]

##SUBSET DES ÉLÉMENTS À AFFICHER
BASIN<-subset(oLANDUSE,type=="basin")
CAFE<-subset(oPOINTS2,type=="cafe")
GRASS<-subset(oLANDUSE,type=="grass")
HOSPITALS<-subset(oPOINTS,type=="hospital")
INTERNET<-subset(oPOINTS,type=="internet_cafe")
METRO<-subset(RAILWAYS,type=="subway")
METRO1<-subset(METRO,name=="Métro 1")
METRO5<-subset(METRO,name=="Métro 5")
METRO7<-subset(METRO,name=="Métro 7")
METRO8<-subset(METRO,name=="Métro 8")
PARK<-subset(oNATURAL,type=="park")
PUB<-subset(oPOINTS2,type=="pub")
RIVERBANK<-subset(oNATURAL,type=="riverbank")
STATIONS<-subset(oPOINTS,type=="station")
ST1<-STATIONS[which(STATIONS@data$name=="Bastille"),][2,]
ST2<-STATIONS[which(STATIONS@data$name %in% c("Pont Marie","Sully - Morland","Saint-Paul")),]
NSTATIONS<-spRbind(ST1,ST2)
WATER<-subset(oNATURAL,type=="water")

#PLOTS DIAPOSITIVES 7-16
plot(UBLOCS,col="grey",lwd=5,border="grey")
plot(RIVERBANK,add=TRUE,col="#9bbffe",lwd=1,border="#8db7fd")
plot(WATER,add=TRUE,col="#9bbffe",lwd=1,border="#8db7fd")
plot(BASIN,add=TRUE,col="#9bbffe",lwd=1,border="#8db7fd")
plot(PARK,add=TRUE,col="#cadfaa",lwd=0.5)
plot(GRASS,add=TRUE,col="#cadfaa",lwd=0.5)
plot(OROADS,add=T,col="white",border="black")
# Places (sans Paris et les arrondissements)
PLACES_clean = subset(PLACES, !grepl("Arrondissement", PLACES@data$name))
PLACES_clean = subset(PLACES_clean, PLACES_clean@data$name != "Paris")
coordPLACESdecalees = coordinates(PLACES_clean)
coordPLACESdecalees[,2] = coordPLACESdecalees[,2]+0.0005
plot(PLACES_clean,add=TRUE,col="brown",lwd=5,pch=0)
text(coordPLACESdecalees, offset = 20, labels=PLACES_clean@data$name, cex = 0.7, col = "brown")
plot(scible,add=TRUE,lwd=2,cex=2.5)
plot(scible,add=TRUE,pch=21,cex=1.7,lwd=2,bg="orange")
plot(METRO7,add=TRUE,col="#ff82b4",lwd=5)
plot(METRO1,add=TRUE,col="#ffcd09",lwd=5)
plot(METRO8,add=TRUE,col="#cfacd4",lwd=5)
plot(METRO5,add=TRUE,col="#f28e42",lwd=5)
plot(NSTATIONS,col="white",cex=1.1,pch=19,add=T)
plot(NSTATIONS,col="#1a3c90",cex=1.2,pch="O",add=T,lwd=2)
plot(NSTATIONS,col="#1a3c90",cex=0.5,lwd=5,pch="M",add=T)
plot(CAFE,col="white",cex=1.1,pch=19,add=T)
plot(CAFE,col="brown",cex=1.2,pch="O",add=T,lwd=2)
plot(CAFE,col="brown",cex=0.5,lwd=5,pch="C",add=T)
plot(PUB,col="white",cex=1.1,pch=19,add=T)
plot(PUB,col="#075933",cex=1.2,pch="O",add=T,lwd=2)
plot(PUB,col="#075933",cex=0.5,lwd=5,pch="P",add=T)
plot(INTERNET,col="white",cex=1.4,pch=19,add=T)
plot(INTERNET,col="orange",cex=1.5,pch="O",add=T,lwd=2)
plot(INTERNET,col="orange",cex=0.7,lwd=5,pch="IC",add=T)
plot(HOSPITALS,col="white",cex=1.6,border="black",pch=15,add=T)
plot(HOSPITALS,col="grey",cex=1.7,pch=0,add=T,lwd=2)
plot(HOSPITALS,col="red",cex=0.8,lwd=5,pch="H",add=T)

##PLOT DIAPOSITIVE 17
plot(UBLOCS)
plot(oROADS,add=TRUE,lwd=2.5)
plot(scible,add=TRUE,lwd=2,cex=5)
plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

##PLOT DIAPOSITIVE 18
plot(UBLOCS)
plot(oROADS,add=TRUE,lwd=2.5,col="gray")
plot(BLOCS,add=TRUE,border="firebrick",lwd=6)
plot(BLOCS,add=TRUE,border="yellow",lwd=1.5)
plot(scible,add=TRUE,lwd=2,cex=5)
plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

##SUBSET 4 BLOCS
sblocs<-subset(BLOCS,ID%in%c(3,4))
usblocs<-unionSpatialPolygons(sblocs,rep(1,length(sblocs)))

##ROADS dans 4 BLOCS
ov<-over(ROADS,usblocs)
w<-which(!is.na(ov))
oroads<-ROADS[w,]

##POINTS dans 4 BLOCS
ov<-over(POINTS,usblocs)
w<-which(!is.na(ov))
opoints<-POINTS[w,]
copoints<-data.frame(coordinates(opoints)) #coordinates pour la represenation graphique
colnames(copoints)<-c("lon","lat")

##PLOT DIAPOSITIVE 19
plot(usblocs)
plot(sblocs,add=TRUE,border="firebrick",lwd=10)
plot(sblocs,add=TRUE,border="yellow",lwd=3)
plot(oroads,add=TRUE,lwd=2.5)
plot(opoints,add=TRUE,cex=1.8,lwd=3)
plot(opoints,add=TRUE,pch=21,cex=1.8,lwd=3,bg="lightgreen")
plot(scible,add=TRUE,lwd=2,cex=3)
plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

##FONCTIONS AUXILIAIRES
source("YOUR_PATH/functions_raddicts.R")

uu<-1.2
path<-"YOUR_PATH/RADDICTS.png"
PNG.to.rasterGrob(Width=uu*4,Height=uu,Unit="cm",path)

##GGPLOTS DIAPOSITIVE 20,21,22 Il suffit de changer GMAP par GMAP1,GMAP2,GMAP3 dans la ligne 110
GMAP1 <- ggmap(get_map(location=bbox(usblocs),zoom=18, source="google",maptype = "roadmap"), extent = "panel", maprange=FALSE)
GMAP2 <- ggmap(get_map(location=bbox(usblocs),zoom=18, source="stamen",maptype = "toner"), extent = "panel", maprange=FALSE)
GMAP3 <- ggmap(get_map(location=bbox(usblocs),zoom=18, source="google",maptype = "hybrid"), extent = "panel", maprange=FALSE)

  GMAP +stat_density2d(data = copoints, aes(x = lon, y = lat,  fill = ..level..,col=NULL, alpha = ..level..),size = 0.01,h=0.0008, bins = 45, geom = 'polygon')+
  geom_density2d(data = copoints, aes(x = lon, y = lat,col="darkorange4"),h=0.0008)+
  scale_fill_gradient(low = 'yellow', high = 'firebrick')+
  scale_alpha(range = c(0.00, 0.75), guide = FALSE)+
  annotate('custom', x=cible$lon,y=cible$lat, grob=gRADDICTS)+
  theme(axis.text.x=element_blank(),axis.ticks=element_blank(),
        axis.title.x = element_blank(),axis.text.y=element_blank(),axis.title.y = element_blank())+
  theme(legend.position='bottom',legend.title = element_text(colour='black', size=14, face='bold'))+
  theme(legend.text = element_text(colour='black', size = 7, face = 'bold'))+
  theme(plot.title = element_text(size=12,lineheight=.8, face='bold'))+
  theme(plot.background = element_rect(fill='white ', colour='white '))+
  theme(panel.background=element_rect(fill='white ', colour='white '),
        panel.grid.major = element_blank(),panel.grid.minor = element_blank())+
  theme(panel.border=element_rect(colour = 'black',fill='transparent', size=6))+
  theme(legend.position = 'none', axis.title = element_blank(), text = element_text(size = 12))


##PLOTKML DIAPOSITIVE 23
kml(entite_geo1, ".", "zone_rdv.kml", colour = "yellow", altitude = 100, alpha = 0.5)
kml(scible, ".", "point_rdv.kml", colour = "red", width = 2, labels="R addicts Paris", LabelScale = 1.2, size = 1.2)

##RMAPS PAGES 24-26
opoints<-subset(OPOINTS,!is.na(name))
copo<-data.frame(coordinates(opoints))
colnames(copo)<-c("lon","lat")
copo$popup<-opoints@data$name
map = Leaflet$new()
map$setView(c(cible$lat, cible$lon ), 16)
PROVIDER<-"Stamen.Toner" #Autres Stamen.Watercolor (p.25), Esri.WorldImagery (p.26)
map$tileLayer(provider="Stamen.Watercolor")
for(i in sample(1:nrow(copo),20)){
  map$marker(c(copo$lat[i],copo$lon[i]), bindPopup = copo$popup[i])
}


##STRUCTURES DE VOISINAGE
nbpoly1<-poly2nb(BLOCS,queen=TRUE)

##PLOT DIAPOSITIVE 28
plot(BLOCS, border="grey")
plot(oROADS,add=TRUE,lwd=4,col="cornflowerblue")
plot(oROADS,add=TRUE,lwd=2,col="gray")
plot(nbpoly1, coordinates(BLOCS), add=TRUE,bg="gold",pch=21,lwd=2)
plot(UBLOCS, lwd=6,add=TRUE)
plot(scible,add=TRUE,lwd=2,cex=3)
plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

##PLOT DIAPOSITIVE 29
nbpoly2<-knn2nb(knearneigh(coordinates(BLOCS), k = 1), row.names = BLOCS@data$ID2)
plot(BLOCS, border="grey")
plot(oROADS,add=TRUE,lwd=4,col="cornflowerblue")
plot(oROADS,add=TRUE,lwd=2,col="gray")
plot(nbpoly2, coordinates(BLOCS), add=TRUE,bg="gold",pch=21,lwd=2)
plot(UBLOCS, lwd=6,add=TRUE)
plot(scible,add=TRUE,lwd=2,cex=3)
plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

##PLOT DIAPOSITIVE 30
nbpoly2<-knn2nb(knearneigh(coordinates(BLOCS), k = 2), row.names = BLOCS@data$ID2)
plot(BLOCS, border="grey")
plot(oROADS,add=TRUE,lwd=4,col="cornflowerblue")
plot(oROADS,add=TRUE,lwd=2,col="gray")
plot(nbpoly2, coordinates(BLOCS), add=TRUE,bg="gold",pch=21,lwd=2)
plot(UBLOCS, lwd=6,add=TRUE)
plot(scible,add=TRUE,lwd=2,cex=3)
plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")


##ENVIRONNEMENT SOURCE DIAPOSITIVE 31
ov<-over(POINTS,UBLOCS)
w<-which(!is.na(ov))
OPOINTS<-POINTS[w,]
ov<-over(OPOINTS,BLOCS)
ID<-names(table(ov$ID2))
DATA<-data.frame(cbind(ID=ID))
DATA$Y<-table(ov$ID2)
ov<-over(oROADS,BLOCS)
ID<-names(table(ov$ID2))
DATA2<-data.frame(cbind(ID=ID))
DATA2$X<-table(ov$ID2)
DATA<-merge(DATA,DATA2,by="ID");remove(DATA2) #Y nb points du bloc, X nb troncons de rue du bloc


##PLOT DIAPOSITIVE 31
plot(UBLOCS)
plot(oROADS,add=TRUE,lwd=2.5,col="gray")
plot(OPOINTS,add=TRUE,pch=21,bg="lightgreen")
plot(BLOCS,add=TRUE,lwd=6)
plot(scible,add=TRUE,lwd=2,cex=5)
plot(scible,add=TRUE,lwd=2,cex=3)
plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

##ENVIRONNEMENT SKATER
rownames(BLOCS@data)<-BLOCS@data$ID2
temp<-merge(BLOCS@data,DATA,by.x="ID2",by="ID",all.x=TRUE,all.y=FALSE)
rownames(temp)<-temp$ID2
BLOCS@data<-temp[rownames(BLOCS@data),]
BLOCS@data[is.na(BLOCS@data)]<-0
nbpoly1<-poly2nb(BLOCS,queen=FALSE)
lcosts <- nbcosts(nbpoly1, scale(BLOCS@data[,c("Y","X")])) #Definition des poids arbitraires a titre d exemple.
nb.w <- nb2listw(nbpoly1, lcosts, style="B")
mst.bh <- mstree(nb.w,5)

##PLOT DIAPOSITIVE 32
plot(mst.bh, coordinates(BLOCS), col="steelblue",lwd=6,cex.lab=0.01, cex.circles=0.035, fg="blue")
plot(BLOCS, border=gray(.5), add=TRUE,lws=3)
plot(UBLOCS,add=TRUE,lwd=6)
plot(scible,add=TRUE,lwd=2,cex=3)
plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

##FONCTION SKATER
res1 <- skater(edges=mst.bh[,1:2], data=scale(BLOCS@data[,c("Y","X")]),ncuts=5) #Si critere nb groups, a changer ncuts par 5, 15 et 18
res1 <- skater(edges=mst.bh[,1:2], data=scale(BLOCS@data[,c("Y","X")]),vec.crit=BLOCS@data$Y,crit=400) #Si critere variable quantitative, a changer crit par 300 et 400
skBLOCS<-unionSpatialPolygons(BLOCS,res1$groups) #regroupement de mailles geographiques

##PLOT DIAPOSITIVE 33-37
plot(skBLOCS,col=brewer.pal(length(unique(res1$groups)), "PuBu")) #autres palettes: Set3, Pastel1, Pastel2, Blues
plot(oROADS,add=TRUE,lwd=2.5,col="gray")
plot(BLOCS,add=TRUE,border="gray34")
plot(UBLOCS,lwd=6,add=TRUE)
plot(skBLOCS,lwd=6,add=TRUE)
plot(OPOINTS,add=TRUE,pch=21,cex=0.5,bg="lightgreen")


##ENVIRONNEMENT MODELES SPATIAUX
# Modèle linéaire simple (non-spatial)
BLOCS@data$X<-as.numeric(BLOCS@data$X)
BLOCS@data$Y<-as.numeric(BLOCS@data$Y)
lmmod=lm(Y~X,data=BLOCS@data)
summary(lmmod)
lm.morantest(lmmod, nb2listw(nbpoly1)) # rejet de H0 : absence d'autocorrélation spatiale dans les résidus
# conversion matrice voisinage
W.mat=nb2mat(nbpoly1)
W.listw=nb2listw(nbpoly1,style="W")
# modèle LAG par max de vraissemblance
lagmodel = lagsarlm(Y~X,data=BLOCS@data,listw=W.listw)
summary(lagmodel) # modèle significatif + rejet d'absence d'autocorrélation spatiale à 0.97% -> bonne prise en compte de l'autocorrélation spatiale
# modèle durbin par max de vraissemblance
durbin = lagsarlm(Y~X,data=BLOCS@data,listw=W.listw,type="mixed")
summary(durbin)
# modèle MESS
semmod = lagmess(Y~X,data=BLOCS@data,listw=W.listw)
summary(semmod)
# modèle SEM
semmod = spautolm(Y ~ X,data=BLOCS@data,listw=W.listw)
summary(semmod)
# comparaison de modèles
# méthode 1 : les tests
lm.LMtests(lmmod,listw=W.listw, test="all")
# méthode 2 : sélection sur critère
AIC(lmmod,lagmodel,durbin,semmod)
# impacts
impacts(lagmodel, listw=W.listw)
# erreur relative : LM vs. LAG
FIT<-predict(lagmodel,BLOCS@data,W.listw)
FIT<-data.frame(cbind(LAG=fitted(lagmodel),LM=fitted(lmmod),Y=BLOCS@data$Y))
head(FIT)
FIT$ELAG<-abs(FIT$Y-FIT$LAG)/FIT$Y
FIT$ELM<-abs(FIT$Y-FIT$LM)/FIT$Y
FIT$Y<-NULL
FIT[FIT>10^10]<-0
BLOCS@data<-cbind(BLOCS@data,FIT)
BLOCS@data$ERREUR_LAG<-BLOCS@data$ELAG
BLOCS@data$ERREUR_LM<-BLOCS@data$ELM

#SPPLOT DIAPOSITIVE 38
spplot(BLOCS,c("ERREUR_LAG","ERREUR_LM"),col.regions = brewer.pal(9, "PuBu"))
	library("ggmap")
	library("stringr")
	library("plotGoogleMaps")
	library("pgirmess")
	library("png")
	library("proto")
	library("gridExtra")
	library("spdep")
	library("RColorBrewer")
	library("rMaps")
	library("leafletR")
	library("rCharts")
	library("rgdal")

	# Fonctions auxiliaires ###################

	geom_custom <- function (mapping = NULL, data = NULL, stat = "identity", position = "identity", show_guide = FALSE, ...) {
	GeomCustom$new(mapping = mapping, data = data, stat = stat,
	position = position, show_guide = show_guide, ...)
	}
	GeomCustom <- proto(ggplot2:::Geom, {

	objname <- "custom"

	new <- function(., data = NULL, mapping = NULL, grob =NULL, ...) {
	.super$new(., data = data, mapping = mapping, grob=grob, inherit.aes = FALSE, ...)
	}

	draw <- function(., data, scales, coordinates, grob, just=c("centre", "centre"), ...) {
	with(coord_transform(coordinates, data, scales), {
	width <- grobWidth(grob)
	height <- grobHeight(grob)
	grob <- editGrob(grob, vp=viewport(x = x, y = y,
	width=width, height=height, just = just), ...)
	ggname(.$my_name(), grob)
	})
	}

	default_aes <- function(.)
	aes(x=0.5, y=0.5)
	default_stat <- function(.) StatIdentity

	})
	PNG.to.rasterGrob<-function(Width=0.6,Height=0.6,Unit="cm",path){
	NAME<-toupper(gsub("^./(.).png$", "\\1", path))
	assign(paste('i',NAME,sep=''),readPNG(path),inherits = TRUE)
	assign(paste('g',NAME,sep=''),rasterGrob(eval(parse(text=paste('i',NAME,sep=''))),width = unit(Width, Unit), height = unit(Height, Unit), interpolate=TRUE),inherits = TRUE)
	print(paste("L'objet du type raster s'appel:",paste('g',NAME,sep=''),sep=' '))
	}

	# Début script ###################

	##PROJECTIONS
	LIIE<-"+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs"
	WGS84<-"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

	##GADM ADMINISTRATIVE BOUNDARIES
	load(url("http://biogeo.ucdavis.edu/data/gadm2/R/FRA_adm5.RData"))
	cible<-geocode("25 Rue du Petit Musc 75004 Paris",output="latlon")
	scible<-SpatialPointsDataFrame(coords=cible, data=cible, coords.nrs = numeric(0), proj4string = CRS(WGS84), match.ID = TRUE)
	id_entite_geo<-over(scible,gadm)$ID_5
	entite_geo<-subset(gadm,ID_5%in%id_entite_geo)

	##OSM LAYERS
	path.osm<-"http://download.geofabrik.de/europe/france/ile-de-france-latest.shp.zip"
	path.out<-"YOUR_PATH"
	download.file(path.osm, "parisosm.zip")
	dir.create(path.out)
	unzip("parisosm.zip")
	unzip("parisosm.zip",exdir=path.out)
	files<-list.files(path.out,pattern=".shp")
	ifiles<-substr(files,1,str_locate(files,".shp")[,1]-1)
	lapply(1:length(ifiles),function(i){assign(toupper(ifiles[i]),readOGR(dsn=path.out,layer=ifiles[i]),inherits=TRUE)})

	##LA GRILLE
	NB_SIDE_BLOCS<-10
	start_point<-spTransform(scible,CRS(LIIE))
	SIDE<-NB_SIDE_BLOCS/2
	coord<-matrix(c(coordinates(start_point)[1],coordinates(start_point)[2]+100SIDE,coordinates(start_point)[1],coordinates(start_point)[2]-100SIDE),nr=2,byrow=TRUE)
	TEMP1<-pave(coord,round(SIDE*1.25,0),SIDE,output="spdf")
	TEMP2<-pave(coord,round(SIDE*1.25,0),SIDE,output="spdf",ydown=FALSE)
	TEMP1<-spChFIDs(TEMP1,as.character(paste("S1",rownames(TEMP1@data),sep="")))
	TEMP2<-spChFIDs(TEMP2,as.character(paste("S2",rownames(TEMP2@data),sep="")))
	BLOCS<-spRbind(TEMP1,TEMP2)
	co<-coordinates(BLOCS)
	BLOCS@data$ID2<-1:length(BLOCS)
	BLOCS<-spChFIDs(BLOCS,as.character(BLOCS@data$ID2))
	proj4string(BLOCS)<-CRS(LIIE)
	BLOCS<-spTransform(BLOCS,CRS(proj4string(LANDUSE)))
	UBLOCS<-unionSpatialPolygons(BLOCS,rep(1,length(BLOCS)))

	##OVER
	oLANDUSE<-over(LANDUSE,UBLOCS)
	oLANDUSE<-LANDUSE[which(!is.na(oLANDUSE)),]
	oNATURAL<-over(NATURAL,UBLOCS)
	oNATURAL<-NATURAL[which(!is.na(oNATURAL)),]
	oPLACES<-over(PLACES,UBLOCS)
	oPLACES<-PLACES[which(!is.na(oPLACES)),]
	oPOINTS<-over(POINTS,UBLOCS)
	oPOINTS<-POINTS[which(!is.na(oPOINTS)),]
	oPOINTS2<-over(POINTS,UBLOCS_BIS)
	oPOINTS2<-POINTS[which(!is.na(oPOINTS2)),]
	oRAILWAYS<-over(RAILWAYS,UBLOCS)
	oRAILWAYS<-RAILWAYS[which(!is.na(oRAILWAYS)),]
	oROADS<-over(ROADS,UBLOCS)
	oROADS<-ROADS[which(!is.na(oROADS))]

	##SUBSET DES ÉLÉMENTS À AFFICHER
	BASIN<-subset(oLANDUSE,type=="basin")
	CAFE<-subset(oPOINTS2,type=="cafe")
	GRASS<-subset(oLANDUSE,type=="grass")
	HOSPITALS<-subset(oPOINTS,type=="hospital")
	INTERNET<-subset(oPOINTS,type=="internet_cafe")
	METRO<-subset(RAILWAYS,type=="subway")
	METRO1<-subset(METRO,name=="Métro 1")
	METRO5<-subset(METRO,name=="Métro 5")
	METRO7<-subset(METRO,name=="Métro 7")
	METRO8<-subset(METRO,name=="Métro 8")
	PARK<-subset(oNATURAL,type=="park")
	PUB<-subset(oPOINTS2,type=="pub")
	RIVERBANK<-subset(oNATURAL,type=="riverbank")
	STATIONS<-subset(oPOINTS,type=="station")
	ST1<-STATIONS[which(STATIONS@data$name=="Bastille"),][2,]
	ST2<-STATIONS[which(STATIONS@data$name %in% c("Pont Marie","Sully - Morland","Saint-Paul")),]
	NSTATIONS<-spRbind(ST1,ST2)
	WATER<-subset(oNATURAL,type=="water")

	#PLOTS DIAPOSITIVES 7-16
	plot(UBLOCS,col="grey",lwd=5,border="grey")
	plot(RIVERBANK,add=TRUE,col="#9bbffe",lwd=1,border="#8db7fd")
	plot(WATER,add=TRUE,col="#9bbffe",lwd=1,border="#8db7fd")
	plot(BASIN,add=TRUE,col="#9bbffe",lwd=1,border="#8db7fd")
	plot(PARK,add=TRUE,col="#cadfaa",lwd=0.5)
	plot(GRASS,add=TRUE,col="#cadfaa",lwd=0.5)
	plot(OROADS,add=T,col="white",border="black")
	# Places (sans Paris et les arrondissements)
	PLACES_clean = subset(PLACES, !grepl("Arrondissement", PLACES@data$name))
	PLACES_clean = subset(PLACES_clean, PLACES_clean@data$name != "Paris")
	coordPLACESdecalees = coordinates(PLACES_clean)
	coordPLACESdecalees[,2] = coordPLACESdecalees[,2]+0.0005
	plot(PLACES_clean,add=TRUE,col="brown",lwd=5,pch=0)
	text(coordPLACESdecalees, offset = 20, labels=PLACES_clean@data$name, cex = 0.7, col = "brown")
	plot(scible,add=TRUE,lwd=2,cex=2.5)
	plot(scible,add=TRUE,pch=21,cex=1.7,lwd=2,bg="orange")
	plot(METRO7,add=TRUE,col="#ff82b4",lwd=5)
	plot(METRO1,add=TRUE,col="#ffcd09",lwd=5)
	plot(METRO8,add=TRUE,col="#cfacd4",lwd=5)
	plot(METRO5,add=TRUE,col="#f28e42",lwd=5)
	plot(NSTATIONS,col="white",cex=1.1,pch=19,add=T)
	plot(NSTATIONS,col="#1a3c90",cex=1.2,pch="O",add=T,lwd=2)
	plot(NSTATIONS,col="#1a3c90",cex=0.5,lwd=5,pch="M",add=T)
	plot(CAFE,col="white",cex=1.1,pch=19,add=T)
	plot(CAFE,col="brown",cex=1.2,pch="O",add=T,lwd=2)
	plot(CAFE,col="brown",cex=0.5,lwd=5,pch="C",add=T)
	plot(PUB,col="white",cex=1.1,pch=19,add=T)
	plot(PUB,col="#075933",cex=1.2,pch="O",add=T,lwd=2)
	plot(PUB,col="#075933",cex=0.5,lwd=5,pch="P",add=T)
	plot(INTERNET,col="white",cex=1.4,pch=19,add=T)
	plot(INTERNET,col="orange",cex=1.5,pch="O",add=T,lwd=2)
	plot(INTERNET,col="orange",cex=0.7,lwd=5,pch="IC",add=T)
	plot(HOSPITALS,col="white",cex=1.6,border="black",pch=15,add=T)
	plot(HOSPITALS,col="grey",cex=1.7,pch=0,add=T,lwd=2)
	plot(HOSPITALS,col="red",cex=0.8,lwd=5,pch="H",add=T)

	##PLOT DIAPOSITIVE 17
	plot(UBLOCS)
	plot(oROADS,add=TRUE,lwd=2.5)
	plot(scible,add=TRUE,lwd=2,cex=5)
	plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

	##PLOT DIAPOSITIVE 18
	plot(UBLOCS)
	plot(oROADS,add=TRUE,lwd=2.5,col="gray")
	plot(BLOCS,add=TRUE,border="firebrick",lwd=6)
	plot(BLOCS,add=TRUE,border="yellow",lwd=1.5)
	plot(scible,add=TRUE,lwd=2,cex=5)
	plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

	##SUBSET 4 BLOCS
	sblocs<-subset(BLOCS,ID%in%c(3,4))
	usblocs<-unionSpatialPolygons(sblocs,rep(1,length(sblocs)))

	##ROADS dans 4 BLOCS
	ov<-over(ROADS,usblocs)
	w<-which(!is.na(ov))
	oroads<-ROADS[w,]

	##POINTS dans 4 BLOCS
	ov<-over(POINTS,usblocs)
	w<-which(!is.na(ov))
	opoints<-POINTS[w,]
	copoints<-data.frame(coordinates(opoints)) #coordinates pour la represenation graphique
	colnames(copoints)<-c("lon","lat")

	##PLOT DIAPOSITIVE 19
	plot(usblocs)
	plot(sblocs,add=TRUE,border="firebrick",lwd=10)
	plot(sblocs,add=TRUE,border="yellow",lwd=3)
	plot(oroads,add=TRUE,lwd=2.5)
	plot(opoints,add=TRUE,cex=1.8,lwd=3)
	plot(opoints,add=TRUE,pch=21,cex=1.8,lwd=3,bg="lightgreen")
	plot(scible,add=TRUE,lwd=2,cex=3)
	plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

	##FONCTIONS AUXILIAIRES
	source("YOUR_PATH/functions_raddicts.R")

	uu<-1.2
	path<-"YOUR_PATH/RADDICTS.png"
	PNG.to.rasterGrob(Width=uu*4,Height=uu,Unit="cm",path)

	##GGPLOTS DIAPOSITIVE 20,21,22 Il suffit de changer GMAP par GMAP1,GMAP2,GMAP3 dans la ligne 110
	GMAP1 <- ggmap(get_map(location=bbox(usblocs),zoom=18, source="google",maptype = "roadmap"), extent = "panel", maprange=FALSE)
	GMAP2 <- ggmap(get_map(location=bbox(usblocs),zoom=18, source="stamen",maptype = "toner"), extent = "panel", maprange=FALSE)
	GMAP3 <- ggmap(get_map(location=bbox(usblocs),zoom=18, source="google",maptype = "hybrid"), extent = "panel", maprange=FALSE)

	GMAP +stat_density2d(data = copoints, aes(x = lon, y = lat, fill = ..level..,col=NULL, alpha = ..level..),size = 0.01,h=0.0008, bins = 45, geom = 'polygon')+
	geom_density2d(data = copoints, aes(x = lon, y = lat,col="darkorange4"),h=0.0008)+
	scale_fill_gradient(low = 'yellow', high = 'firebrick')+
	scale_alpha(range = c(0.00, 0.75), guide = FALSE)+
	annotate('custom', x=cible$lon,y=cible$lat, grob=gRADDICTS)+
	theme(axis.text.x=element_blank(),axis.ticks=element_blank(),
	axis.title.x = element_blank(),axis.text.y=element_blank(),axis.title.y = element_blank())+
	theme(legend.position='bottom',legend.title = element_text(colour='black', size=14, face='bold'))+
	theme(legend.text = element_text(colour='black', size = 7, face = 'bold'))+
	theme(plot.title = element_text(size=12,lineheight=.8, face='bold'))+
	theme(plot.background = element_rect(fill='white ', colour='white '))+
	theme(panel.background=element_rect(fill='white ', colour='white '),
	panel.grid.major = element_blank(),panel.grid.minor = element_blank())+
	theme(panel.border=element_rect(colour = 'black',fill='transparent', size=6))+
	theme(legend.position = 'none', axis.title = element_blank(), text = element_text(size = 12))


	##PLOTKML DIAPOSITIVE 23
	kml(entite_geo1, ".", "zone_rdv.kml", colour = "yellow", altitude = 100, alpha = 0.5)
	kml(scible, ".", "point_rdv.kml", colour = "red", width = 2, labels="R addicts Paris", LabelScale = 1.2, size = 1.2)

	##RMAPS PAGES 24-26
	opoints<-subset(OPOINTS,!is.na(name))
	copo<-data.frame(coordinates(opoints))
	colnames(copo)<-c("lon","lat")
	copo$popup<-opoints@data$name
	map = Leaflet$new()
	map$setView(c(cible$lat, cible$lon ), 16)
	PROVIDER<-"Stamen.Toner" #Autres Stamen.Watercolor (p.25), Esri.WorldImagery (p.26)
	map$tileLayer(provider="Stamen.Watercolor")
	for(i in sample(1:nrow(copo),20)){
	map$marker(c(copo$lat[i],copo$lon[i]), bindPopup = copo$popup[i])
	}


	##STRUCTURES DE VOISINAGE
	nbpoly1<-poly2nb(BLOCS,queen=TRUE)

	##PLOT DIAPOSITIVE 28
	plot(BLOCS, border="grey")
	plot(oROADS,add=TRUE,lwd=4,col="cornflowerblue")
	plot(oROADS,add=TRUE,lwd=2,col="gray")
	plot(nbpoly1, coordinates(BLOCS), add=TRUE,bg="gold",pch=21,lwd=2)
	plot(UBLOCS, lwd=6,add=TRUE)
	plot(scible,add=TRUE,lwd=2,cex=3)
	plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

	##PLOT DIAPOSITIVE 29
	nbpoly2<-knn2nb(knearneigh(coordinates(BLOCS), k = 1), row.names = BLOCS@data$ID2)
	plot(BLOCS, border="grey")
	plot(oROADS,add=TRUE,lwd=4,col="cornflowerblue")
	plot(oROADS,add=TRUE,lwd=2,col="gray")
	plot(nbpoly2, coordinates(BLOCS), add=TRUE,bg="gold",pch=21,lwd=2)
	plot(UBLOCS, lwd=6,add=TRUE)
	plot(scible,add=TRUE,lwd=2,cex=3)
	plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

	##PLOT DIAPOSITIVE 30
	nbpoly2<-knn2nb(knearneigh(coordinates(BLOCS), k = 2), row.names = BLOCS@data$ID2)
	plot(BLOCS, border="grey")
	plot(oROADS,add=TRUE,lwd=4,col="cornflowerblue")
	plot(oROADS,add=TRUE,lwd=2,col="gray")
	plot(nbpoly2, coordinates(BLOCS), add=TRUE,bg="gold",pch=21,lwd=2)
	plot(UBLOCS, lwd=6,add=TRUE)
	plot(scible,add=TRUE,lwd=2,cex=3)
	plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")


	##ENVIRONNEMENT SOURCE DIAPOSITIVE 31
	ov<-over(POINTS,UBLOCS)
	w<-which(!is.na(ov))
	OPOINTS<-POINTS[w,]
	ov<-over(OPOINTS,BLOCS)
	ID<-names(table(ov$ID2))
	DATA<-data.frame(cbind(ID=ID))
	DATA$Y<-table(ov$ID2)
	ov<-over(oROADS,BLOCS)
	ID<-names(table(ov$ID2))
	DATA2<-data.frame(cbind(ID=ID))
	DATA2$X<-table(ov$ID2)
	DATA<-merge(DATA,DATA2,by="ID");remove(DATA2) #Y nb points du bloc, X nb troncons de rue du bloc


	##PLOT DIAPOSITIVE 31
	plot(UBLOCS)
	plot(oROADS,add=TRUE,lwd=2.5,col="gray")
	plot(OPOINTS,add=TRUE,pch=21,bg="lightgreen")
	plot(BLOCS,add=TRUE,lwd=6)
	plot(scible,add=TRUE,lwd=2,cex=5)
	plot(scible,add=TRUE,lwd=2,cex=3)
	plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

	##ENVIRONNEMENT SKATER
	rownames(BLOCS@data)<-BLOCS@data$ID2
	temp<-merge(BLOCS@data,DATA,by.x="ID2",by="ID",all.x=TRUE,all.y=FALSE)
	rownames(temp)<-temp$ID2
	BLOCS@data<-temp[rownames(BLOCS@data),]
	BLOCS@data[is.na(BLOCS@data)]<-0
	nbpoly1<-poly2nb(BLOCS,queen=FALSE)
	lcosts <- nbcosts(nbpoly1, scale(BLOCS@data[,c("Y","X")])) #Definition des poids arbitraires a titre d exemple.
	nb.w <- nb2listw(nbpoly1, lcosts, style="B")
	mst.bh <- mstree(nb.w,5)

	##PLOT DIAPOSITIVE 32
	plot(mst.bh, coordinates(BLOCS), col="steelblue",lwd=6,cex.lab=0.01, cex.circles=0.035, fg="blue")
	plot(BLOCS, border=gray(.5), add=TRUE,lws=3)
	plot(UBLOCS,add=TRUE,lwd=6)
	plot(scible,add=TRUE,lwd=2,cex=3)
	plot(scible,add=TRUE,pch=21,cex=3,lwd=2,bg="orange")

	##FONCTION SKATER
	res1 <- skater(edges=mst.bh[,1:2], data=scale(BLOCS@data[,c("Y","X")]),ncuts=5) #Si critere nb groups, a changer ncuts par 5, 15 et 18
	res1 <- skater(edges=mst.bh[,1:2], data=scale(BLOCS@data[,c("Y","X")]),vec.crit=BLOCS@data$Y,crit=400) #Si critere variable quantitative, a changer crit par 300 et 400
	skBLOCS<-unionSpatialPolygons(BLOCS,res1$groups) #regroupement de mailles geographiques

	##PLOT DIAPOSITIVE 33-37
	plot(skBLOCS,col=brewer.pal(length(unique(res1$groups)), "PuBu")) #autres palettes: Set3, Pastel1, Pastel2, Blues
	plot(oROADS,add=TRUE,lwd=2.5,col="gray")
	plot(BLOCS,add=TRUE,border="gray34")
	plot(UBLOCS,lwd=6,add=TRUE)
	plot(skBLOCS,lwd=6,add=TRUE)
	plot(OPOINTS,add=TRUE,pch=21,cex=0.5,bg="lightgreen")


	##ENVIRONNEMENT MODELES SPATIAUX
	# Modèle linéaire simple (non-spatial)
	BLOCS@data$X<-as.numeric(BLOCS@data$X)
	BLOCS@data$Y<-as.numeric(BLOCS@data$Y)
	lmmod=lm(Y~X,data=BLOCS@data)
	summary(lmmod)
	lm.morantest(lmmod, nb2listw(nbpoly1)) # rejet de H0 : absence d'autocorrélation spatiale dans les résidus
	# conversion matrice voisinage
	W.mat=nb2mat(nbpoly1)
	W.listw=nb2listw(nbpoly1,style="W")
	# modèle LAG par max de vraissemblance
	lagmodel = lagsarlm(Y~X,data=BLOCS@data,listw=W.listw)
	summary(lagmodel) # modèle significatif + rejet d'absence d'autocorrélation spatiale à 0.97% -> bonne prise en compte de l'autocorrélation spatiale
	# modèle durbin par max de vraissemblance
	durbin = lagsarlm(Y~X,data=BLOCS@data,listw=W.listw,type="mixed")
	summary(durbin)
	# modèle MESS
	semmod = lagmess(Y~X,data=BLOCS@data,listw=W.listw)
	summary(semmod)
	# modèle SEM
	semmod = spautolm(Y ~ X,data=BLOCS@data,listw=W.listw)
	summary(semmod)
	# comparaison de modèles
	# méthode 1 : les tests
	lm.LMtests(lmmod,listw=W.listw, test="all")
	# méthode 2 : sélection sur critère
	AIC(lmmod,lagmodel,durbin,semmod)
	# impacts
	impacts(lagmodel, listw=W.listw)
	# erreur relative : LM vs. LAG
	FIT<-predict(lagmodel,BLOCS@data,W.listw)
	FIT<-data.frame(cbind(LAG=fitted(lagmodel),LM=fitted(lmmod),Y=BLOCS@data$Y))
	head(FIT)
	FIT$ELAG<-abs(FIT$Y-FIT$LAG)/FIT$Y
	FIT$ELM<-abs(FIT$Y-FIT$LM)/FIT$Y
	FIT$Y<-NULL
	FIT[FIT>10^10]<-0
	BLOCS@data<-cbind(BLOCS@data,FIT)
	BLOCS@data$ERREUR_LAG<-BLOCS@data$ELAG
	BLOCS@data$ERREUR_LM<-BLOCS@data$ELM

	#SPPLOT DIAPOSITIVE 38
	spplot(BLOCS,c("ERREUR_LAG","ERREUR_LM"),col.regions = brewer.pal(9, "PuBu"))