Section 21 Measure Classification

RUN_THIS_CHAPTER = TRUE
if(RUN_THIS_CHAPTER){print(paste("Chapter last run:", Sys.Date()))}

## [1] "Chapter last run: 2024-11-15"

21.1 Introduction

Our measures are too loosely grouped, making it difficult for our optimization. We are aiming to have around 500 measures, currently we have more.

We will need our HRU ids spatially, as well as our Measures, erosion classes, and property register.

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

## [1] "sf 1.0.18"

## [1] "dplyr 1.1.4"

## [1] "DT 0.33"

Get our maps:

# For measure IDs and location:
lu_map <- read_sf("model_data/input/land/CS10_LU.shp")
# For HRU IDs:
hru_map <- read_sf("model_data/cs10_setup/optain-cs10/data/vector/hru.shp")
# For Erosion Risk Levels:
erosion_map <- read_sf("model_data/input/soil/erosionriskclasses_kraakstad.shp")
# For Farm ownership IDs:
farm_map <- read_sf("model_data/input/property/matrikkel.shp")
farm_map <- st_transform(farm_map, st_crs(lu_map))

Filter data

lu_filter <- lu_map %>% select("buffer_6m_", "gully", "wetland", "dam", "type", "geometry")
hru_filter <-  hru_map  %>% select("name", "type")
erosion_filter <- erosion_map %>% select("A_HPKL", "geometry")
farm_filter <- farm_map %>% select("GARDSNUMME", "geometry")
farm_dissolved <- farm_filter %>%
  group_by(GARDSNUMME) %>%
  summarise()
erosion_dissolved <- erosion_filter %>%
  group_by(A_HPKL) %>%
  summarise()

Map Farm ID

farm_dissolved$GARDSNUMME<-as.factor(farm_dissolved$GARDSNUMME)
farmcolors <- farm_dissolved$GARDSNUMME %>% unique() %>% length()
# samples randomizes legend
myfill <- hcl.colors(farmcolors, palette = "Dark 3") %>% sample()
mapview(farm_dissolved, zcol = "GARDSNUMME", legend = FALSE, col.regions = myfill)

Map Erosion Classes

myfill <- hcl.colors(4, palette = "Fall")
erosion_dissolved$A_HPKL <- as.factor(erosion_dissolved$A_HPKL)
mapview(erosion_dissolved, zcol = "A_HPKL", col.regions = myfill)

Ideally the above 4 maps would have their attributes spatially joined in R, but since I cannot figure out a nice way to do that, I did it quick and nicely in QGIS:

mea_map <- read_sf("model_data/input/measure/measure_comola_map.shp")
myfill <- hcl.colors(5, palette = "Fall")
NA_ahpkl <- mea_map$A_HPKL %>% is.na() %>% which()
mea_map$A_HPKL[NA_ahpkl] <- 0
mapview(mea_map, zcol = "A_HPKL", col.regions = myfill)

21.2 Determining Classification

We plan to do the following farm based classification:

70 farms in the catchment could implement gullies:

gully_farms <- mea_map %>% filter(gully > 0) %>% pull(GARDSNUMME) %>% unique() %>% length()
gully_farms

## [1] 70

68 farms can implement buffers:

buffer_farms <- mea_map %>% filter(buffer_6m_ > 0) %>% pull(GARDSNUMME) %>% unique() %>% length()
buffer_farms

## [1] 68

And we can implement 3 ponds.

dams <- mea_map %>% filter(dam > 0) %>% pull(dam) %>% unique() %>% length()
dams

## [1] 3

We have 501 agricultural fields to apply no-till to.

non_agri <- c("frst", "past", "rngb", "urml", "utrn", "wetf")
mea_map %>% filter((type %in% non_agri) == FALSE) %>% pull(type) %>% unique() %>% length()

## [1] 501

If we group these by farm level, we have 82

mea_map %>% filter((type %in% non_agri) == FALSE) %>% pull(GARDSNUMME) %>% unique() %>% length()

