Section 6 Landuse.lum update

This section covers the modifications made to the landuse file. For some reason it was written in a tutorial fashion, as if the reader were new to R.

require(tidyverse)
require(reshape2)
require(sf)
require(DT)
require(dplyr) # require dplyr last to overwrite plyr count()
require(ggplot2)

ggplot theme(s):

sf_theme <- theme(axis.text.x=element_blank(), #remove x axis labels
        axis.ticks.x=element_blank(), #remove x axis ticks
        axis.text.y=element_blank(),  #remove y axis labels
        axis.ticks.y=element_blank()  #remove y axis ticks
        )

6.1 Data preperation

We read in the land use file. We want to set skip = 1 to ignore the SWAT+editor text, and set header = T. We then convert it to a tibble format for better printing to console

landuse_lum <-
  read.table("model_data/cs10_setup/swat_input/landuse.lum", skip = 1, header = T) %>% tibble::as_tibble()

head(landuse_lum)

## # A tibble: 6 × 14
##   name       cal_group plnt_com mgt   cn2   cons_prac urban urb_ro ov_mann tile 
##   <chr>      <chr>     <chr>    <chr> <chr> <chr>     <chr> <chr>  <chr>   <chr>
## 1 a_001f_1_… null      nopl_co… null  null  null      null  null   null    mw24…
## 2 a_001f_2_… null      nopl_co… null  null  null      null  null   null    mw24…
## 3 a_001f_3_… null      nopl_co… null  null  null      null  null   null    mw24…
## 4 a_001f_4_… null      nopl_co… null  null  null      null  null   null    mw24…
## 5 a_001f_5_… null      nopl_co… null  null  null      null  null   null    mw24…
## 6 a_001f_6_… null      nopl_co… null  null  null      null  null   null    mw24…
## # ℹ 4 more variables: sep <chr>, vfs <chr>, grww <chr>, bmp <chr>

Our field_id for our cropland does not match name of landuse_lum so we need to parse it out, we can do this many ways, but a safe way is to split it by “_” and combine the first 3 splits with that same underscore:

splitted <- landuse_lum$name %>% str_split("_")

landuse_lum$field_id <-
  paste(splitted %>% map(1), splitted %>% map(2), splitted %>% map(3), sep = "_")

head(landuse_lum$field_id)

## [1] "a_001f_1" "a_001f_2" "a_001f_3" "a_001f_4" "a_001f_5" "a_001f_6"

But wait, this does not work for our “non-fields”. So lets find out which ones they are, and set them to NA

not_fields <- which(!grepl(x=landuse_lum$field_id, "a_"))
landuse_lum$field_id[not_fields] <- NA

So, what non-fields do we have?

landuse_lum$name[not_fields]

## [1] "frst_lum" "past_lum" "rngb_lum" "urml_lum" "utrn_lum" "wetf_lum"

Good to know. We’ll keep that in mind.

6.2 Assigning landuse values

6.2.1 Setting `cal_group`

There is no info on this column, so we are going to leave it as null. Dicussion in Issue #18

6.2.2 Setting `plnt_com`

This step will be done by SWATFarmR in Section 11. Dicussed in Issue #15

6.2.3 Setting `mgt`

This step will be done by SWATFarmR in Section 11.

Discussed in Issue #14

6.2.4 Setting `urb_ro`

We want to set the urb_ro column for all urban land uses (in our case this would be urml and utrn) to usgs_reg.

This has been discussed in Issue #6

We can do this like so:

landuse_lum$urb_ro[which(landuse_lum$name %in% c("urml_lum", "utrn_lum"))] <-
  "usgs_reg"

## # A tibble: 2 × 15
##   name     cal_group plnt_com mgt   cn2   cons_prac urban urb_ro   ov_mann tile 
##   <chr>    <chr>     <chr>    <chr> <chr> <chr>     <chr> <chr>    <chr>   <chr>
## 1 urml_lum null      null     null  null  null      null  usgs_reg null    null 
## 2 utrn_lum null      null     null  null  null      null  usgs_reg null    null 
## # ℹ 5 more variables: sep <chr>, vfs <chr>, grww <chr>, bmp <chr>,
## #   field_id <chr>

