Load packages
library(sp)
library(raster)
library(rpart)
library(DT)
Indicate input folder
InDir <- 'C:/Users/Linde/Desktop/Scriptie/Data per script/30MaScriptFinal (3)/'
load("C:/Users/Linde/Desktop/Scriptie/Data per script/30MaScriptFinal (3)/EurekaAll_Brick_df.RData")
Check reconstruction difference if removing precipitation 30Ma
See head of dataframe
datatable(head(Paleo30Good_Brick_df[1:5,],20))
Paleo30Good_Brick_df_NOPRECIP <- Paleo30Good_Brick_df
Paleo30Good_Brick_df_NOPRECIP$Precip[Paleo30Good_Brick_df_NOPRECIP$Precip >= 319] <- NA
summary(Paleo30Good_Brick_df_NOPRECIP$Precip);
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## NA NA NA NaN NA NA 180000
Predict30_RpartAll_NOPRECIP <- predict(Rpart.EurekaAll.Final, Paleo30Good_Brick_df_NOPRECIP, type = "class");
Predict30_RpartAll_df_NOPRECIP <- as.data.frame(Predict30_RpartAll_NOPRECIP);
summary(Predict30_RpartAll_df)
## Predict30_RpartAll
## 0:176144
## 1: 3856
summary(Predict30_RpartAll_df_NOPRECIP)
## Predict30_RpartAll_NOPRECIP
## 0:176371
## 1: 3629
RECONSTRUCTION 30MA NO PRECIPITATION
Predict30_RpartAll_NOPRECIP <- predict(Rpart.EurekaAll.Final, Paleo30Good_Brick_df_NOPRECIP, type = "class")
Predict30_RpartAll_df_NOPRECIP <- as.data.frame(Predict30_RpartAll_NOPRECIP)
Get the rownames of data subsets
rownames_pred30_RpartAll_NOPRECIP <- rownames(Predict30_RpartAll_df_NOPRECIP)
rownames_pred30_numeric_RpartAll_NOPRECIP <- as.numeric(rownames_pred30_RpartAll_NOPRECIP)
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_NOPRECIP <- Coordinates_df$x[rownames_pred30_numeric_RpartAll_NOPRECIP]
YCor_pred30_RpartAll_NOPRECIP <- Coordinates_df$y[rownames_pred30_numeric_RpartAll_NOPRECIP]
Get the coordinates into one file.
XYCor_pred30_RpartAll_NOPRECIP <- data.frame(XCor_pred30_RpartAll_NOPRECIP, YCor_pred30_RpartAll_NOPRECIP)
Get the coordinates and the 0s and 1s into one file.
XYZ_pred30_RpartAll_NOPRECIP <- data.frame(XYCor_pred30_RpartAll_NOPRECIP,Predict30_RpartAll_NOPRECIP)
Into raster
RasterPlot_pred30_RpartAll_NOPRECIP <- rasterFromXYZ(XYZ_pred30_RpartAll_NOPRECIP, res=c(0.1,0.1), crs="+proj=longlat +datum=WGS84 +no_defs ")
plot both
plot(RasterPlot_pred30_RpartAll, main="Predicted30")
plot(RasterPlot_pred30_RpartAll_NOPRECIP, main="Predicted30 NOPRECIP")