#script to compare to the "Alertas Tempranas" system #data from "Alertas Tempranas" are in a folder called /raster_alertas #takes file fires.Rdata from predict.R and compares those results to results from #the national deforestation alert system library(sp) library(spdep) library(INLA) library(raster) library(ggplot2) library(SpatialEpi) library(maptools) library(spatialEco) library(ROCR) library(OptimalCutpoints) library(rgdal) library(viridis) library(plyr) rm(list = ls()) load("fires.Rdata") a13i<-raster("raster_alertas/alert2013i.tif") a13ii<-raster("raster_alertas/alert2013ii.tif") a14i<-raster("raster_alertas/alert2014i.tif") a14ii<-raster("raster_alertas/alert2014ii.tif") sr <- "+proj=longlat +ellps=clrk66 +no_defs" a14iipj <- projectRaster(a14ii, crs = sr) a14ii <- resample(a14iipj, a14i, "bilinear") alert2014 <- cover(a14i, a14ii) a13ipj <- projectRaster(a13i, crs = sr) a13iipj <- projectRaster(a13ii, crs = sr) alert2013 <- cover(a13ipj, a13iipj) def13<-subset(defo, year=="y13") coordinates(def13) = ~ x + y def14<-subset(defo, year=="y14") coordinates(def14) = ~ x + y ras <- raster(ncol = length(unique(coordinates(def13)[,1])), nrow = length(unique(coordinates(def13)[,2]))) extent(ras) <- extent(def13) sra13<-rasterize(def13, ras, "defo", fun = max) ras <- raster(ncol = length(unique(coordinates(def14)[,1])), nrow = length(unique(coordinates(def14)[,2]))) extent(ras) <- extent(def14) sra14<-rasterize(def14, ras, "defo", fun = max) aras1<-resample(alert2013, sra13, method="bilinear") aras2<-resample(alert2014, sra14, method="bilinear") aras1[aras1>0]=1 aras2[aras2>0]=1 aras1[is.na(aras1)]=0 aras2[is.na(aras2)]=0 s1<-as.data.frame(stack(sra13, aras1)) s2<-as.data.frame(stack(sra14, aras2)) s1<-na.omit(s1) s2<-na.omit(s2) colnames(s1)<-c("def","ale") colnames(s2)<-c("def","ale") aler1.rocr<- prediction(s1$ale,s1$def) aler2.rocr<- prediction(s2$ale,s2$def) aler1.auc<-performance(aler1.rocr, measure = "auc") aler2.auc<-performance(aler2.rocr, measure = "auc") aler1.perf<-performance(aler1.rocr, measure = "tpr", x.measure = "fpr") aler2.perf<-performance(aler2.rocr, measure = "tpr", x.measure = "fpr") oc.ale13 <- optimal.cutpoints(X = ale ~ def, tag.healthy = 0, methods = "PrevalenceMatching", data = s1, pop.prev = NULL, control = control.cutpoints(), ci.fit = FALSE, conf.level = 0.95, trace = FALSE) oc.ale14 <- optimal.cutpoints(X = ale ~ def, tag.healthy = 0, methods = "PrevalenceMatching", data = s2, pop.prev = NULL, control = control.cutpoints(), ci.fit = FALSE, conf.level = 0.95, trace = FALSE) sink("alert_auc.txt") print(aler1.auc@y.values) print(aler2.auc@y.values) print(oc.ale13) print(oc.ale14) sink() pdf("alerta2013.pdf") plot(aras1, col=rev(inferno(256))) dev.off() tiff(file = "figure4_alerta2013.tiff", width = 2400, height = 2400, units = "px", res = 600) plot(aras1, col=rev(inferno(256))) dev.off() pdf("alerta2014.pdf") plot(aras2, col=rev(inferno(256))) dev.off() tiff(file = "figure4_alerta2014.tiff", width = 2400, height = 2400, units = "px", res = 600) plot(aras2, col=rev(inferno(256))) dev.off() pdf("auc_ale2013.pdf") plot(aler1.perf) dev.off() tiff(file = "figure4_auc_ale2013.tiff", width = 2400, height = 2400, units = "px", res = 600) plot(aler1.perf) dev.off() pdf("auc_ale2014.pdf") plot(aler2.perf) dev.off() tiff(file = "figure4_auc_ale2014.tiff", width = 2400, height = 2400, units = "px", res = 600) plot(aler2.perf) dev.off() save.image("fires.Rdata")