Very good. In our case this was only two – could be done by hand, but that will not be the case for all of our land uses.

6.2.5 Setting `urban`

This one is easy, we set our urban column to an urban parameter set of the same name. The rest we leave as null

This has been discussed in Issue #5

landuse_lum$urban[which(landuse_lum$name =="utrn_lum")] <- "utrn"
landuse_lum$urban[which(landuse_lum$name =="urml_lum")] <- "urml"

## # A tibble: 1 × 15
##   name     cal_group plnt_com mgt   cn2   cons_prac urban urb_ro   ov_mann tile 
##   <chr>    <chr>     <chr>    <chr> <chr> <chr>     <chr> <chr>    <chr>   <chr>
## 1 urml_lum null      null     null  null  null      urml  usgs_reg null    null 
## # ℹ 5 more variables: sep <chr>, vfs <chr>, grww <chr>, bmp <chr>,
## #   field_id <chr>

Lets make sure other land uses still have null:

## # A tibble: 1 × 15
##   name     cal_group plnt_com  mgt   cn2   cons_prac urban urb_ro ov_mann tile 
##   <chr>    <chr>     <chr>     <chr> <chr> <chr>     <chr> <chr>  <chr>   <chr>
## 1 frst_lum null      nopl_comm null  null  null      null  null   null    null 
## # ℹ 5 more variables: sep <chr>, vfs <chr>, grww <chr>, bmp <chr>,
## #   field_id <chr>

Looks good.

6.2.6 Setting Manning’s n (`ovn`)

This has been discussed in Issue #13

Lets get the easy ones out of the way

landuse_lum$ov_mann[which(landuse_lum$name =="past_lum")] <- "densegrass"
landuse_lum$ov_mann[which(landuse_lum$name =="rngb_lum")] <- "rangeland_20cover"
landuse_lum$ov_mann[which(landuse_lum$name =="urml_lum")] <- "urban_rubble"
landuse_lum$ov_mann[which(landuse_lum$name =="urtn_lum")] <- "urban_asphalt"

Don’t fall asleep yet! For wetf we want forest_heavy but with a higher value. this means we need to add a new entry. And since we are doing this fancy Rmarkdown stuff, we’re going to do it with code! Lets jump in and get at this file. (REMEBER TO REMOVE THE TEMP)

ovn_table_path <- "model_data/cs10_setup/swat_input/ovn_table.lum"
ovn_table_path_new <- "model_data/cs10_setup/swat_input_mod/ovn_table.lum"

ovn_table <- readLines(ovn_table_path)

Good, now where is this forest heavy entry? and what does the format look like?

index <- grepl(x=ovn_table, "forest_heavy") %>% which %>% min()

ovn_table[c(2, index)] %>% print()

## [1] "name                  ovn_mean       ovn_min       ovn_max  description" 
## [2] "forest_heavy           0.80000       0.70000       0.90000  Forest_heavy"

Now, lets make our own and add it in. But only if it doesn’t exist it (Like if the script has been run before…)

if(grepl(x = ovn_table, "forest_heavy_cs10") %>% which %>% length() == 0) {
  entry <-
    "forest_heavy_cs10      0.90000       0.80000       0.95000  Forest_heavy_mod"
  ovn_table <- c(ovn_table, entry)
}

Did it work? Yes:

ovn_table %>% tail()

## [1] "forest_med             0.60000       0.50000       0.70000  Forest_medimum_good"
## [2] "forest_heavy           0.80000       0.70000       0.90000  Forest_heavy"       
## [3] "urban_asphalt          0.01100       0.01100       0.01100  Urban_asphalt"      
## [4] "urban_concrete         0.01200       0.01200       0.01200  Urban_concrete"     
## [5] "urban_rubble           0.02400       0.02400       0.02400  Urban_rubble"       
## [6] "forest_heavy_cs10      0.90000       0.80000       0.95000  Forest_heavy_mod"

Now lets write this new table

writeLines(ovn_table, con = ovn_table_path_new)

