import wflow.wflow_lib as tr
%pylab inline

demfilename = "example_dem.map"

"""
Explanation: 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.
End of explanation
"""


dem = tr.readmap(demfilename)

"""
Explanation: 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.
End of explanation
"""


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

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


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

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

"""
Explanation: 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
End of explanation
"""


!aguila ldd.map

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

"""
Explanation: 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:
End of explanation
"""


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

"""
Explanation: 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:
End of explanation
"""


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

"""
Explanation: 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):
End of explanation
"""


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

"""
Explanation: 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..
End of explanation
"""


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

"""
Explanation: Now use this river network to adjust the DEM. This step (the averaging) takes some time....
End of explanation
"""


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


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

"""
Explanation: Recreate the LDD and see how it looks now...
End of explanation
"""


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

"""
Explanation: Now we read in the gauge location to get a couple of subcatchments
End of explanation
"""


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

"""
Explanation: Instead of defining subcatchments using gauges we can also use river junctions to find them.
End of explanation
"""


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

"""
Explanation: 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
End of explanation
"""


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

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

"""
Explanation: 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:
End of explanation
"""


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

"""
Explanation: Make a wflow_river file by including only rivers with stralerorder of 6 and larger
End of explanation
"""


