%load_ext autoreload
%autoreload 2

pywatemsedem API

Introduction

This tutorial describes the use of the classes and functions of the pywatemsedem Python package. This python package functions as a model API to prepare and process data to create all input-files for the WaTEM/SEDEM. The Python package also contains functions for post-processing the modeloutput, yet this is not handled in this notebook.

Note:

The example data available in the subfolder data/langegracht is preclipped to the catchment shape. Do note that out-of-bound rasters and shapes are automatically clipped by the inputted catchment vector or raster data layer, as long as the clips fully overlay with the catchment outline (i.e. not missing values).

Imports and example data

import os
from pathlib import Path

The example data are located in the test data folder

inputdata_folder = Path(r"..") / ".." / "tests" /  "data"

inputdata_folder.exists()

Define WaTEM/SEDEM exe

Pick it up the latest version from https://github.com/watem-sedem/watem-sedem/releases and put the exe in the folder of this notebook (or in a .env-file)

from dotenv import load_dotenv, find_dotenv
load_dotenv(find_dotenv())

watem_sedem_binary = Path(os.environ.get("WATEMSEDEM"))

watem_sedem_binary

watem_sedem_binary.exists()

watem_sedem_binary

Generate catchment

from pywatemsedem.catchment import Catchment

Give your catchment a name

name_catchment = 'langegracht'

Initialize the catchment with a catchment vector definition and a dtm raster.

rst_dtm = inputdata_folder  / "dtm.tif"

vct_catchment = inputdata_folder / "catchm_langegracht.shp"

vct_catchment.exists()

str(vct_catchment.resolve())

Feed the name of the catchment, the outline vector, the dtm raster, the desired resolution, the desired coordinate definition, the nodata-value and the year.

catchment = Catchment(name_catchment, vct_catchment, rst_dtm, 20, 31370, -9999)

Plot the mask

catchment.mask.plot()

Plot DTM

catchment.dtm.plot()

Input rasters

Clip the land-use, K-factor and CN soil type maps and convert to the right data format by the functions below

rst_kfactor = inputdata_folder / "K3.tif"
rst_landuse = inputdata_folder / "basemap_landuse.tif"

catchment.kfactor = rst_kfactor

catchment.kfactor.plot(nodata=-9999)

catchment.landuse = rst_landuse
catchment.landuse.plot(nodata=-9999)

catchment.landuse.write("landuse.rst")

Input vectors

Rivers and infrastructure

vct_river = inputdata_folder / "river.shp"
vct_infra_poly = inputdata_folder / "infrastructure.shp"
vct_infra_line = inputdata_folder / "roads.shp"
vct_water = inputdata_folder / "water.shp"

Assigning the river vector will generate:

• A river raster (with values -1 for river, else nodata/0)

• A river segments (values for every river segment)

catchment.vct_river = vct_river

catchment.vct_river.plot()

By assigning the river vector, the river raster and segments are automatically created. You can access these rasters by:

catchment.river.plot(nodata=-9999)

You can also have a look at the extracted segment

catchment.segments.plot(nodata=0)

Assign and plot infrastructure polygons and lines

catchment.vct_infrastructure_buildings = vct_infra_poly
catchment.infrastructure_buildings.plot(nodata=-9999)

catchment.vct_infrastructure_roads = vct_infra_line
catchment.infrastructure_roads.plot(nodata=-9999)

catchment.infrastructure.plot(nodata=-9999)

catchment.vct_water = vct_water
catchment.water.plot(nodata=-9999)

The Choices-object

The second element we have to defines is the user choices, variables and other parameters. We have chosen to define these choices in a different class/object. In this way the defined choices can be re-used in calculations for e.g. a different catchment.

First, we need to initiate the user choices by importing the necessary classes

from pywatemsedem.choices import Choices, Options, Parameters, Extensions, ExtensionsParameters, Output

Make use of a preset default values for choices

default_choices = inputdata_folder / "userchoices.ini"

default_choices.exists()

options = Options()
options.read_values_from_ini(default_choices)

parameters = Parameters()
parameters.read_values_from_ini(default_choices)

extensions = Extensions()
extensions.read_values_from_ini(default_choices)

output = Output()
output.read_values_from_ini(default_choices)

extensionparameters = ExtensionsParameters(extensions)
extensionparameters.apply_defaults()
extensionparameters.ktc_low.read_value_from_ini(default_choices)
extensionparameters.ktc_high.read_value_from_ini(default_choices)
extensionparameters.ktc_limit.read_value_from_ini(default_choices)

choices = Choices(options, parameters, extensions, extensionparameters, output)

The Scenario-object

With the defined catchment and all the choices made by the user we can start to prepare all the necessary inputdata for a scenario.

A scenario-instance needs a valid catchment-instance and a valid choices-instance to be initiated. Here, we will use the catchment created above. Also as scenario-number is needed.

from pywatemsedem.scenario import Scenario

scenario_nr = 1
scenario = Scenario(catchment, 2019, scenario_nr, choices)

Add parcels information

scenario.vct_parcels = inputdata_folder / r"parcels.shp"

scenario.vct_parcels.geodata

scenario.parcels.plot(nodata=-9999)

Run model without any measures

scenario.composite_landuse = scenario.create_composite_landuse()

scenario.composite_landuse.plot()

scenario.cfactor = scenario.create_cfactor()

scenario.cfactor.plot()

scenario.ktc = scenario.create_ktc(choices.extensionparameters.ktc_low.value,
                                   choices.extensionparameters.ktc_high.value,
                                   choices.extensionparameters.ktc_limit.value,
                                   not choices.extensions.create_ktc_map.value)

