Install packages

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

Load packages

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

Load Data

load("C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/EurekaAll_Brick_df.RData")

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/'

Import paleo-variables

Land

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

RasterPlotLand35 <- raster(paste(InDir, "Land35.asc", sep = ""))

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

RasterPlotLand35
## 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/land35.asc 
## names      : land35 
## values     : -2147483648, 2147483647  (min, max)
summary(RasterPlotLand35)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         land35
## Min.         0
## 1st Qu.      0
## Median       0
## 3rd Qu.      0
## Max.         1
## NA's         0
#head(RasterPlotLandNow)
#View(RasterPlotLandNow)
RasterPlotLand35[RasterPlotLand35 > 1] <- 1
RasterPlotLand35[RasterPlotLand35 < 0] <- 0

RasterPlotLand35
## 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      : land35 
## values     : 0, 1  (min, max)
summary(RasterPlotLand35)
##         land35
## Min.         0
## 1st Qu.      0
## Median       0
## 3rd Qu.      0
## Max.         1
## NA's         0

Lithology

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

RasterPlotLith35 <- raster(paste(InDir, "Lithology35.asc", sep = ""))

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

RasterPlotLith35
## 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/Lithology35.asc 
## names      : Lithology35 
## values     : -2147483648, 2147483647  (min, max)
summary(RasterPlotLith35)
## Warning in sampleInt(length(x), size): size changed to n because it cannot be
## larger than n when replace is FALSE
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         Lithology35
## Min.              1
## 1st Qu.           1
## Median            2
## 3rd Qu.           5
## Max.              5
## NA's              0
#head(RasterPlotLandNow)
#View(RasterPlotLandNow)
RasterPlotLith35[RasterPlotLith35 > 5] <- NA
RasterPlotLith35[RasterPlotLith35 < 1] <- NA

RasterPlotLith35
## 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      : Lithology35 
## values     : 1, 5  (min, max)
summary(RasterPlotLith35)
##         Lithology35
## Min.              1
## 1st Qu.           1
## Median            2
## 3rd Qu.           5
## Max.              5
## NA's         155480

GeoMean 1 = sedimentary rocks = 96.7 2 = volcanic rocks = 67.2 3 = metamorphic rocks = 216.1 4 = plutonic rocks = 157.1 5 = unconsolidated sediments = 30.2

RasterPlotGEOAGE_Mean35 <- RasterPlotLith35

RasterPlotGEOAGE_Mean35[RasterPlotGEOAGE_Mean35 == 1] <- 96.7
RasterPlotGEOAGE_Mean35[RasterPlotGEOAGE_Mean35 == 2] <- 67.2
RasterPlotGEOAGE_Mean35[RasterPlotGEOAGE_Mean35 == 3] <- 216.1
RasterPlotGEOAGE_Mean35[RasterPlotGEOAGE_Mean35 == 4] <- 157.1
RasterPlotGEOAGE_Mean35[RasterPlotGEOAGE_Mean35 == 5] <- 30.2

RasterPlotGEOAGE_Mean35
## 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      : Lithology35 
## values     : 30.2, 216.1  (min, max)
summary(RasterPlotGEOAGE_Mean35)
##         Lithology35
## Min.           30.2
## 1st Qu.        30.2
## Median         96.7
## 3rd Qu.        96.7
## Max.          216.1
## NA's       155480.0

Temperature

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

RasterPlotTemp35 <- raster(paste(InDir, "temperatuur35.asc", sep = ""))

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

RasterPlotTemp35
## 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/temperatuur35.asc 
## names      : temperatuur35 
## values     : -2147483648, 2147483647  (min, max)
summary(RasterPlotTemp35)
## Warning in sampleInt(length(x), size): size changed to n because it cannot be
## larger than n when replace is FALSE
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         temperatuur35
## Min.               16
## 1st Qu.            26
## Median             29
## 3rd Qu.            29
## Max.               29
## NA's                0
###All temperatures -1C because it was 1 degrees colder


