#!/bin/bash

##############################################################

# Script for generating coastal zones based on 100m CLC data
# Hermann, December 2011

##############################################################

# 0. Preparatory settings

# sfx=523				# Andrus, option A
# sfx=52x				# Andrus, option B

# sfx=523_1_to_10km		# Hermann, option A for EEA SDI
sfx=52x_1_to_10km		# Hermann, option B for EEA SDI

# Before generating Option B: reclassify input layers as follows
# for map in FK90@leac FK00@leac clc06@clcv16 ; do
#	r.reclass $map out=${map%@*} rules=reclass_52x.rules
# done
#
# cat reclass_52x.rules
# 521 = 523
# 522 = 523
# * = *

# For debugging
set -x

# See how long it takes
echo $(date)

# Check if GRASS is running
if [ -z "$GISBASE" ] ; then
	echo "You must be in GRASS GIS to run this program." >&2
	exit 1
fi

# Resolution
res=100

##############################################################

# 1. Generation of the buffer zones

# Zoom to FK00 region = 46000 rows x 59000 cols
g.region -p rast=FK00 res=$res

# A Divine task: one has to devide the Earth into LAND and SEA
# After having done this, one can safely throw the land part away

# Before buffering the SEA, the max. area of the SEA is determined across
# FK90/FK00/CLC06, then masked with the max. LAND area of FK90/FK00/CLC06

# Buffer remainig SEA in appropriate distances. To do this in 100m 
# resolution, you need to divide the region in overlapping windows 
# of ~ 1000 km each. Approach here: windows of 1040 km are used.
# Later, the Northern 20 km and the Southern 20 km are clipped
# in order to avoid buffer artefacts at window (= patch) borders.
# This takes about 3 minutes per window (on whitefish, CPU bound).

for s in 4500000 3500000 2500000 1500000 500000 ; do

	# Remove existing MASK
	r.mask -r

	# Set the region for this window
	echo && g.region -p s=$((  $s - 20000 ))  n=$(( $s + 1020000 ))

	# Find the minimum NULL area, used later for clipping the "wrong" side of the buffer
	time r.mapcalc minarea.NULL.$s = "FK90 == 990 && FK00 == 990 && clc06 == 990 ? 1 : 0"

	# max SEA minus max LAND = new SEA
	time r.mapcalc maxarea.lt.$sfx = "FK90  < 523 || FK00  < 523 || clc06  < 523 ? 1 : 0"
	time r.mapcalc maxarea.of.$sfx = "FK90 == 523 || FK00 == 523 || clc06 == 523 ? 1 : 0"
	time r.mapcalc newarea.of.$sfx = "maxarea.of.$sfx == 1 && maxarea.lt.$sfx == 0 ? 1 : 0"

	# We now buffer the new SEA
	r.mask -o newarea.of.$sfx maskcats=1

	echo && time r.buffer --o in=newarea.of.$sfx out=tmp1.$s \
			distances=1000,2000,3000,4000,5000,6000,7000,8000,9000,10000 units=meters
done

# Remove the raster <MASK>
r.mask -r

# Temporary patches generated in this step
# tmp1.4500000    Distance Zones = Most Northern window
# tmp1.3500000    Distance Zones = Scandinavia, UK, etc.
# tmp1.2500000    Distance Zones = France, Germany, etc.
# tmp1.1500000    Distance Zones = Mediterranean region
# tmp1.500000     Distance Zones = Canary Islands, etc.

##############################################################

# Use r.mapcalc to copy the central part of the patches into
# non-overlapping new patches by changing the region accordingly
# This intermediate step takes about 40 seconds per window

for s in 4500000 3500000 2500000 1500000 500000 ; do
  echo && g.region -p s=$s  n=$(( $s + 1000000 ))
  echo && time r.mapcalc tmp2.$s = tmp1.$s
done

# Temporary patches generated in this step
# tmp2.4500000    Distance Zones = Most Northern window
# tmp2.3500000    Distance Zones = Scandinavia, UK, etc.
# tmp2.2500000    Distance Zones = France, Germany, etc.
# tmp2.1500000    Distance Zones = Mediterranean region
# tmp2.500000     Distance Zones = Canary Islands, etc.

##############################################################

# Patch these non-overlapping pieces into a new layer.
# This takes ~ 10 minutes. Get back to the region first.

g.region -p rast=FK00 res=$res

time r.patch --o tmp2.500000,tmp2.1500000,tmp2.2500000,tmp2.3500000,tmp2.4500000 out=tmp2.patched

# Patch the min NULL area together. Patches are overlapping but this does not matter here.
time r.patch --o minarea.NULL.500000,minarea.NULL.1500000,minarea.NULL.2500000,minarea.NULL.3500000,minarea.NULL.4500000 out=minarea.NULL

##############################################################

# In the above steps, the water has grown into 2 directions:
# towards LAND (this is what we want) and towards the SEA, (not wanted).
# The min NULL area is now used to clip the "wrong" buffer side away.

time r.mapcalc buffer_$sfx = "minarea.NULL == 1 || tmp2.patched == 1 ? null() : tmp2.patched - 1"

##############################################################

# Map buffer_$fx now only contains the LAND side of the buffer.
# Its categories are:
#  1 =    0 -  1000m distance
#  2 = 1000 -  2000m distance
# ...
# 10 = 9000 - 10000m distance

##############################################################

# 2. Clean up, make some cosmetics and statistics

# Set the colour table: 1 = red, 2 = magenta, ...
# r.colors buffer_$sfx rules=buffer.colors
r.colors buffer_$sfx color=random

# Don't forget to remove temporary files
g.remove $(g.mlist pat=*00000 sep=,)
g.remove $(g.mlist pat=*area* sep=,)

# Do some statistics
time r.stats --q -c buffer_$sfx > buffer_$sfx.stats

g.region rast=FK00
r.mask -r

r.out.gdal buffer_$sfx out=buffer_$sfx.tif type=Byte nodata=0 createopt=compress=deflate,tiled=yes

echo $(date)

exit 0