## [1] 82

If we expand these to erosion risk, we have:

class1 = mea_map %>% filter((type %in% non_agri) == FALSE) %>% filter(A_HPKL == 1) %>% pull(GARDSNUMME) %>% unique() %>% length()
class2 = mea_map %>% filter((type %in% non_agri) == FALSE) %>% filter(A_HPKL == 2) %>% pull(GARDSNUMME) %>% unique() %>% length()
class3 = mea_map %>% filter((type %in% non_agri) == FALSE) %>% filter(A_HPKL == 3) %>% pull(GARDSNUMME) %>% unique() %>% length()
class4 = mea_map %>% filter((type %in% non_agri) == FALSE) %>% filter(A_HPKL == 4) %>% pull(GARDSNUMME) %>% unique() %>% length()
class0 = mea_map %>% filter((type %in% non_agri) == FALSE) %>% filter(A_HPKL == 0) %>% pull(GARDSNUMME) %>% unique() %>% length()
hrus_df <- data.frame (
  c1 = class1,
  c2 = class2,
  c3 = class3,
  c4 = class4,
  c0 = class0)
DT::datatable(hrus_df)

In total, that would give us 363 measures to optimize: (well under budget)

  class1 + # low erosion risk lowtill farms
  class2 + # medium erosion risk lowtill farms
  class3 + # high erosion risk lowtill farms
  class4 + # very erosion risk lowtill farms
  class0 + # no erosion risk lowtill farms (not mapped)
  gully_farms + # gully farms
  buffer_farms + # buffer strip farms
  dams    # ponds

## [1] 363

21.3 Applying Classification

First we need to allocate incremental IDs to the measures:

class_1_alloc <- c(1:class1)
class_2_alloc <- c((last(class_1_alloc)+1):(last(class_1_alloc)+class2))
class_3_alloc <- c((last(class_2_alloc)+1):(last(class_2_alloc)+class3))
class_4_alloc <- c((last(class_3_alloc)+1):(last(class_3_alloc)+class4))
class_0_alloc <- c((last(class_4_alloc)+1):(last(class_4_alloc)+class0))
gully_alloc <- c((last(class_0_alloc)+1):(last(class_0_alloc)+gully_farms))
buffer_alloc <- c((last(gully_alloc)+1):(last(gully_alloc)+buffer_farms))
dam_alloc <- c((last(buffer_alloc)+1):(last(buffer_alloc)+dams))
full_Series <- c(
  class_1_alloc,
  class_2_alloc,
  class_3_alloc,
  class_4_alloc,
  class_0_alloc,
  gully_alloc,
  buffer_alloc,
  dam_alloc)
full_Series %>% length() # same length?

## [1] 363

(full_Series-c(1:363) == 0) %>% all() # all incremental?

## [1] TRUE

Now we need to translate the gardnummers to the allocated measure IDs..

21.3.1 Gullies (grassed water ways)

Measure range from 223 to 292

# Grabs the measure codes for each gully measure
gully_codes <- mea_map %>% filter(gully > 0) %>% pull(gully) 
# grabs the farm IDs for each gully measure
gully_farm_ids <- mea_map %>% filter(gully > 0) %>% pull(GARDSNUMME)
# Grabs the unique farms which have gully measures
gully_farm_uniq <- gully_farm_ids %>% unique() %>% sort()
# A look up table from gully farm IDs to old measure IDs
mea_to_farm <- data.frame(farm_id = gully_farm_ids, measure_id = gully_codes)
# A look up table from farm IDs which have gullies, to the gully-allocated IDs
farm_to_alloc <- data.frame(farm_id = gully_farm_uniq, gully_alloc = gully_alloc)
# A look up table combined from previous two
lookup <- left_join(farm_to_alloc,mea_to_farm, by = "farm_id", keep = FALSE, )
## Re-coding IDs
# https://stackoverflow.com/questions/66231321/recode-values-based-on-look-up-table-with-dplyr-r
gully_df = mea_map %>% st_drop_geometry %>% filter(gully > 0) %>% select(id, name, GARDSNUMME, gully)
recode_vec <- lookup$gully_alloc
names(recode_vec) <- lookup$measure_id
# reclassifying
gully_reclassed <- gully_df %>% mutate(gully_new = recode(gully, !!!recode_vec))
datatable(gully_reclassed)

