Install packages

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

Load packages

library(sp)
library(raster)
library(rpart)
library(rasclass)
## 
## Attaching package: 'rasclass'
## The following object is masked from 'package:raster':
## 
##     writeRaster
library(rgeos)
## rgeos version: 0.5-3, (SVN revision 634)
##  GEOS runtime version: 3.8.0-CAPI-1.13.1 
##  Linking to sp version: 1.4-2 
##  Polygon checking: TRUE
library(foreign)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:rgeos':
## 
##     intersect, setdiff, union
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
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

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

RasterPlotDEM_Binair30 <- raster(paste(InDir, "dem30_binary_reclassfromlith.asc", sep = ""))

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

#RasterPlotDEM_Binair30

RasterPlotDEM_Binair30[RasterPlotDEM_Binair30 >= 1] <- 1
RasterPlotDEM_Binair30[RasterPlotDEM_Binair30 <= 0] <- 0
RasterPlotDEM_Binair30[is.na(RasterPlotDEM_Binair30)] <- 0

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

Lithology

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

RasterPlotLITOLOGY30 <- raster(paste(InDir, "lithology30.asc", sep = ""))

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

#RasterPlotLITOLOGY30

RasterPlotLITOLOGY30[RasterPlotLITOLOGY30 > 5] <- NA
RasterPlotLITOLOGY30[RasterPlotLITOLOGY30 < 1] <- NA

RasterPlotLITOLOGY30
## 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      : lithology30 
## values     : 1, 5  (min, max)
summary(RasterPlotLITOLOGY30)
##         lithology30
## Min.              1
## 1st Qu.           1
## Median            2
## 3rd Qu.           5
## Max.              5
## NA's         151439
#head(RasterPlotLITOLOGY30)
#View(RasterPlotLITOLOGY30)

GEOAGE Making GeoMean30 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_Mean30 <- RasterPlotLITOLOGY30

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

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

To GeoTIFF

#writeRaster(RasterPlotGEOAGE_Mean30, "RasterPlot_GeoAge_Mean30", format ="GTiff", overwrite = TRUE)

Temperature

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

RasterPlotTEMP30good <- raster(paste(InDir, "temperature30final.asc", sep = ""))

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

#RasterPlotTEMP30good
#summary(RasterPlotTEMP30good)

#> unique(RasterPlotTEMP30good)
# 16 21 24 26 28 29

###All temperatures -2C because it was 2 degrees colder

RasterPlotTEMP30good[RasterPlotTEMP30good == 16] <- 14
RasterPlotTEMP30good[RasterPlotTEMP30good == 21] <- 19
RasterPlotTEMP30good[RasterPlotTEMP30good == 24] <- 22
RasterPlotTEMP30good[RasterPlotTEMP30good == 26] <- 24 
RasterPlotTEMP30good[RasterPlotTEMP30good == 28] <- 26 
RasterPlotTEMP30good[RasterPlotTEMP30good == 29] <- 27 

RasterPlotTEMP30good[RasterPlotTEMP30good > 27] <- NA 
RasterPlotTEMP30good[RasterPlotTEMP30good < 14] <- NA 

RasterPlotTEMP30good
## 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      : temperature30final 
## values     : 14, 27  (min, max)
summary(RasterPlotTEMP30good)
##         temperature30final
## Min.                    14
## 1st Qu.                 24
## Median                  27
## 3rd Qu.                 27
## Max.                    27
## NA's                151873

Precipitation

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

RasterPlotPRECIP30good <- raster(paste(InDir, "precip30good.asc", sep = ""))

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

#RasterPlotPRECIP30good
#summary(RasterPlotPRECIP30good)

RasterPlotPRECIP30good[RasterPlotPRECIP30good > 6388] <- NA
RasterPlotPRECIP30good[RasterPlotPRECIP30good < 319] <- NA

