-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpostProcessing2_weightedHabitat.R
181 lines (148 loc) · 7.81 KB
/
postProcessing2_weightedHabitat.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
#purpose is two-fold:
# convert all focal habitat rasters into single weighted habitat layer for each year
# calculate the difference
# generate a habitat class layer by assigning 1:8 to binary habitats
focalHabitatDir <- "outputs/focalHabitat"
weightedDir <- "outputs/weightedHabitat"
compositeDir <- "outputs/compositeHabitat"
checkPath(weightedDir, create = TRUE)
tiles <- paste0("tile", 1:c(nx*ny))
#Weighted habitat ####
makeWeightedHabitat <- function(tileList, year, outputPath) {
tileList <- tileList[grep(tileList, pattern = year)] #1985 or 2020
if (length(tileList) == 0) {stop("incorrect year")}
#as a float, each individual habitat tile is 20 GB
#this function writes too many temp files so clean up must occur inside lapply loop
#put in function for terra so filename can be provided, allowing for easy deletion
multi <- function(x, by){x * by}
tempFile1 <- tempfile(fileext = ".tif")
tempFile2 <- tempfile(fileext = ".tif")
tempFile3 <- tempfile(fileext = ".tif")
tempFile4 <- tempfile(fileext = ".tif")
fire <- app(rast(tileList[grep("naturalDisturbance", tileList)]),
fun = multi, by = 0.06, filename = tempFile1)
harvestY <- app(rast(tileList[grep("harvest_0to5", tileList)]),
fun = multi, by = 0.04, filename = tempFile2)
harvestO <- app(rast(tileList[grep("harvest_6to20", tileList)]),
fun = multi, by = 0.04, filename = tempFile3)
weightedHabitat <- sum(harvestY, harvestO, fire, na.rm = TRUE, filename = tempFile4)
rm(harvestY, harvestO, fire)
gc()
coniferY <- app(x = rast(tileList[grep("youngConifer", tileList)]),
fun = multi, by = 0.19, filename = tempFile2, overwrite = TRUE)
coniferO <- app(rast(tileList[grep("matureConifer", tileList)]),
fun = multi, by = 0.25, filename = tempFile1, overwrite = TRUE)
weightedHabitat <- sum(coniferY, coniferO, weightedHabitat, na.rm = TRUE,
filename = tempFile3, overwrite = TRUE)
rm(coniferY, coniferO)
gc()
woodlands <- app(x = rast(tileList[grep("openWoodland", tileList)]),
fun = multi, by = 0.22, filename = tempFile1, overwrite = TRUE)
wetlands <- app(rast(tileList[grep("wetland", tileList)]),
fun = multi, by = 0.14, filename = tempFile2, overwrite = TRUE)
weightedHabitat <- sum(weightedHabitat, woodlands, wetlands, filename = tempFile4,
na.rm = TRUE, overwrite = TRUE)
rm(woodlands, wetlands)
gc()
regen <- app(rast(tileList[grep("regenerating", tileList)]),
fun = multi, by = 0.06, filename = tempFile1, overwrite = TRUE)
#prepare final file
tileNum <- stringr::str_extract(tileList[1], pattern = "tile[0-9]+")
outFile <- file.path(outputPath, paste0("weightedHabitat_", year, "_", tileNum, ".tif"))
#TODO: I believe terra needs app(weightedHabitat, sum) to pass args like overwrite correctly
weightedHabitat <- sum(weightedHabitat, regen, na.rm = TRUE, filename = outFile,
overwrite = TRUE, wopt = list(datatype = "INT2U"))
#clean up
#terra now has a tmpFiles function that will do this for you
unlink(x = c(tempFile1, tempFile2, tempFile3, tempFile4))
gc()
message(tileNum, " complete")
}
habitatRasters <- lapply(tiles, list.files, path = focalHabitatDir, full.names = TRUE)
#this could be stored as INT1U - as the theoretical max habitat is 240 (1000 * 0.24 weight)
#this fills up the temp drive so don't run all 8 tiles and 16 years
lapply(habitatRasters, FUN = makeWeightedHabitat, year = 2020, outputPath = weightedDir)
lapply(habitatRasters, FUN = makeWeightedHabitat, year = 1985, outputPath = weightedDir)
#####Composite Habitat#####
#fire = 1, young/old harvest = 2 and 3, young/mature conifer = 4 and 5,
#open woodland 6, wetland 7, regenerating forest = 8
makeCompositeHabitat <- function(tileList, year, outputPath) {
checkPath(outputPath, create = TRUE)
tileList <- tileList[grep(tileList, pattern = year)] #1985 or 2020
if (length(tileList) == 0) {stop("incorrect year")}
#as a float, each individual habitat tile is 20 GB
tempFile1 <- tempfile(fileext = ".tif")
tempFile2 <- tempfile(fileext = ".tif")
tempFile3 <- tempfile(fileext = ".tif")
tempFile4 <- tempfile(fileext = ".tif")
multi <- function(x, by){x * by}
fire <- rast(tileList[grep("naturalDisturbance", tileList)]) #this one is 1 - so it doesn't need to be multiplied
harvestY <- app(x = rast(tileList[grep("harvest_0to5", tileList)]),
fun = multi, by = 2, filename = tempFile1)
harvestO <- app(x = rast(tileList[grep("harvest_6to20", tileList)]),
fun = multi, by = 3, filename = tempFile2)
compositeHabitat <- sum(harvestY, harvestO, fire, filename = tempFile3, na.rm = TRUE)
rm(harvestY, harvestO, fire)
gc()
coniferY <- app(rast(tileList[grep("youngConifer", tileList)]),
fun = multi, by = 4, filename = tempFile1, overwrite = TRUE)
coniferO <- app(rast(tileList[grep("matureConifer", tileList)]),
fun = multi, by = 5, filename = tempFile2, overwrite = TRUE)
compositeHabitat <- sum(coniferY, coniferO, compositeHabitat, na.rm = TRUE,
filename = tempFile4)
rm(coniferY, coniferO)
gc()
woodlands <- app(rast(tileList[grep("openWoodland", tileList)]),
fun = multi, by = 6, filename = tempFile1, overwrite = TRUE)
wetlands <- app(rast(tileList[grep("wetland", tileList)]),
fun = multi, by = 7, filename = tempFile2, overwrite = TRUE)
compositeHabitat <- sum(compositeHabitat, woodlands, wetlands, na.rm = TRUE,
filename = tempFile3, overwrite = TRUE)
rm(woodlands, wetlands)
gc()
regen <- app(rast(tileList[grep("regenerating", tileList)]),
fun = multi, by = 8, filename = tempFile1, overwrite = TRUE)
tileNum <- stringr::str_extract(tileList[1], pattern = "tile[0-9]+")
outFile <- file.path(outputPath, paste0("compositeHabitat_", year, "_", tileNum, ".tif"))
compositeHabitat <- sum(compositeHabitat, regen, na.rm = TRUE, filename = outFile, overwrite = TRUE)
rm(regen)
gc()
if (max(terra::minmax(compositeHabitat)) > 8) {
browser()
}
unlink(c(tempFile1, tempFile2, tempFile3, tempFile4))
message(tileNum, " complete")
}
habitatClasses <- lapply(tiles, list.files, path = "outputs/raw/", full.names = TRUE)
#some memory leakage happens with this function. Not sure why. Don't recommend running all at once
lapply(habitatClasses, makeCompositeHabitat, year = 2020, outputPath = compositeDir)
lapply(habitatClasses, makeCompositeHabitat, year = 1985, outputPath = compositeDir)
gc()
#upload files
#TODO: add the requisite checks for whether to upload, instead of this if block protection
if (FALSE) {
toZip <- list.files("outputs/weightedHabitat", full.names = TRUE)
utils::zip(zipfile = "outputs/weightedHabitat.zip",
files = toZip,
flags = "-j")
thePath <- as_dribble("PFC/Yan/Caribou Hindcast Results V2/weighted habitat")
drive_put("outputs/weightedHabitat.zip", path = thePath)
}
if (FALSE){
# #the composite tiles are 8 GB each, so upload tiles separately)
toZip <- list.files("outputs/focalHabitat1000", full.names = TRUE)
utils::zip(zipfile = "outputs/focalHabitat1000.zip",
files = toZip,
flags = "-j")
thePath <- googledrive::as_dribble("PFC/Yan/Caribou Hindcast Results V2/focal habitat layers X 1000")
drive_put("outputs/focalHabitat1000.zip", path = thePath)
}
if (FALSE){
# #the composite tiles are 8 GB each, so upload tiles separately)
toZip <- list.files("outputs/compositeHabitat", full.names = TRUE)
utils::zip(zipfile = "outputs/compositeHabitat.zip",
files = toZip,
flags = "-j")
thePath <- googledrive::as_dribble("PFC/Yan/Caribou Hindcast Results V2/focal habitat layers X 1000")
drive_put("outputs/focalHabitat1000.zip", path = thePath)
}