Introduction

This notebook shows the steps to prepate the maps needed for the wflow_* set of models. The wflow_prepare scripts have a number of specific features that are not used here. This notebook shows the principles and may be used as a template to create better preprocessing options.

In short the steps to generate input files for the wflow models are:

Create a DEM for the catchment
Create a local crainage network file for the catchment used the DEM and additional information
Create a subarea/catchment map by indicating outflow points (gauges) in the map
Convert/resample landuse and soil maps

The folloing maps are required:

wflow_dem.map
wflow_demmax.map (for sbm only)
wflow_demmin.map (for sbm only)
wflow_ldd.map
wflow_landuse.map
wflow_soil.map
wflow_river.map
wflow_subcatch.map

There are more maps but those are optional (but very useful). See the documentation of the wflow_prepare scripts for details

Example

Process DEM to generate the ldd

Begin by importing the wflow_lib module. This module also imports the pcraster framework. Also set the names of the files we are going to use here.



In [2]:

    
import wflow.wflow_lib as tr
%pylab inline

demfilename = "example_dem.map"









    



Populating the interactive namespace from numpy and matplotlib

start by reading the DEM into memory. Note that a DEM can be rather large, in that case it may be good to resample it first before reading it into memory.



In [3]:

    
dem = tr.readmap(demfilename)

Check the dem by plotting it. Note the use of the pcr2numpy function to convert the pcraster map object to a numpy matrix



In [4]:

    
imshow(tr.pcr2numpy(dem,np.nan))
colorbar()









    Out[4]:





<matplotlib.colorbar.Colorbar instance at 0x7f1d940edf80>

Determine the local drainage direction network from the DEM. This determines how the water flows from one cell to another. This may take some time



In [12]:

    
tr.lddcreate("ldd.map",dem,True,

tr.report('ldd.map',ldd)









    



---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-12-e101aae16147> in <module>()
----> 1 tr.lddcreate("ldd.map",dem,True)
      2 
      3 tr.report('ldd.map',ldd)

TypeError: lddcreate() takes exactly 5 arguments (3 given)

We have now created a first version of the LDD. We can plot it with python but unfortunately it does not show the flow directions as nicely as the specialised pcraster graphics (aguila). We can however, use one of the pcraster functions (streamorder) to show how the main rivers run:



In [6]:

    
!aguila ldd.map



In [ ]:

    
figure(figsize = (10,10))
strorder = tr.streamorder(ldd)
imshow(tr.pcr2numpy(strorder,np.nan))

wflow_lib has a function to search for the largest catchment within the current DEM. It returns a map with a single cell filled at the point of the outlet. Next we use that map to determine the whole catchment area:



In [ ]:

    
outletpoint = tr.find_outlet(ldd)
catch = tr.subcatch(ldd,outletpoint)
imshow(tr.pcr2numpy(catch,np.nan))

In this case we get something that is not how we want it to be. In the northern and easter part of the map the catchment delineation seems to be missing stuff as we can be from the map below (which is more-or-less correct). The reason in this case is that a 1x1km DEM is missing some of the feature to correctly deliniate the DEM. One way to correct this is to use a know river course shapefile and use that to "burn in" the river into the DER and re-create the LDD.

This is how it should look (including a subdivision into subcatchments):



In [ ]:

    
goodcatch = tr.readmap("goodcatch.map")
imshow(tr.pcr2numpy(goodcatch,np.nan))

Use a river network to get better results

In this case we use a grid (converted from a shape-file, best to use qgis to do that) to burn-in the river network (only the higher orders). This will hopefully force the catchment deliniation process to follow the rigth track..



In [ ]:

    
river = tr.cover(tr.readmap("river.map"),0)
imshow(tr.pcr2numpy(river,np.nan))

Now use this river network to adjust the DEM. This step (the averaging) takes some time....



In [ ]:

    
tr.setglobaloption("unitcell")
demregional=tr.windowaverage(dem,30)
ldddem = tr.ifthenelse(tr.scalar(river) >= 1.0, demregional -9000 , dem)

Recreate the LDD and see how it looks now...



In [ ]:

    
ldd = tr.lddcreate_save("wflow_ldd.map",ldddem,True)
strorder = tr.streamorder(ldd)



In [ ]:

    
outletpoint = tr.find_outlet(ldd)
catch = tr.subcatch(ldd,outletpoint)
imshow(tr.pcr2numpy(catch,0))

Now we read in the gauge location to get a couple of subcatchments



In [ ]:

    
gauges=tr.readmap("wflow_gauges.map")
sc = tr.subcatch(ldd,gauges)
imshow(tr.pcr2numpy(sc,0.0),vmin=0,vmax=14)

Instead of defining subcatchments using gauges we can also use river junctions to find them.



In [ ]:

    
subcatch_from_order , virtgauges, b = tr.subcatch_order_a(ldd,7)
imshow(tr.pcr2numpy(subcatch_from_order,0.0))

OK, now use the subcatchments to mask the dem to the catchment we use and save the dem. This will speed-up later calculations. Also save the subcatchment and gauges



In [ ]:

    
dem = tr.ifthen(subcatch_from_order >0, dem)
tr.report(dem,"wflow_dem.map")
tr.report(subcatch_from_order,"wflow_subcatch.map")
# Also save the gauged
tr.report(tr.ifthen(virtgauges>0,virtgauges),"wflow_gauges.map")

Now we need to resample the soil and landuse maps to the same resolution. There are multiple ways to do this. You can use gdal or the pcraster resample program. Here we use resample to convert them to the same resolution as the DEM:



In [ ]:

    
!resample --clone wflow_dem.map soil.map wflow_soil.map



In [ ]:

    
!resample --clone wflow_dem.map landuse.map wflow_landuse.map

Make a wflow_river file by including only rivers with stralerorder of 6 and larger



In [ ]:

    
strorder = tr.streamorder(ldd)
tr.report(tr.ifthen(strorder>6,strorder),"wflow_river.map")



In [ ]: