Section 16 Landuse Scenarios

RUN_THIS_CHAPTER = FALSE

WARNING: This code was rushed and abandoned before it could be properly added to the GitBook. It needs a look-over

First we need to fix the land-use names to be FarmRInput Compatible. The reason this is broken is because we changed our drainage classes after running the buildR, and after messing with our land use lum file.

require(sf)
require(dplyr)
require(mapview)

Update the names based on drainage and LUM

lu_v2 <- sf::read_sf("model_data/input/land/CS10_LU_V2.shp")

get_lum <- which(grepl(x = lu_v2$type, pattern = "a_"))
get_drn <- which((lu_v2$drainage %>% is.na()) == FALSE)

lu_v2$namefixed <- lu_v2$type
lu_v2$namefixed[get_drn] <- paste0(lu_v2$type[get_drn], "_drn")
mapview::mapview(lu_v2)

Now we update our crop rotation with the working names

crop_rotation <- sf::read_sf("model_data/input/crop/land_with_cr.shp")
crop_rotation$lu <- lu_v2$namefixed
crop_rotation <- crop_rotation
write_sf(crop_rotation, "model_data/input/crop/land_with_cr_fixed.shp")

And prepare our FarmR input.

16.1 Scenarios

16.1.1 Status Quo

require(sf)
require(tidyverse)
require(lubridate)
require(reshape2)
require(dplyr)
require(data.table)

rm(list = ls())
source('model_data/code/functions_write_SWATfarmR_input.R')

lu_shp <- 'model_data/input/crop/land_with_cr_fixed.shp'
lu <- st_drop_geometry(read_sf(lu_shp))

mgt_csv <- 'model_data/input/management/mgt_crops_CS10.csv'
mgt_crop <- read.csv(mgt_csv)
lu_generic_csv <- 'model_data/input/management/farmR_generic_CS10_ext.csv' # generic land use management .csv table
mgt_generic <- read.csv(lu_generic_csv)


### TODO  DEFINE THIS
start_y <- 2013 # starting year (consider at least 3 years for warm-up!)
end_y <- 2020 # ending year

hru_crops <- 'a_'
m_yr_sch_existing <- 'n'
crop_myr <- 'past' # name of your farmland grass
max_yr <- 5
crop_s <- c('sgbt','csil','barl')
additional_h_yr_sch_existing <- 'n' # 'y' (yes), 'n' (no)
scen_path = "model_data/input/management/management_scenarios/status_quo/"

check_skip <- check_skip_position(path = scen_path)
check_date_conflicts1()
rota_schedules <- build_rotation_schedules()
check_date_conflicts2()
rota_schedules <- solve_date_conflicts()
check_date_conflicts2()
write_farmR_input(path = paste0(scen_path, "farmR_input_status_quo.csv"))

16.1.1.1 Running farmR

require(SWATfarmR)
dir_swat <- paste0(scen_path, "swat_setup")
unlink(dir_swat, recursive = T)
dir.create(dir_swat, showWarnings = F)
swat_files <- list.files("model_data/cs10_setup/MetNord_MN2/", full.names = T)
file.copy(from = swat_files, to = dir_swat, overwrite = T)

SWATfarmR::load_farmr(file = paste0(dir_swat,"/cs10.farm"))
cs10$reset_files()
unlink(paste0(dir_swat,"/cs10.farm"))
unlink(paste0(dir_swat, "/cs10.mgts"))

fip <- "model_data/input/management/management_scenarios/status_quo/farmR_input_status_quo.csv"
SWATfarmR::new_farmr(project_path = dir_swat, project_name = "status_quo")
status_quo$read_management(file = fip, discard_schedule = T)
pcp <- status_quo$.data$variables$pcp

# using temp API
# TODO add real API
  
# Extract the hydrologic soil group values for all HRUs
hsg <- select(status_quo$.data$meta$hru_attributes, hru, hyd_grp)

  

# Calculate api values for the hsg classes A to D
api_A <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 1)
api_B <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 0.8)
api_C <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 0.7)
api_D <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 0.5)
  
  # Bind the data together into one api table and name them with the hsgs
  api <- bind_cols(api_A, api_B, api_C, api_D)
  names(api) <- c('api_A', 'api_B', 'api_C', 'api_D')
  
  # To add the variable to the farmR you have to tell it which variables are
  # assigned to which HRUs
  hru_asgn <- mutate(hsg, api = paste0('api_', hyd_grp)) %>% select(hru, api)

  # Add the variable api to the farmR project
  status_quo$add_variable(data = api, name = 'api', assign_unit = hru_asgn, overwrite = T)
  
  # Schedule operations
  status_quo$schedule_operations(start_year = 2013,
                         end_year = 2020,
                         n_schedule = 2, replace = "all")
  # Write operations
  status_quo$write_operations(start_year = 2013, end_year = 2020)

