#!/usr/bin/env python
# coding: utf-8

# * **Use the module [`g.extension`]() to install the following `GRASS GIS ADDONS`:**
#     * [`r.geomorphon`](https://grass.osgeo.org/grass72/manuals/addons/r.geomorphon.html)
#     * [`r.area`](https://grass.osgeo.org/grass72/manuals/addons/r.area.html)
#     * [`v.transects`](https://grass.osgeo.org/grass72/manuals/addons/v.transects.html)

# In[1]:


get_ipython().run_cell_magic('file', 'gutil.ipy', 'import base64\nfrom IPython.display import HTML, display\ndef makefigure(layers, label=\'\'):\n    html = "<table><tr>"\n    for i in layers:\n        outname=i+\'.png\'\n        rasterlist = !g.list raster\n        if i in rasterlist:\n            !r.out.png -t -w input={i} output={i}.png --o --qq\n        else:\n            !v.out.png input={i} output={i}.png --o --qq\n        imageFile = open(outname, "rb")\n        imagebyte = base64.b64encode(imageFile.read())\n        imagestr = imagebyte.decode()\n        html+="""<td><p><img src="data:image/png;base64,%s" alt="" width=\'600\'/></p></td>""" % imagestr\n    html += "</tr></table><p><center>%s</center></p>" % label \n    display(HTML(html))\n')


# In[2]:


get_ipython().run_line_magic('run', 'gutil.ipy')


# * **Get the bathymetry data**

# In[3]:


get_ipython().system('wget -O bathy_2015.tif https://nextcloud.epinux.com/index.php/s/4bgMkM8pgtLLNxQ/download -q --show-progress')


# * **Import the data and set**
# 
#     * `GRASS GIS` command: [`r.in.gdal`]()
#     
# * **Set the colortable to `haxby` with histogram equalization**
# 
#     * `GRASS GIS` command: [`r.colors`]()

# In[4]:


get_ipython().system('r.in.gdal input=bathy_2015.tif             output=bathy_2015 -o -e --o --qq')


# In[5]:


get_ipython().system('g.region raster=bathy_2015 -ap')


# In[6]:


get_ipython().system('r.colors color=haxby            map=bathy_2015 -e --q')


# * **Set the `GRASS GIS` working region to the extent of the control unit #1 (visually identified sand dune)**
# 
# ```
# north:      4546380
# south:      4546202
# west:       506185
# east:       506571
# ```

# In[7]:


get_ipython().system('g.region n=4546380            s=4546202            w=506185            e=506571 -apl')


# In[8]:


#thin crest detector
get_ipython().system('r.geomorphon elevation=bathy_2015                forms=geomorphometry                search=9                skip=3                flat=2.0                dist=0                step=0                start=0 --o --qq')
get_ipython().system('r.mapcalc expression="swc=if(geomorphometry==3 | geomorphometry==2,1, null())" --o --q')


# In[9]:


makefigure(['bathy_2015', 
            'geomorphometry'], 
           label='Output of r.geomorphon, input parameters: search=9 skip=3 flat=2.0 dist=0 step=0 start=0 ')


# The whole GRM workflow can be divided into three sections:
# 1. Extract and vectorize sand wave crest (SWC).
# 2. Extract and vectorize sand wave main body (SW).
# 3. Identify lee and stoss side and compute sand wave’s metrics.
# 
# The sequence of each step are summarized and visually described below:
# 
# `1. Extract and vectorize sand wave crest for each bedform`
# 
#     1.1 Extraction of sand wave crest (SWC): From the TFE results (r.geomorphon), SWC are identified by extracted by reclassifying the cells with feature type equal to ridge and summit (feature type class: 2,3) into a new raster feature with category value 1 and setting the remaining cells to null (Figure 7).

# In[10]:


makefigure(['geomorphometry', 
            'swc'], 
           label='<b>Figure 7:</b> Extraction of sand wave crest (SWC).')


#     1.2 SWC thinning: SWC areas are thinned and reduced to a single pixel width (Figure 8).

# In[11]:


get_ipython().system('r.thin input=swc          output=swc_thin          iterations=400 --o --q')


# In[12]:


