Install packages

#install.packages('sp')
#install.packages('raster')
#install.packages('rpart')
#install.packages('DT')

Load packages

library(sp)
library(raster)
library(rpart)
library(DT)

Indicate input folder

InDir <- 'C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/'

Load all present-day variables

RasterPlotGMTED_SAP_BINAIR / Higher than 1 km now

GMTED_SAP <- read.table(paste(InDir,"gmted_sap01.asc", sep = ""), 
                            skip = 6,quote="\"", comment.char="")



RasterPlotGMTED_SAP <- raster(paste(InDir, "gmted_sap01.asc", sep = ""))


# Set extent
bb <- extent(90, 150, -10, 20)
extent(RasterPlotGMTED_SAP) <- bb
RasterPlotGMTED_SAP <- setExtent(RasterPlotGMTED_SAP,bb,keepres=FALSE)

# To binary values to create map with all the places which are higher than 1km

RasterPlotGMTED_SAP_BINAIR <- RasterPlotGMTED_SAP 
RasterPlotGMTED_SAP_BINAIR[RasterPlotGMTED_SAP_BINAIR <= 1000] <- 0
RasterPlotGMTED_SAP_BINAIR[RasterPlotGMTED_SAP_BINAIR >= 1000] <- 1

RasterPlotLandNow/ Land Now

LandNow <- read.table(paste(InDir,"landnu.asc", sep = ""),
                            skip = 6,quote="\"", comment.char="")

RasterPlotLandNow <- raster(paste(InDir, "landnu.asc", sep = ""))

# Set extent
bb <- extent(90, 150, -10, 20)
extent(RasterPlotLandNow) <- bb
RasterPlotLandNow <- setExtent(RasterPlotLandNow,bb,keepres=TRUE)

#RasterPlotLandNow
#summary(RasterPlotLandNow)
#head(RasterPlotLandNow)
#View(RasterPlotLandNow)


RasterPlotLandNow[RasterPlotLandNow >= 1] <- 1
RasterPlotLandNow[RasterPlotLandNow < 1] <- 0

RasterPlotLandNow
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : landnu 
## values     : 0, 1  (min, max)
summary(RasterPlotLandNow)
##         landnu
## Min.         0
## 1st Qu.      0
## Median       0
## 3rd Qu.      0
## Max.         1
## NA's         0

RPEucDistTrans / DistTrans(Euclidean distance)

EucDistTransform <- read.table(paste(InDir,"transbuffer.asc", sep = ""),
                            skip = 6,quote="\"", comment.char="")

RPEucDistTrans <- raster(paste(InDir, "transbuffer.asc", sep = ""))

# Set extent
bb <- extent(90, 150, -10, 20)
extent(RPEucDistTrans) <- bb
RPEucDistTrans <- setExtent(RPEucDistTrans,bb,keepres=TRUE)

RPEucDistTrench / DistTrench(Euclidean distance)

EucDistTrench <- read.table(paste(InDir,"trenchbuffer.asc", sep = ""),
                            skip = 6,quote="\"", comment.char="")

RPEucDistTrench <- raster(paste(InDir, "trenchbuffer.asc", sep = ""))

# Set extent
bb <- extent(90, 150, -10, 20)
extent(RPEucDistTrench) <- bb
RPEucDistTrench <- setExtent(RPEucDistTrench,bb,keepres=TRUE)

RPEucDistTrench
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/trenchbuffer.asc 
## names      : trenchbuffer
summary(RPEucDistTrench)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         trenchbuffer
## Min.             0.0
## 1st Qu.     125369.2
## Median      313902.4
## 3rd Qu.     625254.1
## Max.       1468607.0
## NA's             0.0