And now we can enter our wetland class:

landuse_lum$ov_mann[which(landuse_lum$name =="wetf_lum")] <- "forest_heavy_cs10"

For the fields we are going to need the crop rotation info which we created in section 5.

crop_rotation <- read_csv("model_data/input/crop/cs10_crop_rotation.csv",
    show_col_types = F)

head(crop_rotation)

## # A tibble: 6 × 5
##   field    y_2016 y_2017 y_2018 y_2019
##   <chr>    <chr>  <chr>  <chr>  <chr> 
## 1 a_025f_1 wwht   wwht   wwht   wwht  
## 2 a_025f_3 wwht   wwht   wwht   wwht  
## 3 a_025f_4 wwht   wwht   wwht   wwht  
## 4 a_025f_2 wwht   wwht   wwht   wwht  
## 5 a_105f_1 wwht   wwht   wwht   wwht  
## 6 a_105f_4 wwht   wwht   wwht   wwht

We have decided to classify ovn based on the degree to which a crop rotation contains meadow, conventional crops (wwht, pota), and conservation crops (all others). To do this we need to analyze the crop rotation.

Now it gets a little complicated, but trust me its not actually as bad as it looks. We need to count how many times certain crops show up in the crop rotation of a certain field. Lets do it for conventional crops first.

Now, lets get into it

conv_crop_count <- crop_rotation %>% melt(. ,"field") %>% group_by(field) %>%
  filter(value %in% c("wwht", "pota")) %>% dplyr::count() %>% dplyr::rename(conv = n)

head(conv_crop_count)

## # A tibble: 6 × 2
## # Groups:   field [6]
##   field     conv
##   <chr>    <int>
## 1 a_001f_1     4
## 2 a_001f_2     4
## 3 a_001f_3     4
## 4 a_001f_4     4
## 5 a_001f_5     4
## 6 a_001f_6     4

What happened here? well we took our crop_rotation, and melted it by field using melt. This function converts the data into tidy format (Wickham 2014). This format makes it easier to apply the following operations on a field basis. What does this look like?

crop_rotation %>% melt(. ,"field") %>% arrange(field) %>% head()

##      field variable value
## 1 a_001f_1   y_2016  wwht
## 2 a_001f_1   y_2017  wwht
## 3 a_001f_1   y_2018  wwht
## 4 a_001f_1   y_2019  wwht
## 5 a_001f_2   y_2016  wwht
## 6 a_001f_2   y_2017  wwht

melt is a function from reshape2 which is why we required it. And what does that “.” mean in there, to the left of "field"?? That is a code word for the object to the left of the “pipe” (in this case crop_rotation). the pipe is “%>%” and passes the object to the left of it, into the function to the right of it. look it up!

When arranged by field, we can see that every field gets one entry per crop per year. This is a good format to count how many times we have a certain type of crop on a field.

Next we group_by the field – this means all the following operations will be done on a field basis. Then we filter our value (which is the crop name). For conventional we needed to filter in any fields with the crop "wwht" or "pota".

What does this look like?

crop_rotation %>% melt(. ,"field") %>% group_by(field) %>%
  filter(value %in% c("wwht", "pota"))

## # A tibble: 1,519 × 3
## # Groups:   field [415]
##    field    variable value
##    <chr>    <fct>    <chr>
##  1 a_025f_1 y_2016   wwht 
##  2 a_025f_3 y_2016   wwht 
##  3 a_025f_4 y_2016   wwht 
##  4 a_025f_2 y_2016   wwht 
##  5 a_105f_1 y_2016   wwht 
##  6 a_105f_4 y_2016   wwht 
##  7 a_105f_2 y_2016   wwht 
##  8 a_105f_3 y_2016   wwht 
##  9 a_187f_1 y_2016   wwht 
## 10 a_206f_1 y_2016   wwht 
## # ℹ 1,509 more rows

Looks good. We only have crops with winter wheat and potatoes, for every year. Exactly what we need, now we just need to count() them.

conv_crop_count <- crop_rotation %>% melt(. ,"field") %>% group_by(field) %>%
  filter(value %in% c("wwht", "pota")) %>% dplyr::count()

