Lisem scripts

kampala.sh

#!/bin/bash -e

./scripts/lisem_mask_prepare.sh --resolution 20 DEM/kampala_demARCs1.tif 19887 --cleanup

./scripts/lisem_resample.sh landuse2016/frBU2016.tif mask20m.map build20m.map
./scripts/lisem_resample.sh landuse2016/frSOIL2016.tif mask20m.map bare20m.map
./scripts/lisem_resample.sh landuse2016/frVEG2016.tif mask20m.map veg20m.map


./scripts/lisem_resample.sh soils/bulkdens_sl2r.tif mask20m.map bulkdens_sl2r.map
./scripts/lisem_resample.sh soils/clay_sl2r.tif mask20m.map clay_sl2r.map
./scripts/lisem_resample.sh soils/orgcarb_sl2r.tif mask20m.map orgcarb_sl2r.map
./scripts/lisem_resample.sh soils/sand_sl2r.tif mask20m.map sand_sl2r.map

./scripts/lisem_resample.sh soils/bulkdens_sl5r.tif mask20m.map bulkdens_sl5r.map
./scripts/lisem_resample.sh soils/clay_sl5r.tif mask20m.map clay_sl5r.map
./scripts/lisem_resample.sh soils/orgcarb_sl5r.tif mask20m.map orgcarb_sl5r.map
./scripts/lisem_resample.sh soils/sand_sl5r.tif mask20m.map sand_sl5r.map

lisem_mask_prepare.sh

#!/bin/bash -e

set -o errexit -o pipefail -o noclobber -o nounset

! getopt --test > /dev/null
if [[ ${PIPESTATUS[0]} -ne 4 ]]; then
    echo "I’m sorry, `getopt --test` failed in this environment."
    exit 1
fi

# A POSIX variable
# Reset in case getopts has been used previously in the shell.
SHORTOPTS=l:r:ch
LONGOPTS=ldd-value:,resolution:,cleanup,help

# Initialize our own variables:
INPUT_MAP_DEM=""
DISCHARGE_POINT=0
LDD_VALUE=1e9
OUTPUT_RESOLUTION=20
SHOULD_CLEANUP=0

function myHelp () {
# Using a here doc with standard out.
cat <<-END
Tool for preparing DEM and mask data.

USAGE:

    ./lisem_mask_prepare.sh DEM_FILE.tif [DISCHARGE_POINT] [PARAMS]

EXAMPLE:
./lisem_mask_prepare.sh --resolution 20 DEM/kampala_demARCs1.tif 19744 --cleanup

ARGUMENTS:

    DEM_FILE.tif - should be raster file with DEM data, currently only tif supported
    DISCHARGE_POINT - if DISCHARGE_POINT not provided, then the script stops before creating watershed.
    -c|--cleanup - delete all the temp files and keep only the resampled DEM and mask
    -r|--resolution - desired resolution. Please make sure the value is multiple of input resolution of the DEM.
    -l|--ldd-value - nobody cares...
END
}

# -use ! and PIPESTATUS to get exit code with errexit set
# -temporarily store output to be able to check for errors
# -activate quoting/enhanced mode (e.g. by writing out “--options”)
# -pass arguments only via   -- "$@"   to separate them correctly
! PARSED=$(getopt --options=$SHORTOPTS --longoptions=$LONGOPTS --name "$0" -- "$@")
if [[ ${PIPESTATUS[0]} -ne 0 ]]; then
    # e.g. return value is 1
    #  then getopt has complained about wrong arguments to stdout
    exit 2
fi

# read getopt’s output this way to handle the quoting right:
eval set -- "$PARSED"

# now enjoy the options in order and nicely split until we see --
while true; do
    case "$1" in
        -l|--ldd-value)
            LDD_VALUE="$2"
            shift 2
            ;;
        -r|--resolution)
            OUTPUT_RESOLUTION="$2"
            shift 2
            ;;
        -c|--cleanup)
            SHOULD_CLEANUP=1
            shift
            ;;
        -h|--help)
            myHelp
            exit 0
            ;;
        --)
            shift
            break
            ;;
        *)
            echo "Programming error"
            exit 3
            ;;
    esac
done

