%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 geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pathlib import Path
import os

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

import sys

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()

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()

catchment.river.plot()

You can also have a look at the extracted segment

catchment.segments.plot()

Assign and plot infrastructure polygons and lines

catchment.vct_infrastructure_buildings = vct_infra_poly
catchment.infrastructure_buildings.plot()

catchment.vct_infrastructure_roads = vct_infra_line
catchment.infrastructure_roads.plot()

catchment.infrastructure.plot()

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

The UserChoices-object

The second element we have to defines is the user choices, variables and other paramters. 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 a UserChoices-object

from pywatemsedem.userchoices import UserChoices

choices = UserChoices()

Make use of a preset default values for choices

default_choices = inputdata_folder / "userchoices.ini"

default_choices.exists()

choices.set_ecm_options(default_choices)
choices.dict_ecm_options

choices.set_model_version("WS")

choices.set_model_options(default_choices)

choices.set_model_variables(default_choices)
choices.dict_variables

choices.set_output(default_choices)
choices.dict_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 userchoices-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(bool(choices.dict_ecm_options["UseTeelttechn"]))

scenario.cfactor.plot()

scenario.ktc = scenario.create_ktc(choices.dict_variables["ktc low"],
                                   choices.dict_variables["ktc high"],
                                   choices.dict_variables["ktc limit"],
                                   choices.dict_model_options["UserProvidedKTC"])

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.choices.dict_ecm_options["UseGras"] = 1

scenario.composite_landuse = scenario.create_composite_landuse()

scenario.composite_landuse.plot()

scenario.cfactor = scenario.create_cfactor(bool(choices.dict_ecm_options["UseTeelttechn"]))

scenario.ktc = scenario.create_ktc(choices.dict_variables["ktc low"],
                                   choices.dict_variables["ktc high"],
                                   choices.dict_variables["ktc limit"],
                                   choices.dict_model_options["UserProvidedKTC"])

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.vct_buffers = inputdata_folder / r"buffers.shp"

scenario.choices.dict_ecm_options["Include buffers"] = 1

scenario.composite_landuse = scenario.create_composite_landuse()
scenario.cfactor = scenario.create_cfactor(bool(choices.dict_ecm_options["UseTeelttechn"]))
scenario.ktc = scenario.create_ktc(choices.dict_variables["ktc low"],
                                   choices.dict_variables["ktc high"],
                                   choices.dict_variables["ktc limit"],
                                   choices.dict_model_options["UserProvidedKTC"])

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

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})\]

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

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

scenario.choices.dict_ecm_options['UseTeelttechn'] = 1

scenario.composite_landuse = scenario.create_composite_landuse()
scenario.cfactor = scenario.create_cfactor(bool(choices.dict_ecm_options["UseTeelttechn"]))
scenario.ktc = scenario.create_ktc(choices.dict_variables["ktc low"],
                                   choices.dict_variables["ktc high"],
                                   choices.dict_variables["ktc limit"],
                                   choices.dict_model_options["UserProvidedKTC"])

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

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.dict_model_options["Force Routing"] = 1

scenario.vct_force_routing = inputdata_folder / "force_routing.shp"

scenario.composite_landuse = scenario.create_composite_landuse()
scenario.cfactor = scenario.create_cfactor(bool(choices.dict_ecm_options["UseTeelttechn"]))
scenario.ktc = scenario.create_ktc(choices.dict_variables["ktc low"],
                                   choices.dict_variables["ktc high"],
                                   choices.dict_variables["ktc limit"],
                                   choices.dict_model_options["UserProvidedKTC"])

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

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.dict_model_options

scenario.choices.dict_model_options["Manual outlet selection"] = 1

scenario.vct_outlets = inputdata_folder / "outlets.shp"

scenario.composite_landuse = scenario.create_composite_landuse()
scenario.cfactor = scenario.create_cfactor(bool(choices.dict_ecm_options["UseTeelttechn"]))
scenario.ktc = scenario.create_ktc(choices.dict_variables["ktc low"],
                                   choices.dict_variables["ktc high"],
                                   choices.dict_variables["ktc limit"],
                                   choices.dict_model_options["UserProvidedKTC"])

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(bool(choices.dict_ecm_options["UseTeelttechn"]))
scenario.ktc = scenario.create_ktc(choices.dict_variables["ktc low"],
                                   choices.dict_variables["ktc high"],
                                   choices.dict_variables["ktc limit"],
                                   choices.dict_model_options["UserProvidedKTC"])

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

from pywatemsedem.geo.rasters import AbstractRaster

d = AbstractRaster(None,None)