21.3.2 Buffer Strips

Same procedure, measure IDs 293 to 360

# Grabs the measure codes for each buffer_6m_  measure
buffer_codes <- mea_map %>% filter(buffer_6m_  > 0) %>% pull(buffer_6m_ ) 
# grabs the farm IDs for each buffer measure
buffer_farm_ids <- mea_map %>% filter(buffer_6m_  > 0) %>% pull(GARDSNUMME)
# Grabs the unique farms which have buffer_6m_  measures
buffer_farm_uniq <- buffer_farm_ids %>% unique() %>% sort()
# A look up table from buffer farm IDs to old measure IDs
mea_to_farm <- data.frame(farm_id = buffer_farm_ids, measure_id = buffer_codes)
# A look up table from farm IDs which have buffers, to the buffer-allocated IDs
farm_to_alloc <- data.frame(farm_id = buffer_farm_uniq, buffer_alloc = buffer_alloc)
# A look up table combined from previous two
lookup <- left_join(farm_to_alloc,mea_to_farm, by = "farm_id", keep = FALSE, )
## Recoding IDs
# https://stackoverflow.com/questions/66231321/recode-values-based-on-look-up-table-with-dplyr-r
buffer_df = mea_map %>% st_drop_geometry %>% filter(buffer_6m_ > 0) %>% select(id, name, GARDSNUMME, buffer_6m_)
recode_vec <- lookup$buffer_alloc
names(recode_vec) <- lookup$measure_id
# reclassifiying
buffer_reclassed <- buffer_df %>% mutate(buffer_new = recode(buffer_6m_, !!!recode_vec))
datatable(buffer_reclassed)

21.3.3 Ponds

# Grabs the measure codes for each dam measure
pond_codes <- mea_map %>% filter(dam > 0) %>% pull(dam) 
# grabs the farm IDs for each pond measure
pond_farm_ids <- mea_map %>% filter(dam > 0) %>% pull(GARDSNUMME)
# Grabs the unique farms which have pond measures
pond_farm_uniq <- pond_farm_ids %>% unique() %>% sort()
# A look up table from pond farm IDs to old measure IDs
mea_to_farm <- data.frame(farm_id = pond_farm_ids, measure_id = pond_codes)
# A look up table from farm IDs which have gullies, to the pond-allocated IDs
farm_to_alloc <- data.frame(farm_id = pond_farm_uniq, pond_alloc = dam_alloc)
# A look up table combined from previous two
lookup <- left_join(farm_to_alloc,mea_to_farm, by = "farm_id", keep = FALSE, )
## Re-coding IDs
# https://stackoverflow.com/questions/66231321/recode-values-based-on-look-up-table-with-dplyr-r
pond_df = mea_map %>% st_drop_geometry() %>% filter(dam > 0) %>% select(id, name, GARDSNUMME, dam)
recode_vec <- lookup$pond_alloc
names(recode_vec) <- lookup$measure_id
# reclassifying
pond_reclassed <- pond_df %>% mutate(pond_new = recode(dam, !!!recode_vec))
datatable(pond_reclassed)

21.3.4 NoTill (management measure)

Re-classify for each erosion risk level

