Change detection with time series and machine learning

Change detection.R

#' @author Yonatan Tarazona Coronel
#'
# DETECCION DE CAMBIOS CON SERIES TEMPORALES Y MACHINE LEARNING - TELEDETECCION
## 1. APRENDIZAJE SUPERVISADO
### PASOS

#### i. Imágenes y librerías en R
#### ii. Datos de entrenamiento
#### iii. Training-Testing y aprendizaje del modelo
#### iv. Evaluamos el modelo y Aplicar el modelo a individuos nuevos

# Ubicamos nuestro directorio y leemos los paquetes
setwd("G:/Webinar_TimeSeries_Machine Learning R-QGIS/Images")
suppressMessages(library(raster))
suppressMessages(library(rgdal))
suppressMessages(library(kknn))

# Leemos las imágenes a usar
img2011<-brick('L5003069_20110903.tif')
img2016<-brick('L8003069_20160916.tif')
dim(img2011);dim(img2016)

## le asigamos nombres óptimos
# Imagen 2011
names(img2011) <- c(paste("B", 1:6, sep = ""))
names(img2011)
# Imagen 2016
names(img2016) <- c(paste("B", 1:6, sep = ""))
names(img2016)

## Cargamos el vector de muestras (testing y training)
points <- readOGR("G:/Webinar_TimeSeries_Machine Learning R-QGIS/Vector","Points_2011_value")
tabla <- points@data
names<-c("B1","B2","B3","B4","B5","B6","Class")
head(tabla)

# Extraemos los valores de las bandas para cada clase
vegt<-extract(img2011,points)
tabla<-data.frame(vegt,Class=tabla[,1])
colnames(tabla)<-names
head(tabla);dim(tabla)

# Separamos los puntos de entrenamiento en testing y training
Muestra<-sample(1:dim(tabla)[1], 100)
testing<-tabla[Muestra,]
training<-tabla[-Muestra,]
dim(testing);dim(training)

# Aplicamos el modelo para training
model<-train.kknn(Class~.,data=training, kmax = 15)

# Predicción de individuos nuevos y ver la precisión espacial de la clasificación
prediccion<- predict(model, testing[,-7])
MC<-table(testing[,7],prediccion)
kappa<-sum(diag(MC))/sum(MC)
kappa; MC

beginCluster()
Output2011<- clusterR(img2011, raster::predict, args = list(model = model))
endCluster()

plot(Output2011)
writeRaster(Output2011, "Clasificacion-kknn-2011.tif",drivername="GTiff")

## Clasificacion para el 2016
points <- readOGR("G:/Webinar_TimeSeries_Machine Learning R-QGIS/Vector","Points_2016_value")
tabla <- points@data
names<-c("B1","B2","B3","B4","B5","B6","Class")
head(tabla)

vegt<-extract(img2016,points)
tabla<-data.frame(vegt,Class=tabla[,1])
colnames(tabla)<-names
head(tabla);dim(tabla)

Muestra<-sample(1:dim(tabla)[1], 100)
testing<-tabla[Muestra,]
training<-tabla[-Muestra,]
dim(testing);dim(training)

model<-train.kknn(Class~.,data=training, kmax = 15)

prediccion<- predict(model, testing[,-7])
MC<-table(testing[,7],prediccion)
kappa<-sum(diag(MC))/sum(MC)
kappa; MC

beginCluster()
Output2016<- clusterR(img2016, raster::predict, args = list(model = model))
endCluster()

plot(Output2016)
writeRaster(Output2016, "Clasificacion-kknn-2016.tif",drivername="GTiff")

## Cambios entre el 2011-2016

Cambio2016.2011 <- Output2016-Output2011
plot(Cambio2016.2011)
writeRaster(Cambio2016.2011, "Cambios_2016-2011.tif",drivername="GTiff")


## 2. SERIES DE TIEMPO

### 2.1 Directorio y librerias
setwd("G:/Webinar_TimeSeries_Machine Learning R-QGIS/Images/TimeSeriesEVI")
suppressMessages(library(raster))
suppressMessages(library(strucchange))
suppressMessages(library(bfast))
suppressMessages(library(forecast))
suppressMessages(library(rgdal))

## 2.2 Análisis preliminares de series de tiempo

lista.evi<-list.files("G:/Webinar_TimeSeries_Machine Learning R-QGIS/Images/TimeSeriesEVI","_Area.tif")
apilar <- stack(lista.evi)
areglo <- as.array(apilar)
band<-raster(lista.evi[[25]]);plot(band)

## 2.2.1 Descomponer una serie de tiempo
data("CO2")
d.st<-stl(co2, "per")
plot(d.st)

## 2.2.2 Visualización de areas de cambio
### Obtener un raster con valores de filas y columnas
b<-raster(lista.evi[[1]])
for(i in 1:dim(apilar)[1]){
  for(j in 1:dim(apilar)[2]){
    pixel<-paste(i,j,sep = ".")
    b[i,j]<-as.numeric(pixel)
  }
}
writeRaster(b,'Raster_posicion.tif',drivername="GTiff") # Guardar el raster

### Zona de bosque intacto
ts<-ts(areglo[100,110,], start = c(1990),freq = 1)
plot(ts, main="Bosque Tropical", ylim=c(0,1))

### Zona de cambio
ts <- ts(areglo[282,142,], start=c(1990), freq = 1)
plot(ts, type="o")

## 2.2.3 Deteccioin de cambios en Series temporales
### bfast y bfastmonitor
### strucchange
### CCDC

## Analizando cambio en el bosque con "strucchange"
s.t <- brick("MODIS_2000-2016.tif")
areglo.m<- as.array(s.t)
ts <- ts(areglo.m[37,19,], start=c(2000,01), freq = 23)
dd<- breakpoints(ts~1, h = 0.15)
plot(ts, type="l", ylim=c(0.78,0.95))
lines(confint(dd))
lines(fitted(dd),col=4,lwd=3,lty=2)

## Analizando cambio en el bosque con "bfast"
ts <- ts(areglo.m[37,19,], start=c(2000,01), freq = 23)
c.bfast <- bfast(ts,h=0.15, season="harmonic", max.iter=1)
plot(c.bfast)