RasterPlotTemp35[RasterPlotTemp35 == 16] <- 15
RasterPlotTemp35[RasterPlotTemp35 == 21] <- 20
RasterPlotTemp35[RasterPlotTemp35 == 24] <- 23
RasterPlotTemp35[RasterPlotTemp35 == 26] <- 25
RasterPlotTemp35[RasterPlotTemp35 == 28] <- 27
RasterPlotTemp35[RasterPlotTemp35 == 29] <- 28

RasterPlotTemp35[RasterPlotTemp35 > 28] <- NA
RasterPlotTemp35[RasterPlotTemp35 < 15] <- NA

RasterPlotTemp35
## 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      : temperatuur35 
## values     : 15, 28  (min, max)
summary(RasterPlotTemp35)
##         temperatuur35
## Min.               15
## 1st Qu.            25
## Median             28
## 3rd Qu.            28
## Max.               28
## NA's           155283

Volcanoes

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

RasterPlotVolc35 <- raster(paste(InDir, "volcanoes35.asc", sep = ""))

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

RasterPlotVolc35
## 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/volcanoes35.asc 
## names      : volcanoes35 
## values     : -2147483648, 2147483647  (min, max)
summary(RasterPlotVolc35) 
## Warning in sampleInt(length(x), size): size changed to n because it cannot be
## larger than n when replace is FALSE
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         volcanoes35
## Min.              1
## 1st Qu.           1
## Median            1
## 3rd Qu.           1
## Max.              1
## NA's              0
#head(RasterPlotLandNow)
#View(RasterPlotLandNow)
RasterPlotVolc35[RasterPlotVolc35 > 1] <- 1
RasterPlotVolc35[RasterPlotVolc35 < 0] <- 0
RasterPlotVolc35[is.na(RasterPlotVolc35)] <- 0

RasterPlotVolc35
## 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      : volcanoes35 
## values     : 0, 1  (min, max)
summary(RasterPlotVolc35)
##         volcanoes35
## Min.              0
## 1st Qu.           0
## Median            0
## 3rd Qu.           0
## Max.              1
## NA's              0

Trenches

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

RasterPlotTrench35 <- raster(paste(InDir, "trenches35.asc", sep = ""))

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

#RasterPlotTrench35
#summary(RasterPlotTrench35)

#head(RasterPlotLandNow)
#View(RasterPlotLandNow)
RasterPlotTrench35[RasterPlotTrench35 > 1] <- 1
RasterPlotTrench35[RasterPlotTrench35 < 0] <- 0
RasterPlotTrench35[is.na(RasterPlotTrench35)] <- 0

RasterPlotTrench35
## 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      : trenches35 
## values     : 0, 1  (min, max)
summary(RasterPlotTrench35)
##         trenches35
## Min.             0
## 1st Qu.          0
## Median           0
## 3rd Qu.          0
## Max.             1
## NA's             0

Faults

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

RasterPlotFaults35 <- raster(paste(InDir, "Faults35.asc", sep = ""))

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

RasterPlotFaults35
## 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/faults35.asc 
## names      : faults35 
## values     : -2147483648, 2147483647  (min, max)
summary(RasterPlotFaults35)
## Warning in sampleInt(length(x), size): size changed to n because it cannot be
## larger than n when replace is FALSE
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         faults35
## Min.           1
## 1st Qu.        1
## Median         1
## 3rd Qu.        1
## Max.           1
## NA's           0
#head(RasterPlotLandNow)
#View(RasterPlotLandNow)
RasterPlotFaults35[RasterPlotFaults35 > 1] <- 1
RasterPlotFaults35[RasterPlotFaults35 < 0] <- 0
RasterPlotFaults35[is.na(RasterPlotFaults35)] <- 0