PEucDistFLTB/ `DistFault(Euclidean distance)

EucDistFLTB <- read.table(paste(InDir,"fltbuffer.asc", sep = ""),
                            skip = 6,quote="\"", comment.char="")

RPEucDistFLTB <- raster(paste(InDir, "fltbuffer.asc", sep = ""))

# Set extent
bb <- extent(90, 150, -10, 20)
extent(RPEucDistFLTB) <- bb
RPEucDistFLTB <- setExtent(RPEucDistFLTB,bb,keepres=TRUE)

RPEucDistFLTB
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/fltbuffer.asc 
## names      : fltbuffer
summary(RPEucDistFLTB)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         fltbuffer
## Min.          0.0
## 1st Qu.  122122.0
## Median   307665.0
## 3rd Qu.  593436.4
## Max.    2683221.0
## NA's          0.0

RPEucDistPlateBoundB/ DistPlateBound(Euclidean distance)

EucDistPlateBoundB <- read.table(paste(InDir,"plateboundbuffer.asc", sep = ""),
                            skip = 6,quote="\"", comment.char="")

RPEucDistPlateBoundB <- raster(paste(InDir, "plateboundbuffer.asc", sep = ""))

# Set extent
bb <- extent(90, 150, -10, 20)
extent(RPEucDistPlateBoundB) <- bb
RPEucDistPlateBoundB <- setExtent(RPEucDistPlateBoundB,bb,keepres=TRUE)

RPEucDistPlateBoundB
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/plateboundbuffer.asc 
## names      : plateboundbuffer
summary(RPEucDistPlateBoundB)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         plateboundbuffer
## Min.                 0.0
## 1st Qu.          69915.4
## Median          193328.4
## 3rd Qu.         451116.8
## Max.           1294251.0
## NA's                 0.0

RPEucDistVolcB/ DistVolc(Euclidean distance)

EucDistVolcB <- read.table(paste(InDir,"volcbuffer.asc", sep = ""),
                            skip = 6,quote="\"", comment.char="")

RPEucDistVolcB <- raster(paste(InDir, "volcbuffer.asc", sep = ""))

# Set extent
bb <- extent(90, 150, -10, 20)
extent(RPEucDistVolcB) <- bb
RPEucDistVolcB <- setExtent(RPEucDistVolcB,bb,keepres=TRUE)

RPEucDistVolcB
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/volcbuffer.asc 
## names      : volcbuffer
summary(RPEucDistVolcB)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         volcbuffer
## Min.           0.0
## 1st Qu.   180219.7
## Median    375689.7
## 3rd Qu.   590064.3
## Max.     1494803.0
## NA's           0.0

RasterPlotLITHOLOGY

Lithology <- read.table(paste(InDir,"lithology_sap01.asc", sep = ""),
                                skip = 6,quote="\"", comment.char="")

RasterPlotLITHOLOGY <- raster(paste(InDir, "lithology_sap01.asc", sep = ""))


# Set extent
extent(RasterPlotLITHOLOGY) <- bb
RasterPlotLITHOLOGY <- setExtent(RasterPlotLITHOLOGY,bb,keepres=TRUE)
RasterPlotLITHOLOGY
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/lithology_sap01.asc 
## names      : lithology_sap01 
## values     : -2147483648, 2147483647  (min, max)
RasterPlotLITHOLOGY[RasterPlotLITHOLOGY <= 1] <- 1
RasterPlotLITHOLOGY[RasterPlotLITHOLOGY >= 14] <- 14
RasterPlotLITHOLOGY
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : lithology_sap01 
## values     : 1, 14  (min, max)

RasterPlotPRECIP

Precip  <- read.table(paste(InDir,"annualprecipitation_sap01_g.asc", sep = ""),
                             skip = 6,quote="\"", comment.char="")


RasterPlotPRECIP  <- raster(paste(InDir, "annualprecipitation_sap01_g.asc", sep = ""))


# Set extent
extent(RasterPlotPRECIP ) <- bb
RasterPlotPRECIP  <- setExtent(RasterPlotPRECIP ,bb,keepres=FALSE)
RasterPlotPRECIP[RasterPlotPRECIP  <= 749] <- 749
RasterPlotPRECIP[RasterPlotPRECIP  >= 6302] <- 6302
RasterPlotPRECIP 
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : annualprecipitation_sap01_g 
## values     : 749, 6302  (min, max)

RasterPlotTEMP

Temp_sap <- read.table(paste(InDir,"annualmeantemp_sap01_g.asc", sep = ""),
                                             skip = 6,quote="\"", comment.char="")
RasterPlotTEMP <- raster(paste(InDir, "annualmeantemp_sap01_g.asc", sep = ""))

# Set extent
bb <- extent(90, 150, -10, 20)
extent(RasterPlotTEMP) <- bb
RasterPlotTEMP<- setExtent(RasterPlotTEMP,bb,keepres=TRUE)
RasterPlotTEMP
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/annualmeantemp_sap01_g.asc 
## names      : annualmeantemp_sap01_g

RasterPlotGEOAGE_Mean

geologicage <- read.table(paste(InDir,"geoage_sap01.asc", sep = ""),
                                  skip = 6,quote="\"", comment.char="")

RasterPlotGEOAGE <- raster(paste(InDir, "geoage_sap01.asc", sep = ""))

# Set extent
extent(RasterPlotGEOAGE) <- bb
RasterPlotGEOAGE <- setExtent(RasterPlotGEOAGE,bb,keepres=TRUE)
#RasterPlotGEOAGE


RasterPlotGEOAGE_Mean <- RasterPlotGEOAGE
#Mean age from geologic era 
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 1] <- NA
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 2] <- NA
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 3] <- 226.6
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 4] <- 250.1
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 5] <- 400.0
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 6] <- 173.2
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 7] <- 1.3
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 8] <- 359.1
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 9] <- 450.0
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 10] <- 492.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 11] <- 198.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 12] <- 328.9
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 13] <- 159.0
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 14] <- 389.1
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 15] <- 438.8
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 16] <- 600
#RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 16] <- 1520.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 17] <- 305.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 18] <- 306.0
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 19] <- 11.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 20] <- NA
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 21] <- 12.8
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 22] <- 133.7
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 23] <- 513.2
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 24] <- 33.0
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 25] <- 105.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 26] <- 303.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 27] <- 34.3
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 28] <- 600
#RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 28] <- 1492.7
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 29] <- 275.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 30] <- 600
#RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 30] <- 2570.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 31] <- 72.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 32] <- 251.9
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 33] <- 44.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 34] <- 159.0
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 35] <- 251.9
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 36] <- 396.5
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 37] <- 600
#RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 37] <- 3250.0
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 38] <- 280.1
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 39] <- 464.6
RasterPlotGEOAGE_Mean[RasterPlotGEOAGE_Mean == 40] <- 431.5
RasterPlotGEOAGE_Mean
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : geoage_sap01 
## values     : 1.3, 600  (min, max)

RasterPlotPLATE3TYPE = PlateBoundary now

Platebound3type <- read.table(paste(InDir,"plateboundtype3_sap01.asc", sep = ""),
                                      skip = 6,quote="\"", comment.char="")

RasterPlotPLATE3TYPE <- raster(paste(InDir, "plateboundtype3_sap01.asc", sep = ""))

# Set extent
extent(RasterPlotPLATE3TYPE) <- bb
RasterPlotPLATE3TYPE <- setExtent(RasterPlotPLATE3TYPE,bb,keepres=TRUE)
RasterPlotPLATE3TYPE
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/plateboundtype3_sap01.asc 
## names      : plateboundtype3_sap01 
## values     : -2147483648, 2147483647  (min, max)
RasterPlotPLATE3TYPE[RasterPlotPLATE3TYPE >= 1] <- 1
RasterPlotPLATE3TYPE
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : plateboundtype3_sap01 
## values     : 1, 1  (min, max)
# Make buffer
RasterPlotPLATE3TYPE <- buffer(RasterPlotPLATE3TYPE, width = 20000)
RasterPlotPLATE3TYPE
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 1  (min, max)
RasterPlotPLATE3TYPE[is.na(RasterPlotPLATE3TYPE)] <- 0
RasterPlotPLATE3TYPE
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)

RasterPlotVOLC_Buffer

Volcanoes <- read.table(paste(InDir,"volcanoes_sap01.asc", sep = ""),
                                skip = 6,quote="\"", comment.char="")

RasterPlotVOLC <- raster(paste(InDir, "volcanoes_sap01.asc", sep = ""))

# Set extent

extent(RasterPlotVOLC) <- bb
RasterPlotVOLC <- setExtent(RasterPlotVOLC,bb,keepres=TRUE)
RasterPlotVOLC
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/volcanoes_sap01.asc 
## names      : volcanoes_sap01 
## values     : -2147483648, 2147483647  (min, max)
RasterPlotVOLC[RasterPlotVOLC >= 1] <- 1
RasterPlotVOLC
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : volcanoes_sap01 
## values     : 1, 1  (min, max)
RasterPlotVOLC_Buffer <- buffer(RasterPlotVOLC, width = 20000)
RasterPlotVOLC_Buffer
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 1  (min, max)
RasterPlotVOLC_Buffer[is.na(RasterPlotVOLC_Buffer)] <- 0
RasterPlotVOLC_Buffer
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)

RasterPlotTRANS_Buffer `

