-
Notifications
You must be signed in to change notification settings - Fork 1
/
0.Prepare_data_dominant_AUV.R
139 lines (101 loc) · 3.65 KB
/
0.Prepare_data_dominant_AUV.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
#################################################
############ 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"
#### DOMINANT DATA ----
### Load data ----
# Habitat data --
df <- read.csv(paste(d.dir, "raw", "dominant-auv_broad.percent.cover.csv", sep='/'))
str(df) # 7 categories
# Bathy 250m all Geo Bay --
bgb <- raster(paste(s.dir, "GB_Bathy_250m.tif", sep='/'))
bgb
plot(bgb)
#bx <- raster(paste(s.dir, "GBmultib_lidarUTM_CMR.tif", sep='/'))
#bx # in UTM
#plot(bx)
# 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='/'))
# Bathy Lidar and Multibeam in latlong --
b <- raster(paste(s.dir, "GBmultib_lidar_CMR.tif", sep='/'))
plot(b)
b
# 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)
str(hab_pred) # 837 observations
names(hab_pred)
## remove unwanted columns ----
auv_dom <- hab_pred[,c(2:4,14, 16:24)] # not sure what to remove for AUV
names(auv_dom)
sp1 <- auv_dom
coordinates(sp1) <- ~longitude+latitude
plot(b)
#plot(sp1, border="white", col="lightgrey", add=TRUE)
plot(sp1, col=rainbow(7), pch=20, fill=sp1$Max_if_2_habitats_have_same, add=TRUE)
# save
writeOGR(sp1, dsn= s.dir, layer= "GB_auv_fine_bathy_habitat_dominant_broad", driver="ESRI Shapefile", overwrite_layer=TRUE)
write.csv(auv_dom, paste(d.dir, "tidy", "GB_auv_fine_bathy_habitat_dominant_broad.csv", sep='/'))