-
Notifications
You must be signed in to change notification settings - Fork 1
/
Pasture preprocessing 2.R
163 lines (132 loc) · 6.92 KB
/
Pasture preprocessing 2.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
library(raster)
library(tidyverse)
library(sp)
bigdata_repo <- "GIS data repository"
projdata_repo <- "DEFRA Clean Growth Project/pH Optimisation/Extension for publication"
# soil data stack
Soil <- stack(find_onedrive(dir = bigdata_repo, path = "SoilGrids 5km/Sand content/Fixed/SNDPPT_M_sl4_5km_ll.tif"),
find_onedrive(dir = bigdata_repo, path = "SoilGrids 5km/Clay content/Fixed/CLYPPT_M_sl4_5km_ll.tif"),
find_onedrive(dir = bigdata_repo, path = "SoilGrids 5km/Depth to bedrock/Fixed/BDRICM_M_5km_ll.tif"),
find_onedrive(dir = bigdata_repo, path = "SoilGrids 5km/OC tonnes per ha/Fixed/OCSTHA_M_30cm_5km_ll.tif"),
find_onedrive(dir = bigdata_repo, path = "SoilGrids 5km/Bulk density/Fixed/BLDFIE_M_sl4_5km_ll.tif"),
find_onedrive(dir = bigdata_repo, path = "soilGrids 5km/Soil pH/Fixed/PHIHOX_M_sl4_5km_ll.tif"))
# precipitation stack
Precip <- raster::stack()
for(i in 4:9){
path <- find_onedrive(dir = bigdata_repo, path = paste("WorldClim data/Precipitation (mm) 30-arc-secs/wc2.0_30s_prec_", formatC(i, width=2, flag="0"), ".tif", sep=""))
x <- raster(path)
Precip <- addLayer(Precip, x)
}
rm(x, i , path)
# uk shapefile
Shp_UK <- find_onedrive(dir = bigdata_repo, path = "DA shapefile/GBR_adm_shp/GBR_adm1.shp") %>% shapefile()
# crop and mask to UK
Soil_UK <- Soil %>% crop(Shp_UK)# %>% mask(Shp_UK)
Precip_UK <- Precip %>% crop(Shp_UK)# %>% mask(Shp_UK)
# resample precipitation to correct extent + resolution
Precip_UK <- Precip_UK %>% resample(Soil_UK[[1]])
# create master stack and convert to df
Master <- stack(Soil_UK, Precip_UK)
Dat_main <- as.data.frame(Master, xy = T) %>%
as_tibble() %>%
drop_na()
# name properly, sum up precipitation and calculate clay fraction
precip_temp <- Dat_main %>% select(wc2.0_30s_prec_04:wc2.0_30s_prec_09)
Dat_main <- Dat_main %>%
rename(Sand = SNDPPT_M_sl4_5km_ll,
Clay = CLYPPT_M_sl4_5km_ll,
Bedrock = BDRICM_M_5km_ll,
OC = OCSTHA_M_30cm_5km_ll,
BD = BLDFIE_M_sl4_5km_ll,
pH = PHIHOX_M_sl4_5km_ll) %>%
select(-(wc2.0_30s_prec_04:wc2.0_30s_prec_09)) %>%
mutate(Precip_mm = rowSums(precip_temp),
Silt = 100 - (Sand + Clay),
pH = pH / 10)
rm(precip_temp)
# assign soil type based on texture data (details here https://cran.r-project.org/web/packages/soiltexture/vignettes/soiltexture_vignette.pdf)
library(soiltexture)
Dat_soil <- Dat_main %>%
dplyr::select(SAND = Sand, SILT = Silt, CLAY = Clay, OC) %>%
mutate(OC = OC / 10) %>%
as.data.frame()
Dat_main <- Dat_main %>%
mutate(Soil_type = TT.points.in.classes(tri.data = Dat_soil,
class.sys = "UK.SSEW.TT",
PiC.type = "t"))
# read in soil category translations (manual) and join to main data
# borrowing from another project's directory here
Soil_cats <- read_csv(find_onedrive(dir = "AgRE Calc PLC/Soil C methodology/Model update/Data preprocessing/Defra RB209 grass model",
path = "Defra soil type translations.csv"))
Dat_main <- Dat_main %>% left_join(Soil_cats %>% rename(Soil_type = cat, RB209_soiltype = trans), by = "Soil_type")
# function to calculate OM fraction based on C (t/ha) and BD (kg / m2)
OM_frac <- function(BD_kg_m2, C_t_ha){
C_kg_m2 <- C_t_ha * 10^3 * 10^-4 * 1 / 0.3
Cfrac <- C_kg_m2 / BD_kg_m2
OMfrac <- Cfrac / 0.58
return(OMfrac)
}
# calculate other variables not based on SSEW soil types and adjust RB209 definitions accordingly
Dat_main <- Dat_main %>%
mutate(RB209_soiltype = ifelse(Bedrock <= 40, "Shallow soils", RB209_soiltype), # shallow soils <= 40cm
OMfrac = OM_frac(BD, OC),
RB209_soiltype = ifelse(OMfrac >= 0.1, "Organic soils", RB209_soiltype), # organic soils
RB209_soiltype = ifelse(OMfrac >= 0.2, "Peat soils", RB209_soiltype)) # peat soils
# water availability class from Defra RB209
Water_availability <- Dat_main %>%
distinct(RB209_soiltype) %>%
arrange(RB209_soiltype) %>%
mutate(SAW = c(2, 3, 1, 2, 3)) # soil available water class — 1 = low, 2 = medium, 3 = high
Dat_main <- Dat_main %>% left_join(Water_availability, by = "RB209_soiltype")
# rainfall classes from Defra RB209 and final pasture growth class
# addition method is something I've inferred — works if you parameterise SAW from 1-3 and rain class from 0-2 or vice versa
# rain classes are <300, 300-400, >400
breaks <- c(min(Dat_main$Precip_mm) - 1, 300, 400, max(Dat_main$Precip_mm) + 1)
Dat_main <- Dat_main %>%
mutate(RC = Precip_mm %>%
cut(breaks = breaks, labels = F) %>%
as.numeric(),
GGC = SAW + RC - 1) # 1 = V poor, 2 = Poor, 3 = Av, 4 = Good, 5 = V good
Dat_main$GGC %>% qplot()
# read in raw data for Defra's grass growth model
# we're borrowing from another project here — so pulling project files directly from that project's repo
Dat_GGM <- find_onedrive(dir = "AgRE Calc PLC/Soil C methodology/Model update/Data preprocessing/Defra RB209 grass model",
path = "Defra RB209 final models.csv") %>%
read_csv()
# gather data and fit models
Dat_GGM <- Dat_GGM %>%
gather(Vpoor:Vgood, key = "GGC", value = "Yield")
GGC_models = list(Vpoor = loess(Yield ~ N_rate, data = Dat_GGM %>% filter(GGC == "Vpoor")),
Poor = loess(Yield ~ N_rate, data = Dat_GGM %>% filter(GGC == "Poor")),
Av = loess(Yield ~ N_rate, data = Dat_GGM %>% filter(GGC == "Av")),
Good = loess(Yield ~ N_rate, data = Dat_GGM %>% filter(GGC == "Good")),
Vgood = loess(Yield ~ N_rate, data = Dat_GGM %>% filter(GGC == "Vgood")))
# make pasture yield predictions
# unable to get spatially specific N estimates as yet so using BFSP N application estimate for pasture
BSFP_Nrate <- 57
Pas_yield <- function(GGC){
predict(GGC_models[[GGC]], newdata = tibble(N_rate = BSFP_Nrate))
}
Dat_main <- Dat_main %>%
mutate(Yield = map_dbl(GGC, Pas_yield))
# adjust pasture yield predictions for pH
Dat_yieldres_pasture <- read_csv(find_onedrive(dir = projdata_repo, path = "Fornara yield response.csv"))
Yield_adjust_fac <- function(pH){
pH_gap <- ifelse(pH < 6, 6 - pH, 0)
Yield_effect <- Dat_yieldres_pasture$mean[1] - 1
Yield_adjust_fac <- 1 / (1 + Yield_effect * pH_gap)
return(Yield_adjust_fac)
}
Dat_main <- Dat_main %>%
mutate(Yield = Yield * Yield_adjust_fac(pH))
# write to raster
Pasture_yield <- rasterFromXYZ(Dat_main %>% select(x, y, Yield),
crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
plot(Pasture_yield)
Pasture_yield
# read in area raster, extend, and mask
Pasture_area <- find_onedrive(dir = bigdata_repo, path = "Created rasters/UK-pasture-area-10km-CLC-based-WGS84.tif") %>% raster()
Pasture_area
Pasture_yield <- extend(Pasture_yield, Pasture_area)
Pasture_yield[is.na(Pasture_area)] <- NA
writeRaster(Pasture_yield, find_onedrive(dir = bigdata_repo, path = "Created rasters/UK-pasture-yield-RB209-10km.tif"), overwrite = T)