RasterPlotPRECIP30good
## 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      : precip30good 
## values     : 319, 6388  (min, max)
summary(RasterPlotPRECIP30good)
##         precip30good
## Min.             319
## 1st Qu.         2395
## Median          3194
## 3rd Qu.         3194
## Max.            6388
## NA's          151445

RPEucDistTransform

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

RPEucDistTransform30 <- raster(paste(InDir, "transform30buffer.asc", sep = ""))

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

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

RPEucDistFaults

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

RPEucDistFaults30 <- raster(paste(InDir, "faults30buffer.asc", sep = ""))

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

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

Faults

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

RasterPlotFaults30 <- raster(paste(InDir, "faults30.asc", sep = ""))

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

#RasterPlotFaults30
#summary(RasterPlotFaults30)
#head(RasterPlotFaults30)
#View(RasterPlotFaults30)


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

RasterPlotFaults30
## 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      : faults30 
## values     : 1, 1  (min, max)
summary(RasterPlotFaults30)
##         faults30
## Min.           1
## 1st Qu.        1
## Median         1
## 3rd Qu.        1
## Max.           1
## NA's      179696
RasterPlotFaults30_Buffer <- buffer(RasterPlotFaults30, width = 20000)
RasterPlotFaults30_Buffer[is.na(RasterPlotFaults30_Buffer)] <- 0

RasterPlotFaults30_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)
summary(RasterPlotFaults30_Buffer)
##         layer
## Min.        0
## 1st Qu.     0
## Median      0
## 3rd Qu.     0
## Max.        1
## NA's        0

RPEucDistPlateBound

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

RPEucDistPlateBound30 <- raster(paste(InDir, "plateboundary30buffer.asc", sep = ""))

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

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

RPEucDistTrench

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

RPEucDistTrench30 <- raster(paste(InDir, "trenches30buffer.asc", sep = ""))

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

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

RPEucDistVolc

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

RPEucDistVolc30 <- raster(paste(InDir, "volc30buffer.asc", sep = ""))

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

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

Volcanoes

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

RasterPlotVOLC30_Buffer <- raster(paste(InDir, "volcanoes30buffer2.asc", sep = ""))

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

RasterPlotVOLC30_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     : C:/Users/Linde/Desktop/Scriptie/Schrijven begin/Textuele conceptversie/Final Conceptversion MSc Linde Berbers/Supplementary Information/Supplementary Data 2. RStudio Data/Script Data/volcanoes30buffer2.asc 
## names      : volcanoes30buffer2 
## values     : -2147483648, 2147483647  (min, max)
summary(RasterPlotVOLC30_Buffer)
## 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)
##         volcanoes30buffer2
## Min.                     1
## 1st Qu.                 11
## Median                  22
## 3rd Qu.                 32
## Max.                    42
## NA's                     0
#head(RasterPlotFaults30)
#View(RasterPlotFaults30)


RasterPlotVOLC30_Buffer[RasterPlotVOLC30_Buffer >= 1] <- 1
RasterPlotVOLC30_Buffer[is.na(RasterPlotVOLC30_Buffer)] <- 0

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

Transform faults

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

RasterPlotTransform30 <- raster(paste(InDir, "transform30.asc", sep = ""))

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

#RasterPlotTransform30
#summary(RasterPlotTransform30)
#head(RasterPlotTransform30)
#View(RasterPlotTransform30)


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

RasterPlotTransform30
## 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      : transform30 
## values     : 1, 1  (min, max)
summary(RasterPlotTransform30)
##         transform30
## Min.              1
## 1st Qu.           1
## Median            1
## 3rd Qu.           1
## Max.              1
## NA's         179948
RasterPlotTransform30_Buffer <- buffer(RasterPlotTransform30, width = 20000)
RasterPlotTransform30_Buffer[is.na(RasterPlotTransform30_Buffer)] <- 0

RasterPlotTransform30_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)
summary(RasterPlotTransform30_Buffer)
##         layer
## Min.        0
## 1st Qu.     0
## Median      0
## 3rd Qu.     0
## Max.        1
## NA's        0

Trenches

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