custom_reclass <- function(mea_map, erosion_class, mea_alloc) {
  non_agri <- c("frst", "past", "rngb", "urml", "utrn", "wetf")
  # farm codes of farms with class 1 erosion agricultural fields.
  relevant_data <- mea_map %>% st_drop_geometry %>% filter((type %in% non_agri) == FALSE) %>% filter(A_HPKL == erosion_class) %>% select(id, name, GARDSNUMME, type)
  class_X_farms <- relevant_data$GARDSNUMME %>% unique() %>% sort()
  class_X_hrus <- relevant_data$id %>% unique() %>% sort()
  recode_vec <- mea_alloc
  names(recode_vec) <- class_X_farms
  # reclassifying
  class_X_reclass <- relevant_data %>% mutate(notill = recode(GARDSNUMME, !!!recode_vec))
  return(class_X_reclass)
}
notill_1 <- custom_reclass(mea_map, erosion_class = 1, mea_alloc = class_1_alloc)
notill_2 <- custom_reclass(mea_map, erosion_class = 2, mea_alloc = class_2_alloc)
notill_3 <- custom_reclass(mea_map, erosion_class = 3, mea_alloc = class_3_alloc)
notill_4 <- custom_reclass(mea_map, erosion_class = 4, mea_alloc = class_4_alloc)
notill_0 <- custom_reclass(mea_map, erosion_class = 0, mea_alloc = class_0_alloc)
# notill_4$measure_id %>% unique() %>% length() == 10 # check
datatable(notill_3)

Merge into single dataframe

full_notill <- rbind((notill_1 %>% select(name, notill)),
                     (notill_2 %>% select(name, notill)),
                     (notill_3 %>% select(name, notill)),
                     (notill_4 %>% select(name, notill)),
                     (notill_0 %>% select(name, notill)))
full_notill <- full_notill[order(full_notill$name),]
datatable(full_notill)

21.4 HRU Scenario Table

We now need to generate the scenario HRU table for COMOLA.

Find source code here: UFZ-CLOUD

library(sf)
library(data.table)
library(stringr)
library(tidyverse)

## [1] "sf 1.0.18"

## [1] "data.table 1.16.2"

## [1] "stringr 1.5.1"

## [1] "tidyverse 2.0.0"

# all OPTAIN measures
measures_optain <- data.frame('nswrm_code' = c('buffer',
                                               'edgefilter',
                                               'hedge',
                                               'grassslope',
                                               'grassrchrg',
                                               'pond',
                                               'afforest',
                                               'floodres',
                                               'channres',
                                               'swale',
                                               'wetland',
                                               'cdrain',
                                               'terrace',
                                               'notill',
                                               'lowtill',
                                               'lowtillcc',
                                               'mulching',
                                               'subsoiling',
                                               'rotation',
                                               'intercrop',
                                               'covercrop',
                                               'earlysow',
                                               'droughtplt'),
                              'description' = c('Riparian buffers',
                                                'Edge-of-field filter strips',
                                                'Hedges/Field division',
                                                'Grasland cover on erosive slopes',
                                                'Grasland cover in recharge area',
                                                'Retention/detention ponds',
                                                'Afforestation',
                                                'Floodplain restoration',
                                                'Channel restoration',
                                                'Swales',
                                                'Constructed wetlands',
                                                'Controlled drainage',
                                                'Terracing',
                                                'No-till agriculture',
                                                'Low-till agriculture',
                                                'Low-till agriculture combined with cover crops',
                                                'Mulching',
                                                'Subsoiling',
                                                'Crop rotation',
                                                'Intercropping',
                                                'Green cover/ catch crops',
                                                'Early sowing',
                                                'Drought-resistant plants'))
datatable(measures_optain)

21.4.1 Merging into Landuse Map

“The land use shapefile must include columns for each measure (names according to the list above!). For each measure the id (as integer) of potential sites for implementation must be given. Multiple hrus can have the same measure id, e.g. if an individual grassed waterway consists of multiple hrus. If two individual measures of the same type (e.g. two crossing grassed waterways) occur on the same hru, provide both ids separated by ’_‘, e.g. ’5_6’ (see also example data). Define NA in case there is no potential site for a measure (please study the example lu shapefile)”