RasterPlotFaults35
## 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      : faults35 
## values     : 0, 1  (min, max)
summary(RasterPlotFaults35)
##         faults35
## Min.           0
## 1st Qu.        0
## Median         0
## 3rd Qu.        0
## Max.           1
## NA's           0

Plateboundary

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

RasterPlotPlatebound35 <- raster(paste(InDir, "Platebound35.asc", sep = ""))

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

RasterPlotPlatebound35
## 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/Platebound35.asc 
## names      : Platebound35 
## values     : -2147483648, 2147483647  (min, max)
summary(RasterPlotPlatebound35)
## Warning in sampleInt(length(x), size): size changed to n because it cannot be
## larger than n when replace is FALSE
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         Platebound35
## Min.               1
## 1st Qu.            1
## Median             1
## 3rd Qu.            1
## Max.               1
## NA's               0
RasterPlotPlatebound35[RasterPlotPlatebound35 > 1] <- 1
RasterPlotPlatebound35[RasterPlotPlatebound35 < 0] <- 0
RasterPlotPlatebound35[is.na(RasterPlotPlatebound35)] <- 0

RasterPlotPlatebound35
## 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      : Platebound35 
## values     : 0, 1  (min, max)
summary(RasterPlotPlatebound35)
##         Platebound35
## Min.               0
## 1st Qu.            0
## Median             0
## 3rd Qu.            0
## Max.               1
## NA's               0

Transform faults

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

RasterPlotTransform35 <- raster(paste(InDir, "Transform35.asc", sep = ""))

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

RasterPlotTransform35
## 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/Transform35.asc 
## names      : Transform35 
## values     : -2147483648, 2147483647  (min, max)
summary(RasterPlotTransform35)
## Warning in sampleInt(length(x), size): size changed to n because it cannot be
## larger than n when replace is FALSE
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         Transform35
## Min.              1
## 1st Qu.           1
## Median            1
## 3rd Qu.           1
## Max.              1
## NA's              0
RasterPlotTransform35[RasterPlotTransform35 > 1] <- 1
RasterPlotTransform35[RasterPlotTransform35 < 0] <- 0
RasterPlotTransform35[is.na(RasterPlotTransform35)] <- 0

RasterPlotTransform35
## 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      : Transform35 
## values     : 0, 1  (min, max)
summary(RasterPlotTransform35)
##         Transform35
## Min.              0
## 1st Qu.           0
## Median            0
## 3rd Qu.           0
## Max.              1
## NA's              0

EucDistVolc

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

RPEucDistVolc35 <- raster(paste(InDir, "Eucdistvolc35.asc", sep = ""))

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

RPEucDistVolc35
## 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/Eucdistvolc35.asc 
## names      : Eucdistvolc35
summary(RPEucDistVolc35)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         Eucdistvolc35
## Min.              0.0
## 1st Qu.      367100.5
## Median       801977.3
## 3rd Qu.     1436154.2
## Max.        3120683.0
## NA's              0.0

EucDistTrench

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

RPEucDistTrench35 <- raster(paste(InDir, "Eucdisttrench35.asc", sep = ""))

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

RPEucDistTrench35
## 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/Eucdisttrench35.asc 
## names      : Eucdisttrench35
summary(RPEucDistTrench35)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         Eucdisttrench35
## Min.                0.0
## 1st Qu.        267709.3
## Median         554284.3
## 3rd Qu.        873824.6
## Max.          1849989.0
## NA's                0.0

EucDistTransform

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

RPEucDistTransform35 <- raster(paste(InDir, "Eucdisttransform35.asc", sep = ""))

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

RPEucDistTransform35
## 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/Eucdisttransform35.asc 
## names      : Eucdisttransform35
summary(RPEucDistTransform35)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         Eucdisttransform35
## Min.                     0
## 1st Qu.            1019413
## Median             1663219
## 3rd Qu.            2637281
## Max.               4491634
## NA's                     0

EucDistPlateBound

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