head(conv_crop_count)

## # A tibble: 6 × 2
## # Groups:   field [6]
##   field        n
##   <chr>    <int>
## 1 a_001f_1     4
## 2 a_001f_2     4
## 3 a_001f_3     4
## 4 a_001f_4     4
## 5 a_001f_5     4
## 6 a_001f_6     4

And we are back where we started! See, not that complicated. One last thing we need to do is rename n to conv, we do that like so:

conv_crop_count <- crop_rotation %>% melt(. ,"field") %>% group_by(field) %>% 
  filter(value %in% c("wwht", "pota")) %>% dplyr::count() %>% dplyr::rename(conv = n)

head(conv_crop_count)

## # A tibble: 6 × 2
## # Groups:   field [6]
##   field     conv
##   <chr>    <int>
## 1 a_001f_1     4
## 2 a_001f_2     4
## 3 a_001f_3     4
## 4 a_001f_4     4
## 5 a_001f_5     4
## 6 a_001f_6     4

Now lets go ahead and do the same thing for the two other categories: cons andmeadow.

meadow_crop_count <- crop_rotation %>% melt(. ,"field") %>% group_by(field) %>%
  filter(value == "meadow" ) %>% dplyr::count() %>%
  dplyr::rename(meadow = n)

cons_crop_count <- crop_rotation %>% melt(. ,"field") %>% group_by(field) %>%
  filter(!(value %in% c("wwht", "pota", "meadow"))) %>% dplyr::count() %>%
  dplyr::rename(cons = n)

cons_crop_count %>% head()

## # A tibble: 6 × 2
## # Groups:   field [6]
##   field     cons
##   <chr>    <int>
## 1 a_002f_1     2
## 2 a_002f_2     2
## 3 a_002f_3     4
## 4 a_002f_4     4
## 5 a_002f_5     2
## 6 a_002f_6     4

meadow_crop_count %>% head()

## # A tibble: 6 × 2
## # Groups:   field [6]
##   field    meadow
##   <chr>     <int>
## 1 a_012f_3      4
## 2 a_012f_6      4
## 3 a_016f_3      4
## 4 a_017f_2      4
## 5 a_017f_3      4
## 6 a_046f_1      4

meadow was a simple filter, all we have to do was grab crops with the meadow name. For cons crops, we just grabbed the remaining crops that were not conv or meadow.

Great. We have these 3 separate, we need to combine them. we can do that with left_join and join by the field column which contains our IDs

# create a base dataframe to join to
crop_fractions <- crop_rotation %>% dplyr::select(field) %>% distinct()

# join our 3 data frames
crop_fractions <- left_join(crop_fractions, conv_crop_count, by = "field")
crop_fractions <- left_join(crop_fractions, cons_crop_count, by = "field")
crop_fractions <- left_join(crop_fractions, meadow_crop_count, by = "field")

# lets have a look
crop_fractions %>% head()

## # A tibble: 6 × 4
##   field     conv  cons meadow
##   <chr>    <int> <int>  <int>
## 1 a_025f_1     4    NA     NA
## 2 a_025f_3     4    NA     NA
## 3 a_025f_4     4    NA     NA
## 4 a_025f_2     4    NA     NA
## 5 a_105f_1     4    NA     NA
## 6 a_105f_4     4    NA     NA

That does not look good! when it seems like when the count is 0, it is returned as NA. Lets fix that…

crop_fractions <- crop_fractions %>%
  mutate(cons = ifelse(is.na(cons), 0, cons))
crop_fractions <- crop_fractions %>%
  mutate(conv = ifelse(is.na(conv), 0, conv))
crop_fractions <- crop_fractions %>%
  mutate(meadow = ifelse(is.na(meadow), 0, meadow))

What are we doing here? we are mutating the the 3 columns using an ifelse statement. The statement is simple. If the value is.na then we set it to 0. Otherwise, we set it to the same value it had before.

did it work?

crop_fractions %>% head()