scenario.ktc.plot(nodata=-9999)

scenario.prepare_input_files()

scenario.create_ini_file()

scenario.run_model(watem_sedem_binary)

Run model with grass strips

scenario_nr = 2
scenario = Scenario(catchment, 2019, scenario_nr, choices)
scenario.vct_parcels = inputdata_folder / r"parcels.shp"

Assign grass strips

scenario.vct_grass_strips = inputdata_folder / r"grass_strips.shp"

scenario.vct_grass_strips.plot(column="width")

scenario.composite_landuse = scenario.create_composite_landuse()

scenario.composite_landuse.plot()

scenario.cfactor = scenario.create_cfactor()

scenario.ktc = scenario.create_ktc(choices.extensionparameters.ktc_low.value,
                                   choices.extensionparameters.ktc_high.value,
                                   choices.extensionparameters.ktc_limit.value,
                                   not choices.extensions.create_ktc_map.value)

Prepare run and execute

scenario.prepare_input_files()
scenario.create_ini_file()
scenario.run_model(watem_sedem_binary)

Run model with buffers

scenario_nr = 3
scenario = Scenario(catchment, 2019, scenario_nr, choices)
scenario.vct_parcels = inputdata_folder / r"parcels.shp"

scenario.choices.extensions.include_buffers = True

scenario.vct_buffers = inputdata_folder / r"buffers.shp"

scenario.composite_landuse = scenario.create_composite_landuse()
scenario.cfactor = scenario.create_cfactor()
scenario.ktc = scenario.create_ktc(choices.extensionparameters.ktc_low.value,
                                   choices.extensionparameters.ktc_high.value,
                                   choices.extensionparameters.ktc_limit.value,
                                   not choices.extensions.create_ktc_map.value)

scenario.prepare_input_files()
scenario.create_ini_file()
# scenario.run_model(watem_sedem_binary)

scenario.buffers.plot(nodata=0)

Run model with technical tillage measures

scenario_nr = 4
scenario = Scenario(catchment, 2019, scenario_nr, choices)
scenario.vct_parcels = inputdata_folder / r"parcels.shp"

Technical tillage measures are implemented at the level of a parcel, for which we can define a “reduction”. The column ‘C_reduct’ is used to reduce the final C-factor:

\[C_{factor,reduced}=C_{factor}∗(1-C_{reduction})\]

In this example we add a reduction of 80% to all parcels with LANDUSE code -9999.

# We take all parcels with a crop code of 311 to reduce with 80 %.
scenario.vct_parcels.geodata.loc[scenario.vct_parcels.geodata["LANDUSE"]==-9999, "C_reduct"] = 0.8

Enable source-oriented measures by setting ‘UseTeelttechn’ to one.

scenario.composite_landuse = scenario.create_composite_landuse()
scenario.cfactor = scenario.create_cfactor(use_source_oriented_measures=True)
scenario.ktc = scenario.create_ktc(choices.extensionparameters.ktc_low.value,
                                   choices.extensionparameters.ktc_high.value,
                                   choices.extensionparameters.ktc_limit.value,
                                   not choices.extensions.create_ktc_map.value)

scenario.prepare_input_files()
scenario.create_ini_file()
# scenario.run_model(watem_sedem_binary)

Run model with forced routing

scenario_nr = 5
scenario = Scenario(catchment, 2019, scenario_nr, choices)
scenario.vct_parcels = inputdata_folder / r"parcels.shp"

scenario.choices.extensions.force_routing = True

scenario.vct_force_routing = inputdata_folder / "force_routing.shp"

scenario.composite_landuse = scenario.create_composite_landuse()
scenario.cfactor = scenario.create_cfactor()
scenario.ktc = scenario.create_ktc(choices.extensionparameters.ktc_low.value,
                                   choices.extensionparameters.ktc_high.value,
                                   choices.extensionparameters.ktc_limit.value,
                                   not choices.extensions.create_ktc_map.value)

scenario.prepare_input_files()
scenario.create_ini_file()
scenario.run_model(watem_sedem_binary)

Run model with specific outlets

scenario_nr = 6
scenario = Scenario(catchment, 2019, scenario_nr, choices)
scenario.vct_parcels = inputdata_folder / r"parcels.shp"

scenario.choices.extensions.manual_outlet_selection = True

scenario.vct_outlets = inputdata_folder / "outlets.shp"

scenario.composite_landuse = scenario.create_composite_landuse()
scenario.cfactor = scenario.create_cfactor()
scenario.ktc = scenario.create_ktc(choices.extensionparameters.ktc_low.value,
                                   choices.extensionparameters.ktc_high.value,
                                   choices.extensionparameters.ktc_limit.value,
                                   not choices.extensions.create_ktc_map.value)

scenario.prepare_input_files()
scenario.create_ini_file()
# scenario.run_model(watem_sedem_binary)

Run model without parcels

scenario_nr = 7
scenario = Scenario(catchment, 2019, scenario_nr, choices)

scenario.composite_landuse = scenario.create_composite_landuse()
scenario.cfactor = scenario.create_cfactor()
scenario.ktc = scenario.create_ktc(choices.extensionparameters.ktc_low.value,
                                   choices.extensionparameters.ktc_high.value,
                                   choices.extensionparameters.ktc_limit.value,
                                   not choices.extensions.create_ktc_map.value)

scenario.prepare_input_files()
scenario.create_ini_file()
# scenario.run_model(watem_sedem_binary)

scenario.composite_landuse.plot()