# handle non-option arguments
if [[ $# -le 2 && $# -ge 1 ]]; then
  INPUT_MAP_DEM=$1

  if [[ $# -eq 2 ]]; then
    DISCHARGE_POINT=$2
  fi
else
  echo "$0: A single input DEM file and discharge point value are required."
  exit 4
fi


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

# Selecting a catchment is done in PCRaster, so we have to convert this file
# to PCRaster format (from .tif to .map), type in (GDAL is integrated in PCRaster to facilitate conversion):
pcrcalc "demarc.map = \"$INPUT_MAP_DEM\""

# Get the current grid cell size
CURRENT_RESOLUTION=$(mapattr -p demarc.map  | grep cell_length | tr -d -c 0-9)

if (( OUTPUT_RESOLUTION % CURRENT_RESOLUTION )); then
  echo "Unable to resample, curent resolution is $CURRENT_RESOLUTION, output resolution set to $OUTPUT_RESOLUTION. Must be divisable without a remainder"
  exit 1
fi

RESAMPLE_ARG=-r$(($OUTPUT_RESOLUTION / $CURRENT_RESOLUTION))
OUTPUT_MAP_DEM="demarc"$OUTPUT_RESOLUTION"m.map"
OUTPUT_MAP_MASK="mask"$OUTPUT_RESOLUTION"m.map"

# We will resample this to a lower resolution, from 5m to 20m, for the modelling, this takes some time!
resample demarc.map $OUTPUT_MAP_DEM $RESAMPLE_ARG

# We also have to force the coordinate system from “bottom to top” because of a bug in PCRaster:
mapattr -s $OUTPUT_MAP_DEM -P yb2t

# To determine the boundaries of a watershed, we first create a flow network, using the steepest
# slope. In PCRaster this is called a Local Drain Direction network. The command to create a LDD
# network is:
pcrcalc "ldd.map = lddcreate( $OUTPUT_MAP_DEM, $LDD_VALUE, $LDD_VALUE, $LDD_VALUE, $LDD_VALUE )"

# We will now find a location that is a likely outlet for the catchment under investigation. First make
# a cumulative flow network. This map has the value ‘1’ at the tops mof hills where the flow starts,
# and each time you go a cell downslope a value ‘1’ is added. In the valley at the end of a network,
# the value is the sum of all cells upstream.
pcrcalc "ups.map = accuflux(ldd.map, 1)"

if [[ $DISCHARGE_POINT = 0 ]]; then
  echo "ldd.map created, but no DISCHARGE_POINT provided"
  exit 0
fi

# Check the discharge point in the end of your study area in ups.map and set it as a DISCHARGE_POINT value.
# It means that X upstream cells are connected to this cell.Creating the watershed (or catchment)
# for this point is done as follows:
pcrcalc "ws.map = catchment(ldd.map, if(ups.map eq $DISCHARGE_POINT then nominal(1)))"

# and to select the watershed only and fill the surrounding are with missing value:
pcrcalc "ws1.map = if(ws.map eq 1 then scalar(1))"

# finally we cut the research area down to size, removing all missing values to create the mask for the database:
resample ws1.map $OUTPUT_MAP_MASK -C

# and we cut the dem down to size:
resample --clone $OUTPUT_MAP_MASK $OUTPUT_MAP_DEM $OUTPUT_MAP_DEM

# and multiply the $OUTPUT_MAP_DEM with the mask to see only the catchment:
pcrcalc "$OUTPUT_MAP_DEM *= $OUTPUT_MAP_MASK"

if (($SHOULD_CLEANUP)); then
  echo "Cleaning up..."
  rm ws1.map ws.map ldd.map ups.map demarc.map
fi

lisem_resample.sh

#!/bin/bash -e

INPUT_FILENAME=$1
MASK_FILENAME=$2
OUTPUT_FILENAME=$3

INPUT_FILENAME_TMP=$(basename -- "$1")
INPUT_FILENAME_EXTENSION="${INPUT_FILENAME_TMP##*.}"
FILENAME="${INPUT_FILENAME_TMP%.*}"


pcrcalc "$FILENAME.map = \"$INPUT_FILENAME\""
mapattr -s "$FILENAME.map" -P yb2t
resample --clone $MASK_FILENAME "$FILENAME.map" $OUTPUT_FILENAME

if ! [[ "$FILENAME.map" == $OUTPUT_FILENAME ]]; then
  rm "$FILENAME.map"
fi