Transform <- read.table(paste(InDir,"transform_sap01.asc", sep = ""),
                                skip = 6,quote="\"", comment.char="")

RasterPlotTRANS <- raster(paste(InDir, "transform_sap01.asc", sep = ""))

# Set extent

extent(RasterPlotTRANS) <- bb
RasterPlotTRANS <- setExtent(RasterPlotTRANS,bb,keepres=TRUE)
RasterPlotTRANS
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/transform_sap01.asc 
## names      : transform_sap01 
## values     : -2147483648, 2147483647  (min, max)
RasterPlotTRANS[RasterPlotTRANS >= 1] <- 1
RasterPlotTRANS
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : transform_sap01 
## values     : 1, 1  (min, max)
# Make buffer 

RasterPlotTRANS_Buffer <- buffer(RasterPlotTRANS, width = 20000)
RasterPlotTRANS_Buffer
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 1  (min, max)
RasterPlotTRANS_Buffer[is.na(RasterPlotTRANS_Buffer)] <- 0
RasterPlotTRANS_Buffer
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)

RasterPlotTRENCH_Buffer

Trench <- read.table(paste(InDir,"trench_sap01.asc", sep = ""),
                             skip = 6,quote="\"", comment.char="")

RasterPlotTRENCH <- raster(paste(InDir, "trench_sap01.asc", sep = ""))                                  
# Set extent

