forked from anitas-giraldo/GB_Habitat_Classification
-
Notifications
You must be signed in to change notification settings - Fork 0
/
0.Prepare_data_pres-abs_coarse_BRUV.R
142 lines (103 loc) · 3.88 KB
/
0.Prepare_data_pres-abs_coarse_BRUV.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
#################################################
############ Prepare data for RF ################
### Load libraries --
library(ggplot2)
library(cowplot)
library(sp)
library(rgdal)
library(raster)
# Clear memory ----
rm(list=ls())
### Set directories ----
w.dir <- "Y:/GB_Habitat_Classification"
d.dir <- "Y:/GB_Habitat_Classification/data"
s.dir <- "Y:/GB_Habitat_Classification/spatial_data"
p.dir <- "Y:/GB_Habitat_Classification/plots"
o.dir <- "Y:/GB_Habitat_Classification/outputs"
#### PRESENCE-ABSENCE DATA ----
### Load data ----
# Habitat data --
df <- read.csv(paste(d.dir, "raw", "stereo-BRUVs_broad.presence.absence.csv", sep='/'))
str(df) # 7 categories
# Bathy 250m all Geo Bay --
b <- raster(paste(s.dir, "GB_Bathy_250m.tif", sep='/'))
b
plot(b)
#b <- raster(paste(s.dir, "GBmultib_lidarUTM_CMR.tif", sep='/'))
#b # in UTM
#plot(b)
# Reproject lidat and multibeam data
#b2 <- projectRaster(b, crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
#b2
#plot(b2)
# save
#writeRaster(b2, paste(s.dir, "GBmultib_lidar_CMR.tif", sep='/'))
# Check data points --
dfs <- df
coordinates(dfs) <- ~longitude+latitude
points(dfs)
# Filter data to use just where there is lidar and multib --
dnew <- raster::extract(b, dfs, sp=T)
str(dnew)
dfnew <- as.data.frame(dnew)
dfnew <- na.omit(dfnew)
str(dfnew)
dfsnew <- dfnew
coordinates(dfsnew) <- ~longitude+latitude
points(dfsnew, pch = 20, col="red")
h <- dfsnew
# save points --
#write.csv(dfnew, paste(d.dir, "GB_Bruvs_fine_bathy_habitat_presence_absence_broad.csv"))
#writeOGR(dfsnew, dsn= s.dir, layer= "GB_Bruvs_fine_bathy_habitat_presence_absence_broad", driver="ESRI Shapefile", overwrite_layer=TRUE)
## load shape file with habitat data ----
#h <- readOGR(paste(s.dir, "GB_Bruvs_fine_bathy_habitat_presence_absence_broad.shp", sep='/'))
######################################
#### Calculate Bathy Derivatives ####
slope4 <- raster::terrain(b, "slope", unit="degrees", neighboors=4)
slope4@data@names
slope4@data@names <- "slope4"
slope8 <- raster::terrain(b, "slope", unit="degrees", neighboors=8)
slope8@data@names
slope8@data@names <- "slope8"
aspect4 <- raster::terrain(b, "aspect", unit="degrees", neighboors=4)
aspect4@data@names
aspect4@data@names <- "aspect4"
aspect8 <- raster::terrain(b, "aspect", unit="degrees", neighboors=8)
aspect8@data@names
aspect8@data@names <- "aspect8"
tpi <- raster::terrain(b, "TPI")
tri <- raster::terrain(b, "TRI")
roughness <- raster::terrain(b, "roughness")
flowdir <- raster::terrain(b, "flowdir")
b@data@names <- "depth"
predictors <- stack(b, slope4, slope8, aspect4, aspect8, tpi, tri, roughness, flowdir)
plot(predictors)
names(predictors)
namesp <- names(predictors)
# save
#writeRaster(predictors, paste(s.dir,"predictors.tif", sep='/'), overwrite=T)
#write.csv(namesp, paste(s.dir, "namespredictors.csv", sep='/'))
### Extract predictor values for each observation ----
hp <- raster::extract(predictors, h, sp=T)
hp
hab_pred <- as.data.frame(hp)
head(hab_pred)
## make long data frame where hab classes are factors of one column ----
library(reshape2)
# choose ID vars
ids <- c("X", "campaignid", "sample", "depth", "slope4", "slope8", "aspect4", "aspect8","tpi",
"tri", "roughness", "flowdir", "latitude", "longitude")
Bruv.broad <- melt(hab_pred, id.vars=ids)
head(Bruv.broad)
names(Bruv.broad)[names(Bruv.broad)=="variable"] <- "Class"
names(Bruv.broad)[names(Bruv.broad)=="value"] <- "presence.absence"
str(Bruv.broad)
sp1 <- Bruv.broad
coordinates(sp1) <- ~longitude+latitude
plot(b)
plot(sp1, border="white", col="lightgrey", add=TRUE)
plot(sp1, col=rainbow(7), pch=20, fill=sp1$Class, add=TRUE)
# save
writeOGR(sp1, dsn= s.dir, layer= "GB_Bruvs_fine_bathy_habitat_presence_absence_broad", driver="ESRI Shapefile", overwrite_layer=TRUE)
write.csv(Bruv.broad, paste(d.dir, "tidy", "GB_Bruvs_fine_bathy_habitat_presence_absence_broad.csv", sep='/'))
### DOMINANT CLASS DATA ----