## # A tibble: 6 × 4
##   field     conv  cons meadow
##   <chr>    <dbl> <dbl>  <dbl>
## 1 a_025f_1     4     0      0
## 2 a_025f_3     4     0      0
## 3 a_025f_4     4     0      0
## 4 a_025f_2     4     0      0
## 5 a_105f_1     4     0      0
## 6 a_105f_4     4     0      0

Very good. Now we need to decide what to classify our fields at. lets Define some functions to do that.

is_meadow <- function(conv, cons, meadow) {
  ((meadow > cons) & (meadow > conv)) %>% return()
}

is_cons <- function(conv, cons, meadow) {
  ((cons >= meadow) & (cons > conv)) %>% return()
}

is_conv <- function(conv, cons, meadow) {
  ((conv >= meadow) & (conv >= cons)) %>% return()
}

The & sign means that both conditions need to be met. and >= you should know already. Lets use those functions in action. lets do the meadow first. We will create a new dataframe from crop fractions, named field_class.

field_class <- crop_fractions %>% mutate(meadow = is_meadow(conv,cons,meadow))

head(field_class)

## # A tibble: 6 × 4
##   field     conv  cons meadow
##   <chr>    <dbl> <dbl> <lgl> 
## 1 a_025f_1     4     0 FALSE 
## 2 a_025f_3     4     0 FALSE 
## 3 a_025f_4     4     0 FALSE 
## 4 a_025f_2     4     0 FALSE 
## 5 a_105f_1     4     0 FALSE 
## 6 a_105f_4     4     0 FALSE

Ok, so none of those first 6 fields are meadow dominated. What about cons? (lets stick with our field_class dataframe and just keep adding on)

field_class <- field_class %>% mutate(cons = is_cons(conv,cons,meadow))
head(field_class)

## # A tibble: 6 × 4
##   field     conv cons  meadow
##   <chr>    <dbl> <lgl> <lgl> 
## 1 a_025f_1     4 FALSE FALSE 
## 2 a_025f_3     4 FALSE FALSE 
## 3 a_025f_4     4 FALSE FALSE 
## 4 a_025f_2     4 FALSE FALSE 
## 5 a_105f_1     4 FALSE FALSE 
## 6 a_105f_4     4 FALSE FALSE

Nope, not cons either. Then it must be conv dominant.

field_class <- field_class %>% mutate(conv = is_conv(conv,cons,meadow))

field_class %>% head()

## # A tibble: 6 × 4
##   field    conv  cons  meadow
##   <chr>    <lgl> <lgl> <lgl> 
## 1 a_025f_1 TRUE  FALSE FALSE 
## 2 a_025f_3 TRUE  FALSE FALSE 
## 3 a_025f_4 TRUE  FALSE FALSE 
## 4 a_025f_2 TRUE  FALSE FALSE 
## 5 a_105f_1 TRUE  FALSE FALSE 
## 6 a_105f_4 TRUE  FALSE FALSE

Correct! now lets just double check that we didnt have any cases where all are TRUE or where all are FALSE

field_class %>% dplyr::select(conv, cons, meadow) %>% isTRUE %>% all()

## [1] FALSE

field_class %>% dplyr::select(conv, cons, meadow) %>% isFALSE() %>% all()

## [1] FALSE

Looks good to me! carry on:

We have our field classifications now, they will be very useful to us later. Lets pull them out of our dataframe to store them as a nice list.

cons_fields <- field_class %>% filter(cons) %>% dplyr::select(field) %>% pull()
conv_fields <- field_class %>% filter(conv) %>% dplyr::select(field) %>% pull()
meadow_fields <- field_class %>% filter(meadow) %>% dplyr::select(field) %>% pull()

cons_fields %>% head()

## [1] "a_182f_1" "a_182f_5" "a_182f_3" "a_182f_4" "a_182f_2" "a_182f_7"

conv_fields %>% head()

## [1] "a_025f_1" "a_025f_3" "a_025f_4" "a_025f_2" "a_105f_1" "a_105f_4"

meadow_fields %>% head()

## [1] "a_182f_6" "a_182f_8" "a_017f_3" "a_017f_2" "a_012f_3" "a_012f_6"