RasterPlotTrench30 <- raster(paste(InDir, "trenches30.2.asc", sep = ""))

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

#RasterPlotTrench30
#summary(RasterPlotTrench30)

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

#RasterPlotTrench30
#summary(RasterPlotTrench30)

RasterPlotTrench30_Buffer <- buffer(RasterPlotTrench30, width = 20000)
RasterPlotTrench30_Buffer[is.na(RasterPlotTrench30_Buffer)] <- 0

RasterPlotTrench30_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)
summary(RasterPlotTrench30_Buffer)
##         layer
## Min.        0
## 1st Qu.     0
## Median      0
## 3rd Qu.     0
## Max.        1
## NA's        0

Plateboundaries

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

RasterPlotPlateBoundaries30 <- raster(paste(InDir, "plateboundaries30.asc", sep = ""))

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

#RasterPlotPlateBoundaries30
#summary(RasterPlotPlateBoundaries30)

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

RasterPlotPlateBoundaries30
## 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      : plateboundaries30 
## values     : 1, 1  (min, max)
summary(RasterPlotPlateBoundaries30)
##         plateboundaries30
## Min.                    1
## 1st Qu.                 1
## Median                  1
## 3rd Qu.                 1
## Max.                    1
## NA's               178494
RasterPlotPlateBoundaries30_Buffer <- buffer(RasterPlotPlateBoundaries30, width = 20000)
RasterPlotPlateBoundaries30_Buffer[is.na(RasterPlotPlateBoundaries30_Buffer)] <- 0

RasterPlotPlateBoundaries30_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)
summary(RasterPlotPlateBoundaries30_Buffer)
##         layer
## Min.        0
## 1st Qu.     0
## Median      0
## 3rd Qu.     0
## Max.        1
## NA's        0

Plot all paleo-variables

plot(RasterPlotDEM_Binair30, main = 'Land 30')

plot(RasterPlotGEOAGE_Mean30, main = 'GeoMean30')

plot(RasterPlotLITOLOGY30, main = 'Lithology30')

plot(RasterPlotTEMP30good, main = 'Temp30good')

plot(RPEucDistTransform30, main = 'Distance Transform30')

plot(RPEucDistFaults30, main = 'Distance Faults30')

plot(RasterPlotFaults30_Buffer, main = 'Faults30')

plot(RPEucDistPlateBound30, main = 'Distance Platebound30')

plot(RasterPlotPRECIP30good, main = 'Precip 30 good')

plot(RPEucDistTrench30, main = 'Distance Trench30')

plot(RPEucDistVolc30, main = 'Distance Volc30')

plot(RasterPlotVOLC30_Buffer, main = 'Volc 30')

plot(RasterPlotTransform30_Buffer, main = 'Transform 30')

plot(RasterPlotTrench30_Buffer, main = 'Trench30')

plot(RasterPlotPlateBoundaries30_Buffer, main = 'Plateboundaries 30')

Make dataframe

Make brick of all variables

Paleo30Good_Stack <- stack(RasterPlotDEM_Binair30,RasterPlotGEOAGE_Mean30, RasterPlotLITOLOGY30, RasterPlotTEMP30good, RPEucDistTransform30, RPEucDistFaults30, RasterPlotFaults30_Buffer, RPEucDistPlateBound30, RasterPlotPRECIP30good, RPEucDistTrench30, RPEucDistVolc30, RasterPlotVOLC30_Buffer,RasterPlotTransform30_Buffer,RasterPlotTrench30_Buffer, RasterPlotPlateBoundaries30_Buffer) 


