-
Notifications
You must be signed in to change notification settings - Fork 0
/
Start_here_Import_Rast_Reclass_firescars_edits.R
313 lines (254 loc) · 11.3 KB
/
Start_here_Import_Rast_Reclass_firescars_edits.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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
## Kimberley fire stats scripts
library(sf)
library(tidyverse)
library(fasterize)
library(raster)
library(fs)
library(janitor)
library(here)
library(igraph)
# Setup folders for layers for stats --------------------------------------
st <- here::here("stacks")
ff <- here::here("fireFreq")
ac <- here::here("ageClass")
unburnt <- here::here("unburnt")
acdist <- here::here("ageClass_dist")
tslb <- here::here("tslb")
outstats <- here::here("output_stats")
ifelse(!dir.exists(st),
dir.create(st), FALSE)
ifelse(!dir.exists(ff),
dir.create(ff), FALSE)
ifelse(!dir.exists(ac),
dir.create(ac), FALSE)
ifelse(!dir.exists(unburnt),
dir.create(unburnt), FALSE)
ifelse(!dir.exists(acdist),
dir.create(acdist), FALSE)
ifelse(!dir.exists(tslb),
dir.create(tslb), FALSE)
ifelse(!dir.exists(outstats),
dir.create(outstats), FALSE)
# Set up area of interest (aoi) -------------------------------------------
## define a Area of Interest by importing a shp file of the Kimberley and transforming it to a projected CRS
aoi <- st_read("../Masks/nafi_aoi.shp") %>% ##from GP notes, nafi_aoi.shp is simply kimbreg_ext raster converted to shp. It includes the buffer area
st_transform(3577)
# Set up raster template --------------------------------------------------
#250m resolution. 1 pixel is 6.25ha. 250*250 m 62500 msq. 0.0001ha in a msq.
tmpltrast <- raster(aoi, res = 250, vals = 1)
# Kimberley coast to use as mask ------------------------------------------
##land is 1 and offshore is NA
kimbcoast_alb <- st_read("../Masks/Kimberley_Coast_CLEAN.shp") %>% ##from GP notes, Kimberley_Coast_CLEAN.shp is used to produced kimbcoast
st_transform(3577) %>%
fasterize(tmpltrast, field = "GRIDCODE")
# Set up classification matrices ------------------------------------------
## EDS - early dry season (Jan to June) ###EDS definition may change (Jan - July)
eds <- c(-1,0,0, 0,7,1, 7,12,0)
edsmat <- matrix(eds, ncol=3, byrow=TRUE)
## LDS - late dry season (July to Dec) ###LDS definition may change (Aug - Dec)
lds <- c(-1,0,0, 0,7,0, 7,12,1)
ldsmat <- matrix(lds, ncol=3, byrow=TRUE)
## Binary raster of burnt and not burnt
bin <- c(-1,0,0, 0,12,1)
binmat <- matrix(bin, ncol=3, byrow=TRUE)
##Unburnt
ub <- c(-1,0,1, 0,12,0)
ubmat <- matrix(ub, ncol=3, byrow=TRUE)
# Loop - creates stack of burn by month per year --------------------------
shps <- dir_ls("../../Source/NAFI_SHP", glob = "*.shp$")
yrs <- substr(sapply(str_split(names(shps), pattern = "/"), tail, 1), 1, 4)
fs_stack <- raster::stack()
for(i in seq_along(shps)){
shp <- shps[i]
f_mths <- st_read(shp) %>%
clean_names() %>% # handles mix upper/lower case
filter(month < 13) %>%
st_transform(3577) %>%
fasterize(tmpltrast, field = "month")
f_mths[is.na(f_mths[])] <- 0
f_mths <- mask(f_mths, kimbcoast_alb)
fs_stack <- raster::stack(fs_stack, f_mths)
}
names(fs_stack) <- yrs
# Apply classification matrices to the stack ------------------------------
# classified stacks
fs_eds_st <- reclassify(fs_stack, edsmat)
fs_lds_st <- reclassify(fs_stack, ldsmat)
fs_bin_st <- reclassify(fs_stack, binmat)
fs_ub_st <- reclassify(fs_stack, ubmat)
# writeRaster(fs_eds_st, filename = "./stacks/fs_eds_st.tif", datatype = 'INT1U',
# overwrite = TRUE)
# writeRaster(fs_lds_st, filename = "./stacks/fs_lds_st.tif", datatype = 'INT1U',
# overwrite = TRUE)
# writeRaster(fs_bin_st, filename = "./stacks/fs_bin_st.tif", datatype = 'INT1U',
# overwrite = TRUE)
# writeRaster(fs_ub_st, filename = "./stacks/fs_ub_st.tif", datatype = 'INT1U',
# overwrite = TRUE)
# Create burn year stack --------------------------------------------------
fs_yr_st <- raster::stack()
for(i in seq_along(names(fs_bin_st))){
yr <- as.numeric(substr(yrs[i], 3, 4)) + 2000
r_layer <-fs_bin_st[[i]]
r_layer[r_layer == 1] <- yr
fs_yr_st <- raster::stack(fs_yr_st, r_layer)
}
names(fs_yr_st) <- yrs
# writeRaster(fs_ub_st, filename = "./stacks/fs_yr_st.tif", datatype = 'INT2U',
# overwrite = TRUE)
# Create time since last burn layers -------------------------------------
l <- nlayers(fs_yr_st)
for(i in seq_along(names(fs_yr_st))){
if(i + 1 < l){
ind <- i + 1
mstack <- fs_yr_st[[1:ind]]
on1 <- paste0("./tslb/", gsub("fs", "tslb_", names(mstack)[ind]), ".tif")
mval <- cellStats(mstack[[ind]], 'max', na.rm = TRUE)
cat("doing...", names(mstack)[ind])
beginCluster()
ol <- clusterR(mstack, fun = calc,
args = list(fun = function(x){max(x, na.rm = TRUE)}))
endCluster()
tslb <- calc(ol, fun = function(x){ifelse(x == 0, ind, (mval - x))})# current yr - whatever is there
writeRaster(tslb, filename = on1, datatype = 'INT1U', overwrite = TRUE)
} else {
ind <- l
mstack <- fs_yr_st[[1:ind]]
on1 <- paste0("./tslb/", gsub("fs", "tslb_", names(mstack)[ind]), ".tif")
mval <- cellStats(mstack[[ind]], 'max', na.rm = TRUE)
cat("doing...", names(mstack)[ind])
beginCluster()
ol <- clusterR(mstack, fun = calc,
args = list(fun = function(x){max(x, na.rm = TRUE)}))
endCluster()
tslb <- calc(ol, fun = function(x){ifelse(x == 0, ind, (mval - x))})# current yr - whatever is there
writeRaster(tslb, filename = on1, datatype = 'INT1U', overwrite = TRUE)
}
}
# Create fire frequency since 2000 layers ---------------------------------
## function to output fire frequency since 2000 written to /fireFreq
## fs_bin_st is the binary stack
## s is the initial layer in the stack to calculate back from
fireFreq <- function(fs_bin_st, s = 5){
l <- nlayers(fs_bin_st)
ffval <- c(-1,0,0, 0,1,1, 1,2,2, 2,3,3, 3,4,4, 4,l,5)
ffmat <- matrix(ffval, ncol=3, byrow=TRUE)
for(i in seq_along(names(fs_bin_st))){
ind <- (s - 1) + i
if(ind < l){
mstack <- fs_bin_st[[1:ind]]
on <- paste0("./fireFreq/", gsub("fs", "ff", names(mstack)[ind]), ".tif")
ol <- calc(mstack, fun = sum, na.rm = TRUE)
olc <- reclassify(ol, ffmat)
writeRaster(olc, datatype = 'INT1U', filename = on, overwrite = TRUE)
} else {
mstack <- fs_bin_st[[1:l]]
on <- paste0("./fireFreq/", gsub("fs", "ff", names(mstack)[l]), ".tif")
ol <- calc(mstack, fun = sum, na.rm = TRUE)
olc <- reclassify(ol, ffmat)
writeRaster(olc, datatype = 'INT1U', filename = on, overwrite = TRUE)
}
}
}
fireFreq(fs_bin_st, s = 5)
# Distance to unburnt -----------------------------------------------------
## loop to make cost distance to unburnt (no size limit) for each year, written
## to /unburnt and makes a stack for later use
# unb_dist_stack <- stack()
for(i in seq_along(names(fs_ub_st))){
#first bit
cat("doing...", names(fs_ub_st)[i])
mlayer <- gridDistance(fs_ub_st[[i]], origin =1, omit = NA)
names(mlayer) <- names(fs_ub_st[[i]])
on <- paste0("./unburnt/", gsub("fs", "ubdist_", names(fs_ub_st)[i]), ".tif")
# unb_dist_stack <- stack(unb_dist_stack, mlayer)
writeRaster(mlayer, filename = on, datatype = 'INT2U', overwrite = TRUE)
#second clumpy bit
c <- clump(fs_ub_st[[i]], directions = 4, gaps = FALSE)
c[is.na(c)] <- 0 #burnt areas to 0 (can't be NA)
m_c <- mask(c, kimbcoast_alb) #remask water to NA leaving burnt as 0
f <- freq(c, digits = 0, useNA = 'no', progress = 'text')# on orig clump as don't want burnt to be counted as clump in stats
df <- as_tibble(f) %>%
dplyr::mutate(ha = count * 6.25,
gr20ha = ifelse(ha > 20, 1, 0),
gr20ha = ifelse(value == 0, 0, gr20ha)) %>%
dplyr::select(value, gr20ha)
mat <- data.matrix(df, rownames.force = NA)
g20 <- reclassify(m_c, mat)
g20_out <- gridDistance(g20, origin = 1, omit = NA)
on <- paste0("./unburnt/", gsub("fs", "ubdist20_", names(fs_ub_st)[i]), ".tif")
writeRaster(g20_out, filename = on, datatype = 'INT2U', overwrite = TRUE)
}
# Years since last burn ---------------------------------------------------
## function to output ageclass raster layers for various stat metrics written to
## /ageClass
## fs_yr_st is the year burnt stack
## is parallelised but is slow to run (many hours)
## Parked at this stage - created tslb layer which can be classified to whatever
## year ranges are desired
# ageClass <- function(fs_yr_st){
# l <- nlayers(fs_yr_st)
# ar <- c(-1,0,0, 0,1,1, 1,2,2, 2,3,3, 3,4,4, 4,l,5)
# armat <- matrix(ar, ncol=3, byrow=TRUE)
# y3 <- c(-1,2,0, 2,l,1)
# y3mat <- matrix(y3, ncol=3, byrow=TRUE)
# y5 <- c(-1,4,0, 4,l,1)
# y5mat <- matrix(y5, ncol=3, byrow=TRUE)
#
# for(i in seq_along(names(fs_yr_st))){
# if(i + 1 < l){
# ind <- i + 1
# mstack <- fs_yr_st[[1:ind]]
# on1 <- paste0("./ageClass/", gsub("fs", "ac_annual_", names(mstack)[ind]), ".tif")
# on3 <- gsub("ac_annual_", "ac3_", on1)
# on5 <- gsub("ac_annual_", "ac5_", on1)
# mval <- cellStats(mstack[[ind]], 'max', na.rm = TRUE)
# cat("doing...", names(mstack)[ind])
# beginCluster()
# ol <- clusterR(mstack, fun = calc,
# args = list(fun = function(x){max(x, na.rm = TRUE)}))
# endCluster()
# ysb <- calc(ol, fun = function(x){ifelse(x == 0, ind, (mval - x))})# current yr - whatever is there
# ysb_1 <- reclassify(ysb, armat)
# writeRaster(ysb_1, filename = on1, datatype = 'INT1U', overwrite = TRUE)
# ysb_3 <- reclassify(ysb, y3mat)
# writeRaster(ysb_3, filename = on3, datatype = 'INT1U', overwrite = TRUE)
# ysb_5 <- reclassify(ysb, y5mat)
# writeRaster(ysb_5, filename = on5, datatype = 'INT1U', overwrite = TRUE)
# } else {
# ind <- l
# mstack <- fs_yr_st[[1:ind]]
# on1 <- paste0("./ageClass/", gsub("fs", "ac_annual_", names(mstack)[ind]), ".tif")
# on3 <- gsub("ac_annual_", "ac3_", on1)
# on5 <- gsub("ac_annual_", "ac5_", on1)
# mval <- cellStats(mstack[[ind]], 'max', na.rm = TRUE)
# cat("doing...", names(mstack)[ind])
# beginCluster()
# ol <- clusterR(mstack, fun = calc,
# args = list(fun = function(x){max(x, na.rm = TRUE)}))
# endCluster()
# ysb <- calc(ol, fun = function(x){ifelse(x == 0, ind, (mval - x))})
# ysb_1 <- reclassify(ysb, armat)
# writeRaster(ysb_1, filename = on1, datatype = 'INT1U', overwrite = TRUE)
# ysb_3 <- reclassify(ysb, y3mat)
# writeRaster(ysb_3, filename = on3, datatype = 'INT1U', overwrite = TRUE)
# ysb_5 <- reclassify(ysb, y5mat)
# writeRaster(ysb_5, filename = on5, datatype = 'INT1U', overwrite = TRUE)
# }
# }
# }
# ageClass(fs_yr_st)
# Distance last burn by age class -----------------------------------------
## On hold not naming correctly at present
# yslb <- dir_ls("./ageClass")
# y3s <- yslb[str_detect(yslb, pattern = ".tif$") & str_detect(yslb, pattern = "ac3_")]
# y5s <- yslb[str_detect(yslb, pattern = ".tif$") & str_detect(yslb, pattern = "ac5_")]
# for(i in seq_along(names(fs_yr_st))){
# on3 <- paste0("./ageClass_dist/", gsub("fs", "ac3_dist_", names(fs_yr_st)[i]), ".tif")
# on5 <- paste0("./ageClass_dist/", gsub("fs", "ac5_dist_", names(fs_yr_st)[i]), ".tif")
#
# g3 <- gridDistance(raster(y3s[i]), origin = 1, omit = NA)
# g5 <- gridDistance(raster(y5s[i]), origin = 1, omit = NA)
# writeRaster(g3, filename = on3, datatype = 'INT2U', overwrite = TRUE)
# writeRaster(g5, filename = on5, datatype = 'INT2U', overwrite = TRUE)
#
# }