RPEucDistPlateBound35 <- raster(paste(InDir, "Eucdistplatebound35.asc", sep = ""))

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

RPEucDistPlateBound35
## 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/Eucdistplatebound35.asc 
## names      : Eucdistplatebound35
summary(RPEucDistPlateBound35)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         Eucdistplatebound35
## Min.                    0.0
## 1st Qu.            214560.3
## Median             472205.4
## 3rd Qu.            788811.4
## Max.              1823319.0
## NA's                    0.0

EucDistFaults

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

RPEucDistFaults35 <- raster(paste(InDir, "Eucdistfaults35.asc", sep = ""))

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

RPEucDistFaults35
## 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/Eucdistfaults35.asc 
## names      : Eucdistfaults35
summary(RPEucDistFaults35)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (55.56% of all cells)
##         Eucdistfaults35
## Min.                0.0
## 1st Qu.        907994.5
## Median        1900511.5
## 3rd Qu.       3439708.2
## Max.          5395463.0
## NA's                0.0

Precipitation

RasterPlotPrecip35 <- RasterPlotVolc35
RasterPlotPrecip35[RasterPlotPrecip35 >= 1] <- NA
RasterPlotPrecip35[RasterPlotPrecip35 < 1] <- NA
RasterPlotPrecip35
## 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      : volcanoes35 
## values     : NA, NA  (min, max)
summary(RasterPlotPrecip35)
##         volcanoes35
## Min.             NA
## 1st Qu.          NA
## Median           NA
## 3rd Qu.          NA
## Max.             NA
## NA's         180000

Plot all paleo-variables

plot(RasterPlotLand35, main = 'Land35')

plot(RasterPlotGEOAGE_Mean35, main = 'GeoAge mean 35')

plot(RasterPlotLith35, main = 'Lithology 35')

plot(RasterPlotTemp35, main = 'Temp 35')

plot(RPEucDistTransform35, main = 'EucDist Transform35')

plot(RPEucDistFaults35, main = 'EucDist Faults35')

plot(RasterPlotFaults35, main = 'Faults35')

plot(RPEucDistPlateBound35, main = 'EucDist PlateBound35')

# No precipitation 
plot(RPEucDistTrench35, main = 'EucDist Trench35')

plot(RPEucDistVolc35, main = 'EucDist Volc35')

plot(RasterPlotVolc35, main = 'Volcanoes35')

plot(RasterPlotTransform35, main = 'Transform35')

plot(RasterPlotTrench35, main = 'Trench35')

plot(RasterPlotPlatebound35, main = 'Platebound35')

##Make dataframe

Make brick of all variables

Paleo35Good_Stack <- stack(RasterPlotLand35,RasterPlotGEOAGE_Mean35, RasterPlotLith35, RasterPlotTemp35, RPEucDistTransform35, RPEucDistFaults35, RasterPlotFaults35, RPEucDistPlateBound35,                    RasterPlotPrecip35, RPEucDistTrench35, RPEucDistVolc35, RasterPlotVolc35,RasterPlotTransform35,RasterPlotTrench35, RasterPlotPlatebound35) 

Paleo35Good_Brick <- brick(Paleo35Good_Stack)
Paleo35Good_Brick
## class      : RasterBrick 
## dimensions : 300, 600, 180000, 15  (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/RtmpOu8Wjj/raster/r_tmp_2021-06-30_132155_15648_25222.grd 
## names      :    land35, Lithology35.1, Lithology35.2, temperatuur35, Eucdisttransform35, Eucdistfaults35,  faults35, Eucdistplatebound35, volcanoes35.1, Eucdisttrench35, Eucdistvolc35, volcanoes35.2, Transform35, trenches35, Platebound35 
## min values :       0.0,          30.2,           1.0,          15.0,                0.0,             0.0,       0.0,                 0.0,            NA,             0.0,           0.0,           0.0,         0.0,        0.0,          0.0 
## max values :       1.0,         216.1,           5.0,          28.0,          4495850.0,       5395463.0,       1.0,           1824767.0,            NA,       1849989.0,     3130681.0,           1.0,         1.0,        1.0,          1.0
Paleo35Good_Brick_df <- as.data.frame(Paleo35Good_Brick)
#Paleo35Good_Brick_df

