Section 21 Measure Classification
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.
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:
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_farms68 farms can implement buffers:
buffer_farms <- mea_map %>% filter(buffer_6m_ > 0) %>% pull(GARDSNUMME) %>% unique() %>% length()
buffer_farmsAnd we can implement 3 ponds.
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()If we group these by farm level, we have 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 # ponds21.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?
(full_Series-c(1:363) == 0) %>% all() # all incremental?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
21.4 HRU Scenario Table
We now need to generate the scenario HRU table for COMOLA.
Find source code here: UFZ-CLOUD
# 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))
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=", "))
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 = '')