Fantastic! Now that was a big task, but it makes the next bit very easy. Lets assign the correct ovn to the types of fields we have classified:

6.2.6.1 Agricultural Fields: Meadow

landuse_lum$ov_mann[which(landuse_lum$field_id %in% meadow_fields)] <- "shortgrass"

6.2.6.2 Agricultural Fields: Conventional

landuse_lum$ov_mann[which(landuse_lum$field_id %in% conv_fields)] <- "convtill_nores"

6.2.6.3 Agricultural Fields: Conservational

landuse_lum$ov_mann[which(landuse_lum$field_id %in% cons_fields)] <- "falldisk_res"

6.2.6.4 Forest

Initially we wanted certain classes for the forest land use based on how good the soil quality was. But to do this, we would need multiple land uses, which we do not have. We will just use forest_medium for all frst. We can only change this if we go back to buildR and define more generic forest classes.

landuse_lum$ov_mann[which(landuse_lum$name =="frst_lum")] <- "forest_med"

6.2.6.5 Summary

With all that data processing out of the way, lets have a quick look at what we’ve done:

landuse_lum$ov_mann[which(landuse_lum$name == "past_lum")] <-
  "densegrass"
landuse_lum$ov_mann[which(landuse_lum$name == "rngb_lum")] <-
  "rangeland_20cover"
landuse_lum$ov_mann[which(landuse_lum$name == "urml_lum")] <-
  "urban_rubble"
landuse_lum$ov_mann[which(landuse_lum$name == "urtn_lum")] <-
  "urban_asphalt"
landuse_lum$ov_mann[which(landuse_lum$name == "frst_lum")] <-
  "forest_med"
landuse_lum$ov_mann[which(landuse_lum$field_id %in% meadow_fields)] <-
  "shortgrass"
landuse_lum$ov_mann[which(landuse_lum$field_id %in% conv_fields)] <-
  "convtill_nores"
landuse_lum$ov_mann[which(landuse_lum$field_id %in% cons_fields)] <-
  "falldisk_res"
landuse_lum$ov_mann[which(landuse_lum$name == "wetf_lum")] <-
  "forest_heavy_cs10"

Sweet. Next column..

6.2.7 Setting `cn2`:

This has been discussed in Issue Issue #1

You know the drill – I am sure you know how to read this code by now.

# Brush-brush-weed-grass_mixture_with_brush_the_major_element (poor)
landuse_lum$cn2[which(landuse_lum$name == "rngb_lum")] <- "brush_p"

# Woods (poor)
landuse_lum$cn2[which(landuse_lum$name == "wetf_lum")] <- "wood_p"

# Paved_streets_and_roads;_open_ditches_(incl._right-of-way)
landuse_lum$cn2[which(landuse_lum$name == "utrn_lum")] <- "paveroad"

# Paved_parking_lots_roofs_driveways_etc_(excl_right-of-way)
landuse_lum$cn2[which(landuse_lum$name == "urml_lum")] <- "urban"

# Woods (fair)
landuse_lum$cn2[which(landuse_lum$name == "frst_lum")] <- "wood_f"

# Pasture_grassland_or_range-continuous_forage_for_grazing
landuse_lum$cn2[which(landuse_lum$name == "past_lum")] <- "pastg_g"

# Meadow-continuous_grass_protected_from_grazing_mowed_for_hay
landuse_lum$cn2[which(landuse_lum$field_id %in% meadow_fields)] <-
  "pasth"

# Row_crops
landuse_lum$cn2[which(landuse_lum$field_id %in% cons_fields)] <-
  "rc_conterres_g"

# Row_crops
landuse_lum$cn2[which(landuse_lum$field_id %in% conv_fields)] <-
  "rc_strow_p"

6.2.8 Setting `cons_prac`

This has been discussed in Issue #4

For this, we need to add some custom entries to the database:

cons_prac_path <- "model_data/cs10_setup/swat_input/cons_practice.lum"
cons_prac_path_new <- "model_data/cs10_setup/swat_input_mod/cons_practice.lum"

cons_prac <- readLines(cons_prac_path)
cons_prac[1:3] %>% print()

