#### Data base spaziale
library(maptools)
library(spdep)
setwd("D:\\Lavori\\Lezioni_Analisi_Terr\\PIL")
#IDS<-seq(1,103,by=1)
pilSDF<-readShapePoly("ita_pr00_pil.shp",IDvar="COD_PRO")
COD_PRO<-as(pilSDF,"data.frame")$COD_PRO
summary(pilSDF)
bbox(pilSDF)
names(pilSDF)

### Plot della regione analizzata
plot(pilSDF)
#bbox(pilSDF)
title(main="Province italiane")
invisible(text(getSpPPolygonsLabptSlots(pilSDF), labels=as.character(pilSDF$COD_PRO), cex=0.4))

#### Importazione di matrici spaziali in R


### Matrici di contiguitą ###
prov.w1.nb<-read.gal("ita_pr00.GAL", region.id=COD_PRO)
prov.w1<-nb2listw(prov.w1.nb, style="W",zero.policy=FALSE)
summary(prov.w1)

### Matrici k_Neighbours ###
prov.k3.nb<-read.gwt2nb("ita_pr00_k3.GWT", region.id=COD_PRO)
prov.k3<-nb2listw(prov.k3.nb,style="W",zero.policy=FALSE)
summary(prov.k3)

### Matrici distanza ###
prov.d1.nb<-read.gwt2nb("ita_dis1.GWT",region.id=COD_PRO)
prov.d1<-nb2listw(prov.d1.nb,style="W",zero.policy=FALSE)
summary(prov.d1)


### Moran I ###
MC<-moran.test(as(pilSDF,"data.frame")$NI_PIL_93,prov.w1)
resI <- localmoran(as(pilSDF,"data.frame")$NI_PIL_93,prov.w1,p.adjust.method="none")
fix(resI)

### Scatterplot di Moran ###
moran.plot(as(pilSDF,"data.frame")$NI_PIL_93, prov.w1,labels=as.character(as(pilSDF,"data.frame")$COD_PRO), pch=19)
win.graph()
moran.plot(scale(as(pilSDF,"data.frame")$NI_PIL_93), prov.k3,labels=as.character(as(pilSDF,"data.frame")$NAME), pch=19,xlim=c(-3,3),ylim=c(-3,3),xlab="NI_PIL_93",ylab="W_NI_PIL_93")

x<-as(pilSDF,"data.frame")$NI_PIL_93
zx<-(x-mean(x))/sd(x)

mean(zx)
var(zx)

wfile<-(prov.k3)
wzx<-lag.listw(wfile,zx)

morlm<-lm(wzx ~ zx)
aa<-morlm$coefficients[1]
mori<-morlm$coefficients[2]

aa
mori

par(pty="s") #quadrato
plot(zx,wzx,xlab="NI_PIL_93",ylab="Spatial Lag of NI_PIL_93")
abline(aa,mori)
abline(h=0,lty=2)
abline(v=0,lty=2)
### Moran I ###
str(moran(as(pilSDF,"data.frame")$NI_PIL_93, prov.k3, length(prov.k3.nb), Szero(prov.k3)))

moran.test(as(pilSDF,"data.frame")$NI_PIL_93, prov.k3)# per correlazione spaziale

moran.mc(as(pilSDF,"data.frame")$NI_PIL_93, prov.k3,99)

### Spatially lagged variable ###

xx<-cbind(as(pilSDF,"data.frame")$NI_PIL_93,as(pilSDF,"data.frame")$NI_PIL_95)
Wxx<-data.frame(lag.listw(prov.k3r,xx))
name<-c("W_NI_PIL_93","W_NI_PIL_95")
names(Wxx)<-name

### Moran I multivariata

plot(scale(pilSDF$NI_PIL_95),scale(Wxx$W_NI_PIL_93),xlab="NI_PIL_95",ylab="W_NI_PIL_93",cex=1,xlim=c(-2,2),ylim=c(-2,2))
##identify(pilSDF[,c('pilSDF$NI_PIL_95','Wxx$W_NI_PIL_93')],labels=pilSDF$COD_PRO,cex=0.9)
abline(v=0,lty=3)
abline(h=0,lty=3)
abline(lsfit(scale(pilSDF$NI_PIL_95),scale(Wxx$W_NI_PIL_93)))
ls<-(lsfit(scale(pilSDF$NI_PIL_95),scale(Wxx$W_NI_PIL_93),intercept=F))
#title(main="Moran I spazio-temporale")
options(digits=3)
im<-(ls$coefficients[1])
ims<-toString(im)
title(main=list(paste("Moran I spazio-temporale I=",substr(ims,1,5)),cex=1.7))