For our case, we do not need to consider the 5_6 thing as we have no doubled same-measure groupings.

Lessons learned: In the future if we would want to do scenarios on ie. the length of the grassed water ways, or the width of the buffer strips, we would have to incorporate this 5_6 thing.

lu_mea <- mea_map %>% select(name, type)
lu_mea <- left_join(lu_mea, buffer_reclassed %>% 
                      select(name, buffer = buffer_new), by= "name")
lu_mea <- left_join(lu_mea, gully_reclassed %>% 
                      select(name, grassslope = gully_new), by= "name")
lu_mea <- left_join(lu_mea, pond_reclassed %>% 
                      select(name, wetland = pond_new), by= "name")
lu_mea <- left_join(lu_mea, full_notill, by= "name")
lu_mea <- lu_mea %>%  select(name, type, grassslope, notill, buffer, wetland, geometry)
write_sf(lu_mea, "model_data/input/measure/lu_measures_COMOLA.shp")
lu_mea %>% st_drop_geometry() %>% datatable()

Now to use the provided script to generate the measure_location_cs10.csv

lu_shp <- "model_data/input/measure/lu_measures_COMOLA.shp"
hru_shp <- "model_data/cs10_setup/optain-cs10/data/vector/hru.shp"

# generate scenario hru table (just run the code)--------------------------

is.integer0 <- function(x)
{
  is.integer(x) && length(x) == 0L
}

# read in shapefiles
lu <- read_sf(lu_shp)
hru_swat <- read_sf(hru_shp)

# assign measures to SWAT+ hru names
lu_swat <- st_drop_geometry(st_join(hru_swat, lu %>% select(type, grassslope, notill, buffer, wetland), largest=T))

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

hru <- data.frame('hru'=lu_swat$name)

# generate hru-scenario table
meas.col <- which(colnames(lu_swat) %in% measures_optain$nswrm_code)
hru_scn <- data.frame('hru'=NA, 'nswrm'=NA, 'name'=NA)

for(i in 1:length(meas.col)){
  idx <- which(!is.na(lu_swat[,meas.col[i]]))
  meas_ <- colnames(lu_swat)[meas.col[i]]
  hru_ <- hru
  hru_$nswrm <- NA
  hru_$nswrm[idx] <- meas_
  hru_$name <- NA
  hru_$name[idx] <- paste0(meas_,'_', as.matrix(lu_swat[idx,meas.col[i]]))
  # different grassed waterways on the same hru?
  multiple <- which(str_count(hru_$name, "_") > 1)
  if(!is.integer0(multiple)){
    scn_multiple <- str_split_fixed(hru_$name[multiple],"_",3)
    hru_$name[multiple] <- paste0(meas_,'_',scn_multiple[,2])
    hru_add <- data.frame('hru' = hru_$hru[multiple],
                          'nswrm' = meas_,
                          'name' = paste0(meas_,'_',scn_multiple[,3]))
    hru_ <- rbind.data.frame(hru_, hru_add)
  }
  hru_ <- hru_[-which(is.na(hru_$name)),]
  hru_scn <- rbind.data.frame(hru_scn, hru_)
}

hru_scn <- hru_scn[-1,]
hru_scn$obj_id <- as.numeric(substr(hru_scn$hru,4,nchar(hru_scn$hru)))

hru_scn2 <- hru_scn %>% 
  group_by(name, nswrm) %>%
  summarize(obj_id = paste(sort(unique(obj_id)),collapse=", "))

## `summarise()` has grouped output by 'name'. You can override using the
## `.groups` argument.

hru_scn2$id <- c(1:length(hru_scn2$name))
hru_scn2 <- hru_scn2[,c(4,1,2,3)]

# write hru-scenario table
write_csv(hru_scn2, 'model_data/input/measure/measure_location_CS10.csv', quote = "needed", na = '')

# admire our result:
hru_scn2 %>% datatable()

# clean up for next chapter
rm(list=ls())