## [1] "cons_practice.lum: written by SWAT+ editor v2.1.0 on 2023-05-30 10:39 for SWAT+ rev.60.5.4"
## [2] "name                    usle_p   slp_len_max  description"                                 
## [3] "up_down_slope          1.00000     121.00000  Up_and_down_slope"

Let us do it in a way so that it is only added if it doesn’t exist yet:

if(grepl(x = cons_prac, "agri_conv") %>% which %>% length() == 0) {
  entry <-
    "agri_conv          1.00000      60.00000  no_convervation"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "agri_part_conv") %>% which %>% length() == 0) {
  entry <-
    "agri_part_conv     0.85000      60.00000  75_percent_convential"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "agri_half") %>% which %>% length() == 0) {
  entry <-
    "agri_half          0.70000      50.00000  50_percent_convential"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "agri_part_cons") %>% which %>% length() == 0) {
  entry <-
    "agri_part_cons     0.50000      30.00000  75_percent_consveration"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "agri_cons") %>% which %>% length() == 0) {
  entry <-
    "agri_cons     0.30000      30.00000  100_percent_consveration"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "past_cons") %>% which %>% length() == 0) {
  entry <-
    "past_cons   0.1   60   pasture"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "rngb_cons") %>% which %>% length() == 0) {
  entry <-
    "rngb_cons   0.2   60   rangeland"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "frst_cons") %>% which %>% length() == 0) {
  entry <-
    "frst_cons   0.1   60   forest"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "urml_cons") %>% which %>% length() == 0) {
  entry <-
    "urml_cons   1   60   cs10urban"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "utrn_cons") %>% which %>% length() == 0) {
  entry <-
    "utrn_cons   1   60   road"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "cs10_meadow") %>% which %>% length() == 0) {
  entry <-
    "cs10_meadow   0.2   60   cs10_meadow"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "wetf_cons") %>% which %>% length() == 0) {
  entry <-
    "wetf_cons   0.05   30   wetlands"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "cs10_sed_pond") %>% which %>% length() == 0) {
  entry <-
    "cs10_sed_pond   0.1   30   cs10_sed_pond"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "cs10_cons_wetl") %>% which %>% length() == 0) {
  entry <-
    "cs10_cons_wetl   0.1   30   cs10_cons_wetl"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "cs10_buff_grass") %>% which %>% length() == 0) {
  entry <-
    "cs10_buff_grass   0.25   10   cs10_buff_grass"
  cons_prac <- c(cons_prac, entry)
}

if (grepl(x = cons_prac, "cs_10_buff_wood") %>% which %>% length() == 0) {
  entry <-
    "cs_10_buff_wood   0.2   10   cs_10_buff_wood"
  cons_prac <- c(cons_prac, entry)
}

And write the updated file:

writeLines(cons_prac, con = cons_prac_path_new)

Now we need to define the 0, 25, 50, 75, 100 percent conventional crops. We can use our old “crop_fractions” dataframe:

head(crop_fractions)

## # A tibble: 6 × 4
##   field     conv  cons meadow
##   <chr>    <dbl> <dbl>  <dbl>
## 1 a_025f_1     4     0      0
## 2 a_025f_3     4     0      0
## 3 a_025f_4     4     0      0
## 4 a_025f_2     4     0      0
## 5 a_105f_1     4     0      0
## 6 a_105f_4     4     0      0

Lets derive which fields get classified as what

p100_conv <- crop_fractions %>% filter(conv==4) %>% dplyr::select(field) %>% pull()
p75_conv <- crop_fractions %>% filter(conv==3) %>% dplyr::select(field) %>% pull()
p50_conv <- crop_fractions %>% filter(conv==2) %>% dplyr::select(field) %>% pull()
p25_conv <- crop_fractions %>% filter(conv==1) %>% dplyr::select(field) %>% pull()
p0_conv <- crop_fractions %>% filter(conv==0) %>% dplyr::select(field) %>% pull()

Now we can assign the all the land uses

# Agricultural
landuse_lum$cons_prac[which(landuse_lum$field_id %in% p100_conv)] <-
  "agri_conv"