###Change the names of each variable Make sure they are named the exact same way as in EurekaAll_Brick_df dataframe

names(Paleo35Good_Brick_df)
##  [1] "land35"              "Lithology35.1"       "Lithology35.2"      
##  [4] "temperatuur35"       "Eucdisttransform35"  "Eucdistfaults35"    
##  [7] "faults35"            "Eucdistplatebound35" "volcanoes35.1"      
## [10] "Eucdisttrench35"     "Eucdistvolc35"       "volcanoes35.2"      
## [13] "Transform35"         "trenches35"          "Platebound35"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "land35"] <- "Land"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "Lithology35.1"] <- "GeoMean"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "Lithology35.2"] <- "Lithology"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "temperatuur35"] <- "Temp"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "Eucdisttransform35"] <- "DistTrans"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "Eucdistfaults35"] <- "DistFault"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "faults35"] <- "Faults"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "Eucdistplatebound35"] <- "DistPlateBound"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "volcanoes35.1"] <- "Precip"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "Eucdisttrench35"] <- "DistTrench"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "Eucdistvolc35"] <- "DistVolc"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "volcanoes35.2"] <- "Volc"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "Transform35"] <- "Trans"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "trenches35"] <- "Trench"
names(Paleo35Good_Brick_df)[names(Paleo35Good_Brick_df) == "Platebound35"] <- "PlateBound"

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

See head of dataframe

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

RECONSTRUCTION 35MA

Predict35_RpartAll <- predict(Rpart.EurekaAll.Final, newdata = Paleo35Good_Brick_df, type = "class")
Predict35_RpartAll_df <- as.data.frame(Predict35_RpartAll)

Get the rownames of data subsets

rownames_pred35_RpartAll <- rownames(Predict35_RpartAll_df)
rownames_pred35_numeric_RpartAll <- as.numeric(rownames_pred35_RpartAll)

Get coordinates of each cell in dataframe

Coordinates <- coordinates(EurekaAll_Brick)
Coordinates_df <- as.data.frame(Coordinates)

Get the coordinates of the predicted datasubsets

XCor_pred35_RpartAll <- Coordinates_df$x[rownames_pred35_numeric_RpartAll]
YCor_pred35_RpartAll <- Coordinates_df$y[rownames_pred35_numeric_RpartAll]

Get the coordinates into one file.

XYCor_pred35_RpartAll <- data.frame(XCor_pred35_RpartAll, YCor_pred35_RpartAll)

Get the coordinates and the 0s and 1s into one file.

XYZ_pred35_RpartAll <- data.frame(XYCor_pred35_RpartAll,Predict35_RpartAll)

Into raster

RasterPlot_pred35_RpartAll <- rasterFromXYZ(XYZ_pred35_RpartAll, res=c(0.1,0.1), crs="+proj=longlat +datum=WGS84 +no_defs ")

Plot

plot(RasterPlot_pred35_RpartAll, main="Predicted35")

To GeoTIFF

writeRaster(RasterPlot_pred35_RpartAll, "RasterPlot_Pred35_RpartAll", format ="GTiff", overwrite = TRUE)
## Warning in .gd_SetProject(object, ...): NOT UPDATED FOR PROJ >= 6

Make pdf

pdf('Reconstructie 35Ma')
plot(RasterPlot_pred35_RpartAll, main="Predicted reconstruction 35Ma")
dev.off()
## png 
##   2

Calculate area of 1s

summary(Predict35_RpartAll_df)
##  Predict35_RpartAll
##  0:178643          
##  1:  1357