16.1.2 No Autumn Tillage

For the no-autumn tillage we will simply replace wwht (winter wheat) with swht (summer wheat).

# WARNING
# Micha's code is using global variables, so be careful
rm(list = ls())
source('model_data/code/functions_write_SWATfarmR_input.R')


lu_shp <- 'model_data/input/crop/land_with_cr_fixed.shp'
lu <- st_drop_geometry(read_sf(lu_shp))

lu <- data.frame(lapply(lu, function(x) {
                  gsub("wwht", "swht", x)
}))

lu<-lu %>% tibble()

mgt_csv <- 'model_data/input/management/mgt_crops_CS10.csv'
mgt_crop <- read.csv(mgt_csv)

# remove wwht from mgt_crop
#mgt_crop <- mgt_crop[-which(mgt_crop$crop_mgt == "wwht"),]

lu_generic_csv <- 'model_data/input/management/farmR_generic_CS10_ext.csv' # generic land use management .csv table
mgt_generic <- read.csv(lu_generic_csv)


### TODO  DEFINE THIS
start_y <- 2013 # starting year (consider at least 3 years for warm-up!)
end_y <- 2020 # ending year

hru_crops <- 'a_'
m_yr_sch_existing <- 'n'
crop_myr <- 'past' # name of your farmland grass
max_yr <- 5
crop_s <- c('sgbt','csil','barl')
additional_h_yr_sch_existing <- 'n' # 'y' (yes), 'n' (no)


scen_path = "model_data/input/management/management_scenarios/no_at/"
check_skip <- check_skip_position(path = scen_path)
check_date_conflicts1()
rota_schedules <- build_rotation_schedules()
check_date_conflicts2()
rota_schedules <- solve_date_conflicts()
check_date_conflicts2()
write_farmR_input(path = paste0(scen_path, "farmR_input_no_at.csv"))

16.1.2.1 Running FamR

require(SWATfarmR)
scen_name <- "no_at"
scen_path = paste0("model_data/input/management/management_scenarios/", scen_name, "/")
dir_swat <- paste0(scen_path, "swat_setup")
fip <- paste0("model_data/input/management/management_scenarios/",scen_name, "/farmR_input_", scen_name, ".csv")
dir_swat <- paste0(scen_path, "swat_setup")
unlink(dir_swat, recursive = T)
dir.create(dir_swat, showWarnings = F)
swat_files <- list.files("model_data/cs10_setup/MetNord_MN2/", full.names = T)
file.copy(from = swat_files, to = dir_swat, overwrite = F)

SWATfarmR::load_farmr(file = paste0(dir_swat,"/cs10.farm"))
cs10$reset_files()
unlink(paste0(dir_swat,"/cs10.farm"))
unlink(paste0(dir_swat, "/cs10.mgts"))


SWATfarmR::new_farmr(project_path = dir_swat, project_name = "farmscen")
farmscen$read_management(file = fip, discard_schedule = T)
pcp <- farmscen$.data$variables$pcp

# using temp API
# TODO add real API
  
# Extract the hydrologic soil group values for all HRUs
hsg <- select(farmscen$.data$meta$hru_attributes, hru, hyd_grp)

  

# Calculate api values for the hsg classes A to D
api_A <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 1)
api_B <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 0.8)
api_C <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 0.7)
api_D <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 0.5)
  
  # Bind the data together into one api table and name them with the hsgs
  api <- bind_cols(api_A, api_B, api_C, api_D)
  names(api) <- c('api_A', 'api_B', 'api_C', 'api_D')
  
  # To add the variable to the farmR you have to tell it which variables are
  # assigned to which HRUs
  hru_asgn <- mutate(hsg, api = paste0('api_', hyd_grp)) %>% select(hru, api)

  # Add the variable api to the farmR project
  farmscen$add_variable(data = api, name = 'api', assign_unit = hru_asgn, overwrite = T)
  
  # Schedule operations
  farmscen$schedule_operations(start_year = 2013,
                         end_year = 2020,
                         n_schedule = 2, replace = "all")
  # Write operations
  farmscen$write_operations(start_year = 2013, end_year = 2020)