makefigure(['swc', 
            'swc_thin'], 
           label='<b>Figure 8:</b> SWC thinning.')


#     1.3 SWC clumping: Each SWC is recategorized by grouping cells that form physically discrete areas into unique categories and assign a distinct color to each raster feature, different colors are assigned to each linear feature (Figure 9).

# In[13]:


get_ipython().system('r.clump input=swc_thin           output=swc_thin_clump -d --o --q')


# In[14]:


makefigure(['swc_thin', 
            'swc_thin_clump'], 
           label='<b>Figure 9:</b> SWC clumping.')


#     1.4 SWC filtering by length: Each feature with same category shorter than a given threshold is removed (Figure 10).
# 

# In[15]:


get_ipython().system('r.area input=swc_thin_clump          output=swc_thin_long          lesser=70 --o --q')


# In[16]:


makefigure(['swc_thin_clump', 
            'swc_thin_long'], 
           label='<b>Figure 10:</b> SWC filtering by length.')


#     1.5 SWC vectorization: Conversion from raster to vector to obtain a vector feature of type line representing an approximate sand wave crest (Figure 11).

# In[17]:


get_ipython().system('r.to.vect input=swc_thin_long             output=swc_thin_long_v             type=line --o --qq')


# In[18]:


makefigure(['swc_thin_long', 
            'swc_thin_long_v'], 
           label='<b>Figure 11:</b> SWC vectorization.')


#     1.6 SWC topological cleaning: Line features are cleaned by removing any dangle (Figure 12).

# In[19]:


get_ipython().system('v.generalize input=swc_thin_long_v                output=swc_thin_long_smooth_v                method=douglas                threshold=2 --o --qq')


# In[20]:


makefigure(['swc_thin_long_v', 
            'swc_thin_long_smooth_v'], 
           label='<b>Figure 12:</b> SWC topological cleaning.')


#     1.7 SWC smoothing: A low pass filter is used to obtain a vectorized Sand Wave Crest (Figure 13).

# In[21]:


get_ipython().system('v.clean input=swc_thin_long_smooth_v           output=swc_thin_long_smooth_clean_v           type=line           tool=rmdangle,rmdangle,rmdangle,rmdangle           threshold=5,10,20,30 --o --qq')


# In[22]:


makefigure(['swc_thin_long_smooth_v', 
            'swc_thin_long_smooth_clean_v'], 
           label='<b>Figure 13:</b> SWC smoothing.')


# `2 Identify areas covered by large scale bedforms`
# 
#     2.1 Sand wave (SW) extraction: From the TFE results (r.geomorphon), the entire landform is extracted by reclassifying the cells with feature type equal to summit, ridge, spur, and slope (feature type: class 2, 3, 5, 6) into a new raster feature with category value 1, and setting the remaining cells to null (Figure 14).

# In[23]:


get_ipython().system('g.region n=4546380            s=4546202            w=506185            e=506571 -ap')


# In[24]:


get_ipython().system('r.geomorphon elevation=bathy_2015                forms=geomorphometry_sw                search=30                skip=9                flat=3.7                dist=15                step=0                start=0 --o --qq')


# In[25]:


get_ipython().system('r.mapcalc      expression="sw=if(geomorphometry_sw==3 | geomorphometry_sw==2 | geomorphometry_sw==6 | geomorphometry_sw==5 , 1, null())"      --o --qq')


# In[26]:


makefigure(['geomorphometry_sw', 
            'sw'], 
           label='<b>Figure 14:</b> Sand wave (SW) extraction.')


#     2.2 SW clumping: The SW raster map is recategorized by grouping cells that form physically discrete areas into unique categories and assign a distinct color to each raster feature (Figure 15).

# In[27]:


get_ipython().system('r.clump input=sw           output=sw_clump -d --o --qq')


# In[28]:


makefigure(['sw', 
            'sw_clump'], 
           label='<b>Figure 15:</b> SW clumping.')


#     2.3 SW filtering by area: Each raster feature is reclassified based on its area, all the features having an area smaller than a given threshold are removed (Figure 16).

# In[29]:


get_ipython().system('r.area input=sw_clump          output=sw_clean          lesser=1000 --o --qq')


# In[30]:


makefigure(['sw_clump', 
            'sw_clean'], 
           label='<b>Figure 16:</b> SWC filtering by area.')


#     2.4 SW filling: Null cells within the discrete areas are filled with the same category of the surrounding pixels (Figure 17).

# In[31]:


get_ipython().system('r.neighbors input=sw_clean               output=sw_clean_fill3               method=maximum --o --qq')


# In[32]:


makefigure(['sw_clean', 
            'sw_clean_fill3'], 
           label='<b>Figure 17:</b> SW filling.')


#     2.5 SW vectorization: Conversion from raster to vector to obtain a vector feature of type polygon representing an approximate sand wave body (Figure 18).

# In[33]:


get_ipython().system('r.to.vect input=sw_clean_fill3             output=sw_clean_fill3_v             type=area --o --qq')


# In[34]:


get_ipython().system('g.region n=4546380            s=4546202            w=506185            e=506571 -ap')


# In[35]:


makefigure(['sw_clean_fill3', 
            'sw_clean_fill3_v'], 
           label='<b>Figure 18:</b> SW vectorization.')


# `3 Identify lee and stoss side and compute sand wave’s metrics`
# 
#     3.1 SW and SWC overlay: The vectorized sand wave crest is overlaid on top of the polygonal area representing the sand wave body (Figure 19).

# In[36]:


makefigure(['sw_clean_fill3_v', 
            'swc_thin_long_smooth_clean_v'], 
           label='<b>Figure 19:</b> SW and SWC overlay.')


#     3.2 SWC clipping: The portion of the sand wave crest that is not included in the sand wave body is removed (Figure 20).

# In[37]:


get_ipython().system('v.overlay ainput=swc_thin_long_smooth_clean_v             atype=line             binput=sw_clean_fill3_v             out=sw_ridges_v             operator=and             olayer=0,1,0 --o --qq')


# In[38]:


makefigure(['sw_clean_fill3_v', 
            'sw_ridges_v'], 
           label='<b>Figure 20:</b> SWC clipping.')


#     3.3 SWC buffering: A polygonal area is created by buffering the sand wave crest (buffer distance equal to the DBM cell size) (Figure 21).

# In[39]:


get_ipython().system('v.buffer input=sw_ridges_v            output=sw_buffer_v            distance=1 --o --qq')


# In[40]:


makefigure(['sw_ridges_v', 
            'sw_buffer_v'], 
           label='<b>Figure 21:</b> SWC buffering.')


#     3.4 SW splitting: The buffered sand wave crest is used to split the sand wave body into two parts (Figure 22).

# In[41]:


get_ipython().system('v.overlay ainput=sw_clean_fill3_v             binput=sw_buffer_v             operator=not             output=sw_splitted_v --o --qq')


# In[42]:


makefigure(['sw_clean_fill3_v', 
            'sw_splitted_v'], 
           label='<b>Figure 22:</b> SW splitting.')


#     3.5 Identification of stoss and lee side: the identification is achieved by performing an univariate statistical analysis on the DBM slope, using each sand wave side as a mask. The side with the higher values of the slope is assigned the label of lee side and is colored in grey, while the side with the lower values of the slope is assigned the label of stoss side (Figure 23).

# In[43]:


get_ipython().system('v.split input=swc_thin_long_smooth_clean_v           output=ridges_smooth_split_v           length=70 --o --qq')


# In[44]:


makefigure(['sw_clean_fill3_v', 
            'sw_splitted_v'], 
           label="<b>Figure 23:</b> identification od dune's side.")


#     3.6 SW height: Derivation of SW height by generating a series of vertical profiles along several transects perpendicular to SWC. The spacing between the equidistant vertical profiles is given by an input parameter and set to 10 m as default value, the length of each profile is also variable by the user (by default is set to be the same length of the sand wave ridge) (Figure 24).

# In[45]:


get_ipython().system('v.transects input=ridges_smooth_split_v               output=ridges_smooth_transect_v               transect_spacing=10               dleft=70               dright=70 --o --qq')


# In[46]:


makefigure(['sw_splitted_v', 
            'ridges_smooth_transect_v'], 
           label="<center><b>Figure 24:</b> identification of dune's side.</center>")