extent(RasterPlotTRENCH) <- bb
RasterPlotTRENCH <- setExtent(RasterPlotTRENCH,bb,keepres=TRUE)
RasterPlotTRENCH
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/trench_sap01.asc 
## names      : trench_sap01 
## values     : -2147483648, 2147483647  (min, max)
RasterPlotTRENCH[RasterPlotTRENCH >= 1] <- 1
RasterPlotTRENCH
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : trench_sap01 
## values     : 1, 1  (min, max)
# Make buffer 

RasterPlotTRENCH_Buffer <- buffer(RasterPlotTRENCH, width = 20000)
RasterPlotTRENCH_Buffer
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 1  (min, max)
RasterPlotTRENCH_Buffer[is.na(RasterPlotTRENCH_Buffer)] <- 0
RasterPlotTRENCH_Buffer
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)

RasterPlotFAULTLINETYPE_Buffer

FaultLinesType <- read.table(paste(InDir,"faultlinestype_sap01.asc", sep = ""),
                                     skip = 6,quote="\"", comment.char="")

RasterPlotFAULTLINETYPE <- raster(paste(InDir, "faultlinestype_sap01.asc", sep = ""))                                 
# Set extent

extent(RasterPlotFAULTLINETYPE) <- bb
RasterPlotFAULTLINETYPE <- setExtent(RasterPlotFAULTLINETYPE,bb,keepres=TRUE)
RasterPlotFAULTLINETYPE
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/faultlinestype_sap01.asc 
## names      : faultlinestype_sap01 
## values     : -2147483648, 2147483647  (min, max)
RasterPlotFAULTLINETYPE[RasterPlotFAULTLINETYPE >= 1] <- 1
RasterPlotFAULTLINETYPE
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : faultlinestype_sap01 
## values     : 1, 1  (min, max)
# Make buffer 

RasterPlotFAULTLINETYPE_Buffer <- buffer(RasterPlotFAULTLINETYPE, width = 20000)
RasterPlotFAULTLINETYPE_Buffer
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1, 1  (min, max)
RasterPlotFAULTLINETYPE_Buffer[is.na(RasterPlotFAULTLINETYPE_Buffer)] <- 0
RasterPlotFAULTLINETYPE_Buffer
## class      : RasterLayer 
## dimensions : 300, 600, 180000  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0, 1  (min, max)

Plot all current variables

plot(RasterPlotLandNow, main = 'Land Now')

plot(RasterPlotGMTED_SAP_BINAIR , main = 'Higher than 1km now' )

plot( RasterPlotGEOAGE_Mean, main = 'GeoAge Mean Now' )

plot(RasterPlotLITHOLOGY , main = 'Lithology now' )

plot( RasterPlotTEMP, main = 'Temp now' )

plot(RasterPlotPRECIP , main = 'precip now' )

plot(RPEucDistTrans , main = 'Dist trans now' )

plot(RPEucDistFLTB , main = 'dist faults now' )

plot( RasterPlotFAULTLINETYPE_Buffer, main = 'faults now' )

plot( RPEucDistPlateBoundB, main = 'dist platebound now' )

plot(RPEucDistTrench , main = 'dist trench now' )

plot( RPEucDistVolcB, main = 'dist volc now' )

plot(RasterPlotVOLC_Buffer , main = 'volc now' )