16.1.2.2 Stubbe V2

# WARNING
# Micha's code is using global variables, so be careful
rm(list = ls())
source('model_data/code/functions_write_SWATfarmR_input.R')
mgt_scen_name <- "stubble_v2"


lu_shp <- 'model_data/input/crop/land_with_cr_fixed.shp'
lu <- st_drop_geometry(read_sf(lu_shp))

lu <- data.frame(lapply(lu, function(x) {
                  gsub("wwht", "swht", x)
}))

lu<-lu %>% tibble()

mgt_csv <- 'model_data/input/management/mgt_crops_CS10_Stubble_v2.csv'
mgt_crop <- read.csv(mgt_csv)

# remove wwht from mgt_crop
#mgt_crop <- mgt_crop[-which(mgt_crop$crop_mgt == "wwht"),]

lu_generic_csv <- 'model_data/input/management/farmR_generic_CS10_ext.csv' # generic land use management .csv table
mgt_generic <- read.csv(lu_generic_csv)


### TODO  DEFINE THIS
start_y <- 2013 # starting year (consider at least 3 years for warm-up!)
end_y <- 2020 # ending year

hru_crops <- 'a_'
m_yr_sch_existing <- 'n'
crop_myr <- 'past' # name of your farmland grass
max_yr <- 5
crop_s <- c('sgbt','csil','barl')
additional_h_yr_sch_existing <- 'n' # 'y' (yes), 'n' (no)


scen_path = paste0("model_data/input/management/management_scenarios/", mgt_scen_name)
dir.create(scen_path)
check_skip <- check_skip_position(path = scen_path)
check_date_conflicts1()
rota_schedules <- build_rotation_schedules()
check_date_conflicts2()
rota_schedules <- solve_date_conflicts()
check_date_conflicts2()
write_farmR_input(path = paste0(scen_path, "/farmR_input_", mgt_scen_name, ".csv"))


### Running farmR
require(SWATfarmR)
scen_name <- mgt_scen_name
scen_path = paste0("model_data/input/management/management_scenarios/", scen_name, "/")
dir_swat <- paste0(scen_path, "swat_setup")
fip <- paste0("model_data/input/management/management_scenarios/",scen_name, "/farmR_input_", scen_name, ".csv")
dir_swat <- paste0(scen_path, "swat_setup")
unlink(dir_swat, recursive = T)
dir.create(dir_swat, showWarnings = F)
swat_files <- list.files("model_data/cs10_setup/MetNord_MN2/", full.names = T)
file.copy(from = swat_files, to = dir_swat, overwrite = F)

SWATfarmR::load_farmr(file = paste0(dir_swat,"/cs10.farm"))
cs10$reset_files()
unlink(paste0(dir_swat,"/cs10.farm"))
unlink(paste0(dir_swat, "/cs10.mgts"))


SWATfarmR::new_farmr(project_path = dir_swat, project_name = "farmscen")
farmscen$read_management(file = fip, discard_schedule = T)
pcp <- farmscen$.data$variables$pcp

# using temp API
# TODO add real API
  
# Extract the hydrologic soil group values for all HRUs
hsg <- select(farmscen$.data$meta$hru_attributes, hru, hyd_grp)

  

# Calculate api values for the hsg classes A to D
api_A <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 1)
api_B <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 0.8)
api_C <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 0.7)
api_D <- variable_decay(variable = pcp, n_steps = -5, decay_rate = 0.5)
  
  # Bind the data together into one api table and name them with the hsgs
  api <- bind_cols(api_A, api_B, api_C, api_D)
  names(api) <- c('api_A', 'api_B', 'api_C', 'api_D')
  
  # To add the variable to the farmR you have to tell it which variables are
  # assigned to which HRUs
  hru_asgn <- mutate(hsg, api = paste0('api_', hyd_grp)) %>% select(hru, api)

  # Add the variable api to the farmR project
  farmscen$add_variable(data = api, name = 'api', assign_unit = hru_asgn, overwrite = T)
  
  # Schedule operations
  farmscen$schedule_operations(start_year = 2013,
                         end_year = 2020,
                         n_schedule = 2, replace = "all")
  # Write operations
  farmscen$write_operations(start_year = 2013, end_year = 2020)