Paleo30Good_Brick <- brick(Paleo30Good_Stack)
Paleo30Good_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/RtmpG8GR6S/raster/r_tmp_2021-06-30_140619_17648_92635.grd 
## names      : dem30_bin//ssfromlith, lithology30.1, lithology30.2, temperature30final, transform30buffer, faults30buffer,   layer.1, plateboundary30buffer, precip30good, trenches30buffer, volc30buffer, volcanoes30buffer2,   layer.2,   layer.3,   layer.4 
## min values :                   0.0,          30.2,           1.0,               14.0,               0.0,            0.0,       0.0,                   0.0,        319.0,              0.0,          0.0,                0.0,       0.0,       0.0,       0.0 
## max values :                   1.0,         216.1,           5.0,               27.0,         5319470.0,      5689126.0,       1.0,             1811201.0,       6388.0,        2272025.0,    3310592.0,                1.0,       1.0,       1.0,       1.0
Paleo30Good_Brick_df <- as.data.frame(Paleo30Good_Brick)
#Paleo30Good_Brick_df

Change the names of each variable

Make sure they are named the exact same way as in EurekaAll_Brick_df dataframe

names(Paleo30Good_Brick_df)
##  [1] "dem30_binary_reclassfromlith" "lithology30.1"               
##  [3] "lithology30.2"                "temperature30final"          
##  [5] "transform30buffer"            "faults30buffer"              
##  [7] "layer.1"                      "plateboundary30buffer"       
##  [9] "precip30good"                 "trenches30buffer"            
## [11] "volc30buffer"                 "volcanoes30buffer2"          
## [13] "layer.2"                      "layer.3"                     
## [15] "layer.4"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "dem30_binary_reclassfromlith"] <- "Land"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "lithology30.1"] <- "GeoMean"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "lithology30.2"] <- "Lithology"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "temperature30final"] <- "Temp"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "transform30buffer"] <- "DistTrans"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "faults30buffer"] <- "DistFault"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "layer.1"] <- "Faults"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "plateboundary30buffer"] <- "DistPlateBound"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "precip30good"] <- "Precip"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "trenches30buffer"] <- "DistTrench"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "volc30buffer"] <- "DistVolc"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "volcanoes30buffer2"] <- "Volc"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "layer.2"] <- "Trans"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "layer.3"] <- "Trench"
names(Paleo30Good_Brick_df)[names(Paleo30Good_Brick_df) == "layer.4"] <- "PlateBound"

names(Paleo30Good_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(Paleo30Good_Brick_df,10),options = list(pageLength=10, pageWidth = 5), autoHideNavigation = TRUE)

RECONSTRUCTION 30MA

Predict30_RpartAll <- predict(Rpart.EurekaAll.Final, Paleo30Good_Brick_df, type = "class")
Predict30_RpartAll_df <- as.data.frame(Predict30_RpartAll)

Get the rownames of data subsets

rownames_pred30_RpartAll <- rownames(Predict30_RpartAll_df)
rownames_pred30_numeric_RpartAll <- as.numeric(rownames_pred30_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_pred30_RpartAll <- Coordinates_df$x[rownames_pred30_numeric_RpartAll]
YCor_pred30_RpartAll <- Coordinates_df$y[rownames_pred30_numeric_RpartAll]

Get the coordinates into one file.

XYCor_pred30_RpartAll <- data.frame(XCor_pred30_RpartAll, YCor_pred30_RpartAll)

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

XYZ_pred30_RpartAll <- data.frame(XYCor_pred30_RpartAll,Predict30_RpartAll)

Into raster

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

Plot

plot(RasterPlot_pred30_RpartAll, main="Predicted30")

#1. Ratify your RasterLayer. I assume you'll be wanting some kind of count data associated with each value in your layer, so when you use the ratify function, be sure to set count = TRUE.
px <- ratify(RasterPlot_pred30_RpartAll, count = TRUE)

# 2. Install package foreign. This package will allow you to read and write .dbf files.
# 3 Write the RAT from your RasterLayer object to a .tif.vat.dbf file using the write.dbf function. If px is a single-layer raster (not a RasterBrick), then the RAT you created with the ratify function is the first list element of levels(px), i.e. levels(px)[[1]].
write.dbf(levels(px)[[1]], file = "predict30.tif.vat.dbf")

Calculate area of 1s

summary(Predict30_RpartAll_df)
##  Predict30_RpartAll
##  0:176120          
##  1:  3880