plot(RasterPlotTRANS_Buffer , main = 'trans now' )

plot(RasterPlotTRENCH_Buffer , main = 'trench now' )

plot(RasterPlotPLATE3TYPE , main = 'platebounday now' )

Make brick of all variables

EurekaAll_Stack <- stack(RasterPlotLandNow, RasterPlotGMTED_SAP_BINAIR, RasterPlotGEOAGE_Mean, RasterPlotLITHOLOGY, RasterPlotTEMP, RPEucDistTrans, RPEucDistFLTB, RasterPlotFAULTLINETYPE_Buffer, RPEucDistPlateBoundB, RasterPlotPRECIP, RPEucDistTrench, RPEucDistVolcB, RasterPlotVOLC_Buffer, RasterPlotTRANS_Buffer, RasterPlotTRENCH_Buffer, RasterPlotPLATE3TYPE) 

EurekaAll_Brick <- brick(EurekaAll_Stack)
EurekaAll_Brick
## class      : RasterBrick 
## dimensions : 300, 600, 180000, 16  (nrow, ncol, ncell, nlayers)
## resolution : 0.1, 0.1  (x, y)
## extent     : 90, 150, -10, 20  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Linde/AppData/Local/Temp/RtmpymZ2Ge/raster/r_tmp_2021-06-30_151009_3228_31529.grd 
## names      :       landnu,  gmted_sap01, geoage_sap01, lithology_sap01, annualmeantemp_sap01_g,  transbuffer,    fltbuffer,      layer.1, plateboundbuffer, annualpre//on_sap01_g, trenchbuffer,   volcbuffer,      layer.2,      layer.3,      layer.4, ... 
## min values :     0.000000,     0.000000,     1.300000,        1.000000,               4.108334,     0.000000,     0.000000,     0.000000,         0.000000,            749.000000,     0.000000,     0.000000,     0.000000,     0.000000,     0.000000, ... 
## max values : 1.000000e+00, 1.000000e+00, 6.000000e+02,    1.400000e+01,           2.864583e+01, 1.506163e+06, 2.687642e+06, 1.000000e+00,     1.294251e+06,          6.302000e+03, 1.472936e+06, 1.506296e+06, 1.000000e+00, 1.000000e+00, 1.000000e+00, ...
EurekaAll_Brick_df <- as.data.frame(EurekaAll_Brick)

Change the names of each variable to make sense

names(EurekaAll_Brick_df)
##  [1] "landnu"                      "gmted_sap01"                
##  [3] "geoage_sap01"                "lithology_sap01"            
##  [5] "annualmeantemp_sap01_g"      "transbuffer"                
##  [7] "fltbuffer"                   "layer.1"                    
##  [9] "plateboundbuffer"            "annualprecipitation_sap01_g"
## [11] "trenchbuffer"                "volcbuffer"                 
## [13] "layer.2"                     "layer.3"                    
## [15] "layer.4"                     "layer.5"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "landnu"] <- "Land"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "gmted_sap01"] <- "HigherThan1km"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "geoage_sap01"] <- "GeoMean"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "lithology_sap01"] <- "Lithology"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "annualmeantemp_sap01_g"] <- "Temp"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "annualprecipitation_sap01_g"] <- "Precip"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "transbuffer"] <- "DistTrans"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "fltbuffer"] <- "DistFault"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "layer.1"] <- "Faults"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "plateboundbuffer"] <- "DistPlateBound"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "trenchbuffer"] <- "DistTrench"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "volcbuffer"] <- "DistVolc"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "layer.2"] <- "Volc"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "layer.3"] <- "Trans"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "layer.4"] <- "Trench"
names(EurekaAll_Brick_df)[names(EurekaAll_Brick_df) == "layer.5"] <- "PlateBound"

names(EurekaAll_Brick_df)
##  [1] "Land"           "HigherThan1km"  "GeoMean"        "Lithology"     
##  [5] "Temp"           "DistTrans"      "DistFault"      "Faults"        
##  [9] "DistPlateBound" "Precip"         "DistTrench"     "DistVolc"      
## [13] "Volc"           "Trans"          "Trench"         "PlateBound"

See head of dataframe

datatable(head(EurekaAll_Brick_df,10),autoHideNavigation = TRUE)
## Warning in datatable(head(EurekaAll_Brick_df, 10), autoHideNavigation = TRUE):
## `autoHideNavigation` will be ignored if the `pageLength` option is not provided.