landuse_lum$cons_prac[which(landuse_lum$field_id %in% p75_conv)] <-
  "agri_part_conv"
landuse_lum$cons_prac[which(landuse_lum$field_id %in% p50_conv)] <-
  "agri_half"
landuse_lum$cons_prac[which(landuse_lum$field_id %in% p25_conv)] <-
  "agri_part_cons"
landuse_lum$cons_prac[which(landuse_lum$field_id %in% p0_conv)] <-
  "agri_cons"

# Meadow
landuse_lum$cons_prac[which(field_class$meadow)] <-
  "cs10_meadow"

# Generic
landuse_lum$cons_prac[which(landuse_lum$name == "past_lum")] <-
  "past_cons"
landuse_lum$cons_prac[which(landuse_lum$name == "rngb_lum")] <-
  "rngb_cons"
landuse_lum$cons_prac[which(landuse_lum$name == "frst_lum")] <-
  "frst_cons"
landuse_lum$cons_prac[which(landuse_lum$name == "urml_lum")] <-
  "urml_cons"
landuse_lum$cons_prac[which(landuse_lum$name == "utrn_lum")] <-
  "utrn_cons"
landuse_lum$cons_prac[which(landuse_lum$name == "wetf_lum")] <-
  "wetf_cons"

# Measures
landuse_lum$cons_prac[which(landuse_lum$name == "cs10_sed_pond")] <-
  "cs10_sed_pond"
landuse_lum$cons_prac[which(landuse_lum$name == "cs10_cons_wetl")] <-
  "cs10_cons_wetl"
landuse_lum$cons_prac[which(landuse_lum$name == "cs10_buff_grass")] <-
  "cs10_buff_grass"
landuse_lum$cons_prac[which(landuse_lum$name == "cs_10_buff_wood")] <-
  "cs_10_buff_wood"

6.2.9 Setting `tile`

This column has already been completed by SWATbuildR in 2.4

6.2.10 Setting `sep`

This column has been left as null as discussed in Issue #9

6.2.11 Setting `vfs`

This column has been left as null as discussed in Issue #10

6.2.12 Setting `grww`

This column has been left as null as discussed in Issue #11

6.2.13 Setting `bmp`

This column has been left as null as discussed in Issue #12

6.3 Writing Changes

We are done with our modifications however, we need to remove the field_id column

landuse_lum <- landuse_lum %>% dplyr::select(-field_id)

lets take a look at the final product:

datatable(landuse_lum)

Lets write out changes in a new directory where we will store any SWAT+ input files which we have modified from their source BuildR form.

dir.create("model_data/cs10_setup/swat_input_mod", showWarnings = F)
new_lum_path <- "model_data/cs10_setup/swat_input_mod/landuse.lum"

write.table(landuse_lum, file = new_lum_path, sep = "\t", quote = F, row.names = F)

But wait – we need to keep that pesky header, otherwise the FarmR will be very angry with us.

lum_lines <- readLines(new_lum_path)

header <- "header header HEADER, delete me and you will regret it FOREVER!"

lum_lines2<- c(header,lum_lines)

writeLines(text = lum_lines2, con = new_lum_path)

Done!

Now that We’ve changed another input file, we should test if our model still runs. We will do that with our custom test_swat_mod() function. This function is basically the same as the previous test_swat() function, only that it incorporates the modified SWAT input files

source("model_data/code/test_swat_mod.R")

simout <- test_swat_mod()

cat(tail(simout), sep = "\n")

Related Milestone

landuse.lum update

Related minor issues

Don’t write rownames! - #50

Update the land use file - #40

References

Wickham, Hadley. 2014. “Tidy Data.” Journal of Statistical Software 59 (10). https://doi.org/10.18637/jss.v059.i10.

Section 6 Landuse.lum update

6.1 Data preperation

6.2 Assigning landuse values

6.2.1 Setting cal_group

6.2.2 Setting plnt_com

6.2.3 Setting mgt

6.2.4 Setting urb_ro

6.2.5 Setting urban

6.2.6 Setting Manning’s n (ovn)