Import Present-day Variables
Install packages
#install.packages('sp')
#install.packages('raster')
#install.packages('rpart')Load packages
library(sp)
library(raster)
library(rpart)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 variablesRasterPlotGMTED_SAP_BINAIR
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] <- 1RasterPlotLandNow
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 0RPEucDistTrans
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
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 : trenchbuffersummary(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. 125331.4
## Median 311836.3
## 3rd Qu. 620961.8
## Max. 1471339.0
## NA's 0.0RPEucDistFLTB
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 : fltbuffersummary(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. 121632.3
## Median 307917.9
## 3rd Qu. 594894.0
## Max. 2687642.0
## NA's 0.0RPEucDistPlateBoundB
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 : plateboundbuffersummary(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.00
## 1st Qu. 69778.74
## Median 195527.80
## 3rd Qu. 451845.39
## Max. 1290324.00
## NA's 0.00RPEucDistVolcB
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 : volcbuffersummary(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. 181805.6
## Median 375509.4
## 3rd Qu. 590051.9
## Max. 1496972.0
## NA's 0.0RasterPlotLITHOLOGY
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_gRasterPlotGEOAGE_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
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)RasterPlotSEASEDTHICK
SedThick_SAP <- read.table(paste(InDir,"sedimentthick_sap01_g.asc", sep = ""),
skip = 6,quote="\"", comment.char="")
RasterPlotSEASEDTHICK <- raster(paste(InDir, "sedimentthick_sap01_g.asc", sep = ""))
# Set extent
bb <- extent(90, 150, -10, 20)
extent(RasterPlotSEASEDTHICK) <- bb
RasterPlotSEASEDTHICK <- setExtent(RasterPlotSEASEDTHICK,bb,keepres=TRUE)
RasterPlotSEASEDTHICK## 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 : NA
## source : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/sedimentthick_sap01_g.asc
## names : sedimentthick_sap01_g
## values : -2147483648, 2147483647 (min, max)RasterPlotSEASEDTHICK[RasterPlotSEASEDTHICK < 12] <- NA
RasterPlotSEASEDTHICK[RasterPlotSEASEDTHICK > 15307] <- NA
RasterPlotSEASEDTHICK## 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 : NA
## source : memory
## names : sedimentthick_sap01_g
## values : 12, 15307 (min, max)# Plot Raster
plot(RasterPlotSEASEDTHICK, main = "Sediment thickness (m)") 
RasterPlotSEAFLOORAGE
Seafloorage_sap <- read.table(paste(InDir,"seafloorage_sap01_g.asc", sep = ""),
skip = 6,quote="\"", comment.char="")
RasterPlotSEAFLOORAGE <- raster(paste(InDir,"seafloorage_sap01_g.asc", sep = ""))
# Set extent
bb <- extent(90, 150, -10, 20)
extent(RasterPlotSEAFLOORAGE) <- bb
RasterPlotSEAFLOORAGE <- setExtent(RasterPlotSEAFLOORAGE,bb,keepres=TRUE)
RasterPlotSEAFLOORAGE## 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/seafloorage_sap01_g.asc
## names : seafloorage_sap01_g
## values : -2147483648, 2147483647 (min, max)RasterPlotSEAFLOORAGE[RasterPlotSEAFLOORAGE < 30] <- NA
RasterPlotSEAFLOORAGE[RasterPlotSEAFLOORAGE > 260] <- NA
RasterPlotSEAFLOORAGE## 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 : seafloorage_sap01_g
## values : 30, 260 (min, max)# Plot Raster
plot(RasterPlotSEAFLOORAGE, main = "Seafloor age (My)") 
RasterPlotOPHI
Ophiolites <- read.table(paste(InDir,"ophiolites_sap01.asc", sep = ""),
skip = 6,quote="\"", comment.char="")
RasterPlotOPHI <- raster(paste(InDir,"ophiolites_sap01.asc", sep = ""))
# Set extent
extent(RasterPlotOPHI) <- bb
RasterPlotOPHI <- setExtent(RasterPlotOPHI,bb,keepres=TRUE)
RasterPlotOPHI## 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/ophiolites_sap01.asc
## names : ophiolites_sap01
## values : -2147483648, 2147483647 (min, max)RasterPlotOPHI[RasterPlotOPHI >= 1] <- 1
RasterPlotOPHI[is.na(RasterPlotOPHI)] <- 0
RasterPlotOPHI## 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 : ophiolites_sap01
## values : 0, 1 (min, max)plot(RasterPlotOPHI, main = "Ophiolites")
RasterPlotRIDGE_Buffer
Ridges <- read.table(paste(InDir,"ridge_sap01.asc", sep = ""),
skip = 6,quote="\"", comment.char="")
RasterPlotRIDGE <- raster(paste(InDir,"ridge_sap01.asc", sep = ""))
# Set extent
extent(RasterPlotRIDGE) <- bb
RasterPlotRIDGE <- setExtent(RasterPlotRIDGE,bb,keepres=TRUE)
RasterPlotRIDGE## 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/ridge_sap01.asc
## names : ridge_sap01
## values : -2147483648, 2147483647 (min, max)RasterPlotRIDGE[RasterPlotRIDGE >= 1] <- 1
RasterPlotRIDGE## 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 : ridge_sap01
## values : 1, 1 (min, max)# Make buffer
RasterPlotRIDGE_Buffer <- buffer(RasterPlotRIDGE, width = 20000)
RasterPlotRIDGE_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)RasterPlotRIDGE_Buffer[is.na(RasterPlotRIDGE_Buffer)] <- 0
RasterPlotRIDGE_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 Raster
plot(RasterPlotRIDGE_Buffer, main = "Ridge")
********************************** 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' )
plot(RasterPlotSEASEDTHICK, main = 'SeaSed Thick now')
plot(RasterPlotSEAFLOORAGE, main = 'Seafloor Age Now')
plot(RasterPlotOPHI, main = 'Ophiolites Now')
plot(RasterPlotRIDGE_Buffer, main = 'Ridge Now')
Make brick of all variables
EurekaFirst_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, RasterPlotSEASEDTHICK, RasterPlotSEAFLOORAGE, RasterPlotOPHI, RasterPlotRIDGE_Buffer )
EurekaFirst_Brick <- brick(EurekaFirst_Stack)
#EurekaFirst_Brick
EurekaFirst_Brick_df <- as.data.frame(EurekaFirst_Stack)
#EurekaFirst_Brick_dfChange the names of each variable to make sense
names(EurekaFirst_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"
## [17] "sedimentthick_sap01_g" "seafloorage_sap01_g"
## [19] "ophiolites_sap01" "layer.6"names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "landnu"] <- "Land"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "gmted_sap01"] <- "HigherThan1km"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "geoage_sap01"] <- "GeoMean"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "lithology_sap01"] <- "Lithology"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "annualmeantemp_sap01_g"] <- "Temp"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "annualprecipitation_sap01_g"] <- "Precip"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "transbuffer"] <- "DistTrans"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "fltbuffer"] <- "DistFault"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "layer.1"] <- "Faults"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "plateboundbuffer"] <- "DistPlateBound"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "trenchbuffer"] <- "DistTrench"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "volcbuffer"] <- "DistVolc"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "layer.2"] <- "Volc"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "layer.3"] <- "Trans"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "layer.4"] <- "Trench"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "layer.5"] <- "PlateBound"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "sedimentthick_sap01_g"] <- "SeaSedThick"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "seafloorage_sap01_g"] <- "SeaFloorAge"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "ophiolites_sap01"] <- "Ophiolites"
names(EurekaFirst_Brick_df)[names(EurekaFirst_Brick_df) == "layer.6"] <- "Ridge"
names(EurekaFirst_Brick_df)## [1] "Land" "HigherThan1km" "GeoMean" "Lithology"
## [5] "Temp" "DistTrans" "DistFault" "Faults"
## [9] "DistPlateBound" "Precip" "DistTrench" "DistVolc"
## [13] "Volc" "Trans" "Trench" "PlateBound"
## [17] "SeaSedThick" "SeaFloorAge" "Ophiolites" "Ridge"See head of dataframe
head(EurekaFirst_Brick_df)## Land HigherThan1km GeoMean Lithology Temp DistTrans DistFault Faults
## 1 0 0 NA NA NA 606923.0 449970.3 0
## 2 0 0 NA NA NA 596459.9 439506.4 0
## 3 0 0 NA NA NA 585996.6 429042.5 0
## 4 0 0 NA NA NA 575533.4 418578.6 0
## 5 0 0 NA NA NA 565070.0 408114.6 0
## 6 0 0 NA NA NA 554606.6 397650.6 0
## DistPlateBound Precip DistTrench DistVolc Volc Trans Trench PlateBound
## 1 599569.3 NA 293008.6 832105.9 0 0 0 0
## 2 593449.1 NA 282544.2 826564.6 0 0 0 0
## 3 587456.4 NA 272079.8 821124.1 0 0 0 0
## 4 581595.0 NA 261615.4 815786.4 0 0 0 0
## 5 575533.4 NA 251151.0 810553.5 0 0 0 0
## 6 565070.0 NA 240686.5 805427.5 0 0 0 0
## SeaSedThick SeaFloorAge Ophiolites Ridge
## 1 15307 NA 0 0
## 2 15126 NA 0 0
## 3 14757 NA 0 0
## 4 14578 NA 0 0
## 5 14403 NA 0 0
## 6 14224 NA 0 0