summary(lm(pilSDF$NI_PIL_95~Wxx$W_NI_PIL_93))
ls<-lsfit(scale(pilSDF$NI_PIL_95),scale(Wxx$W_NI_PIL_93)) 
### K-nearest neighbor weights con R ###

provknn<-knearneigh(cbind(as(pilSDF,"data.frame")$XCOO,as(pilSDF,"data.frame")$YCOO),k=3)
provknnb<-knn2nb(provknn,row.names=COD_PRO)
prov.k3r<-nb2listw(provknnb,style="W",zero.policy=FALSE)

### Distance based con R ###

provk1<-knearneigh(cbind(as(pilSDF,"data.frame")$XCOO,as(pilSDF,"data.frame")$YCOO))##restituisce una matrice con gli indici delle regioni delle k regioni pił vicine. Default k=1
provknn1<-knn2nb(provk1)##restituisec un oggetto di classe nb
#summary(provknn1)
provdist1<-nbdists(provknn1,cbind(as(pilSDF,"data.frame")$XCOO,as(pilSDF,"data.frame")$YCOO))##restituisce le distanze
provdistvec<-unlist(provdist1)# da lista a vettore
provmaxdist<-max(provdistvec)
provdnb<-dnearneigh(cbind(as(pilSDF,"data.frame")$XCOO,as(pilSDF,"data.frame")$YCOO),0,provmaxdist,row.names=COD_PRO)
#summary(provdnb)
prov.d1r<-nb2listw(provdnb,style="W",zero.policy=FALSE)


brks <- round(quantile(pilSDF@data$NI_PIL_93, probs=seq(0,1,0.2)), digits=2)
colours <- c("gray94", "gray64", "gray44", "gray34", "gray14")
plot(pilSDF, col=colours[findInterval(pilSDF@data$NI_PIL_93, brks, all.inside=TRUE)])
legend(x=c(320000,350000),y=c(3950000,4200000),legend=leglabs(brks),fill=colours, bty="n",cex=0.8)
title(main="Cluster variabili economiche e sociali")
invisible(text(getSpPPolygonsLabptSlots(pilSDF), labels=as.character(pilSDF$COD_PRO), cex=0.4))

localI<-as.data.frame(resI)
max(localI$Ii)
min(localI$Ii)


#brks<-cut(localI$Pr,breaks=(c(0,0.001,0.01,0.05,max(localI$Pr))),labels=c(1,2,3,4))
#brks.n<-as.numeric(brks)
#localI<-cbind(localI,Breaks=brks.n)
#o<-order(brks.n)
#brks.n[o]

brks<-c(0,0.001,0.01,0.05,max(localI$Pr))
colours <- c("DarkGreen", "green3", "GreenYellow", "white")
plot(pilSDF, col=colours[findInterval(localI$Pr, brks, all.inside=TRUE)])
legend(x=c(320000,350000),y=c(3950000,4200000),legend=leglabs(brks),fill=colours, bty="n",cex=0.8)
title(main="LISA significance map (local I)")
invisible(text(getSpPPolygonsLabptSlots(pilSDF), labels=as.character(pilSDF$COD_PRO), cex=0.4))


###### Local G e G*
Gloc<-localG(pilSDF@data$NI_PIL_93, prov.k3, zero.policy=FALSE, spChk=NULL)

Gloc.s<-localG(pilSDF@data$NI_PIL_93, nb2listw(include.self(provknnb),style="W"), zero.policy=FALSE, spChk=NULL)

win.graph()
brks <- round(quantile(Gloc.s, probs=seq(0,1,0.2)), digits=2)
colours <- c("gray94", "gray64", "gray44", "gray34", "gray14")
plot(pilSDF, col=colours[findInterval(Gloc.s, brks, all.inside=TRUE)])
#bbox(pilSDF)
legend(x=c(320000,350000),y=c(3950000,4200000),legend=leglabs(brks),fill=colours, bty="n",cex=0.8)
title(main="Valori Local G*")
invisible(text(getSpPPolygonsLabptSlots(pilSDF), labels=as.character(pilSDF$COD_PRO), cex=0.4))