g.extension
to install the following GRASS GIS ADDONS
:
%%file gutil.ipy
import base64
from IPython.display import HTML, display
def makefigure(layers, label=''):
html = "<table><tr>"
for i in layers:
outname=i+'.png'
rasterlist = !g.list raster
if i in rasterlist:
!r.out.png -t -w input={i} output={i}.png --o --qq
else:
!v.out.png input={i} output={i}.png --o --qq
imageFile = open(outname, "rb")
imagebyte = base64.b64encode(imageFile.read())
imagestr = imagebyte.decode()
html+="""<td><p><img src="data:image/png;base64,%s" alt="" width='600'/></p></td>""" % imagestr
html += "</tr></table><p><center>%s</center></p>" % label
display(HTML(html))
Overwriting gutil.ipy
%run gutil.ipy
!wget -O bathy_2015.tif https://nextcloud.epinux.com/index.php/s/4bgMkM8pgtLLNxQ/download -q --show-progress
bathy_2015.tif 100%[===================>] 17.91M 6.39MB/s in 2.8s
!r.in.gdal input=bathy_2015.tif \
output=bathy_2015 -o -e --o --qq
!g.region raster=bathy_2015 -ap
projection: 1 (UTM) zone: 19 datum: wgs84 ellipsoid: wgs84 north: 4546448.49987433 south: 4545769.24228443 west: 505691.47804011 east: 509148.78913096 nsres: 1.00037937 ewres: 1.00037937 rows: 679 cols: 3456 cells: 2346624
!r.colors color=haxby \
map=bathy_2015 -e --q
GRASS GIS
working region to the extent of the control unit #1 (visually identified sand dune)north: 4546380
south: 4546202
west: 506185
east: 506571
!g.region n=4546380 \
s=4546202 \
w=506185 \
e=506571 -apl
projection: 1 (UTM) zone: 19 datum: wgs84 ellipsoid: wgs84 north: 4546380 south: 4546202 west: 506185 east: 506571 nsres: 1 ewres: 1 rows: 178 cols: 386 cells: 68708 north-west corner: long: 68:55:34.973432W lat: 41:04:07.123165N north-east corner: long: 68:55:18.433388W lat: 41:04:07.112269N south-east corner: long: 68:55:18.440228W lat: 41:04:01.339863N south-west corner: long: 68:55:34.97987W lat: 41:04:01.350758N center longitude: 68:55:26.70673W center latitude: 41:04:04.231514N
#thin crest detector
!r.geomorphon elevation=bathy_2015 \
forms=geomorphometry \
search=9 \
skip=3 \
flat=2.0 \
dist=0 \
step=0 \
start=0 --o --qq
!r.mapcalc expression="swc=if(geomorphometry==3 | geomorphometry==2,1, null())" --o --q
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:
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).
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).
!r.thin input=swc \
output=swc_thin \
iterations=400 --o --q
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).
!r.clump input=swc_thin \
output=swc_thin_clump -d --o --q
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).
!r.area input=swc_thin_clump \
output=swc_thin_long \
lesser=70 --o --q
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).
!r.to.vect input=swc_thin_long \
output=swc_thin_long_v \
type=line --o --qq
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).
!v.generalize input=swc_thin_long_v \
output=swc_thin_long_smooth_v \
method=douglas \
threshold=2 --o --qq
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).
!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
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).
!g.region n=4546380 \
s=4546202 \
w=506185 \
e=506571 -ap
projection: 1 (UTM) zone: 19 datum: wgs84 ellipsoid: wgs84 north: 4546380 south: 4546202 west: 506185 east: 506571 nsres: 1 ewres: 1 rows: 178 cols: 386 cells: 68708
!r.geomorphon elevation=bathy_2015 \
forms=geomorphometry_sw \
search=30 \
skip=9 \
flat=3.7 \
dist=15 \
step=0 \
start=0 --o --qq
!r.mapcalc \
expression="sw=if(geomorphometry_sw==3 | geomorphometry_sw==2 | geomorphometry_sw==6 | geomorphometry_sw==5 , 1, null())" \
--o --qq
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).
!r.clump input=sw \
output=sw_clump -d --o --qq
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).
!r.area input=sw_clump \
output=sw_clean \
lesser=1000 --o --qq
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).
!r.neighbors input=sw_clean \
output=sw_clean_fill3 \
method=maximum --o --qq
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).
!r.to.vect input=sw_clean_fill3 \
output=sw_clean_fill3_v \
type=area --o --qq
!g.region n=4546380 \
s=4546202 \
w=506185 \
e=506571 -ap
projection: 1 (UTM) zone: 19 datum: wgs84 ellipsoid: wgs84 north: 4546380 south: 4546202 west: 506185 east: 506571 nsres: 1 ewres: 1 rows: 178 cols: 386 cells: 68708
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).
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).
!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
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).
!v.buffer input=sw_ridges_v \
output=sw_buffer_v \
distance=1 --o --qq
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).
!v.overlay ainput=sw_clean_fill3_v \
binput=sw_buffer_v \
operator=not \
output=sw_splitted_v --o --qq
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).
!v.split input=swc_thin_long_smooth_clean_v \
output=ridges_smooth_split_v \
length=70 --o --qq
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).
!v.transects input=ridges_smooth_split_v \
output=ridges_smooth_transect_v \
transect_spacing=10 \
dleft=70 \
dright=70 --o --qq
WARNING: Vector map <ridges_smooth_transect_v> already exists and will be overwritten
makefigure(['sw_splitted_v',
'ridges_smooth_transect_v'],
label="<center><b>Figure 24:</b> identification of dune's side.</center>")