Iss4

depth_to_bedrock/README.md

@@ -0,0 +1,48 @@
+# `Global Depth to Bedrock` Geospatial Dataset
+In this file, the necessary technical details of the dataset is explained.
+
+## Location of the `Global Depth to Bedrock` Dataset Files
+The `Global Depth to Bedrock` geospatial dataset files are located under the following directory accessible from Digital Alliance (formerly Compute Canada) Graham cluster:
+
+```console
+/project/rpp-kshook/Model_Output/DTB
+```
+
+And the structure of the files is as following:
+
+```console
+/project/rpp-kshook/Model_Output/DTB
+├── BDRICM_M_10km_ll.zip
+├── BDRICM_M_1km_ll.zip
+├── BDRICM_M_250m_ll.zip
+├── BDRLOG_M_10km_ll.zip
+├── BDRLOG_M_1km_ll.zip
+├── BDRLOG_M_250m_ll.zip
+├── BDTICM_M_10km_ll.zip
+├── BDTICM_M_1km_ll.zip
+└── BDTICM_M_250m_ll.zip
+```
+
+## Spatial and Temporal Extents
+
+The spatial extent of this dataset covers longitudes from `-180` to `+180` degress and latitudes from `-90` to `+90` degress. This dataset is static and does not vary with time. 
+
+## Dataset Variables
+This variables of this dataset are detailed in the table below:
+
+|#	|Variable Name (used in `gistool`)	|Description				|Comments		|
+|-------|---------------------------------------|---------------------------------------|-----------------------|
+|1	|BDRICM_M_10km_ll			|censored DTB[^1] in `cm`		|5-minute resolution	|
+|2	|BDRLOG_M_10km_ll			|occurrence of R horizon in `%`		|5-minute resolution	|
+|3	|BDTICM_M_10km_ll			|absolute DTB in `cm`			|5-minute resolution	|
+|4	|BDRICM_M_1km_ll			|censored DTB in `cm`			|30-second resolution	|
+|5	|BDRLOG_M_1km_ll			|occurrence of R horizon in `%`		|30-second resolution	|
+|6	|BDTICM_M_1km_ll			|absolute DTB in `cm`			|30-second resolution	|
+|7	|BDRICM_M_250m_ll			|censored DTB in `cm`			|7.5-second resolution	|
+|8	|BDRLOG_M_250m_ll			|occurrence of R horizon in `%`		|7.5-second resolution	|
+|9	|BDTICM_M_250m_ll			|absolute DTB in `cm`			|7.5-second resolution	|
+
+[^1]: DTB: Depth to Bedrock
+
+The following [link](http://globalchange.bnu.edu.cn/research/dtbd.jsp) and [paper](http://onlinelibrary.wiley.com/doi/10.1002/2016MS000686/full) provide extensive details of the mentioned variables.
+

depth_to_bedrock/depth_to_bedrock.sh

@@ -0,0 +1,286 @@
+#!/bin/bash
+# GIS Data Processing Workflow
+# Copyright (C) 2022, University of Saskatchewan
+# Copyright (C) 2021, Wouter Knoben
+#
+# This file is part of GIS Data Processing Workflow
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+# =========================
+# Credits and contributions
+# =========================
+# 1. Parts of the code are taken from https://www.shellscript.sh/tips/getopt/index.html
+# 2. General ideas of GeoTIFF subsetting are taken from https://github.com/CH-Earth/CWARHM
+#    developed mainly by Wouter Knoben (hence the header copyright credit). See the preprint
+#    at: https://www.essoar.org/doi/10.1002/essoar.10509195.1
+
+
+# ================
+# General comments
+# ================
+# * All variables are camelCased for distinguishing from function names;
+# * function names are all in lower_case with words seperated by underscore for legibility;
+# * shell style is based on Google Open Source Projects'
+#   Style Guide: https://google.github.io/styleguide/shellguide.html
+
+
+# ===============
+# Usage Functions
+# ===============
+short_usage() {
+  echo "usage: $(basename $0) -cio DIR -v var1[,var2[...]] [-r INT] [-se DATE] [-ln REAL,REAL] [-f PATH] [-t BOOL] [-a stat1[,stat2,[...]] [-q q1[,q2[...]]]] [-p STR] "
+}
+
+
+# argument parsing using getopt - WORKS ONLY ON LINUX BY DEFAULT
+parsedArguments=$(getopt -a -n depth-to-bedrock -o i:o:v:r:s:e:l:n:f:t:a:q:p:c: --long dataset-dir:,output-dir:,variable:,crs:,start-date:,end-date:,lat-lims:,lon-lims:,shape-file:,print-geotiff:,stat:,quantile:,prefix:,cache: -- "$@")
+validArguments=$?
+if [ "$validArguments" != "0" ]; then
+  short_usage;
+  exit 1;
+fi
+
+# check if no options were passed
+if [ $# -eq 0 ]; then
+  echo "$(basename $0): ERROR! arguments missing";
+  exit 1;
+fi
+
+# check long and short options passed
+eval set -- "$parsedArguments"
+while :
+do
+  case "$1" in
+    -i | --dataset-dir)   geotiffDir="$2"      ; shift 2 ;; # required
+    -o | --output-dir)    outputDir="$2"       ; shift 2 ;; # required
+    -v | --variable)      variables="$2"       ; shift 2 ;; # required
+    -r | --crs)		  crs="$2"	       ; shift 2 ;; # required 
+    -s | --start-date)    startDate="$2"       ; shift 2 ;; # redundant - added for compatibility
+    -e | --end-date)      endDate="$2"         ; shift 2 ;; # redundant - added for compatibility
+    -l | --lat-lims)      latLims="$2"         ; shift 2 ;; # required - could be redundant
+    -n | --lon-lims)      lonLims="$2"         ; shift 2 ;; # required - could be redundant
+    -f | --shape-file)    shapefile="$2"       ; shift 2 ;; # required - could be redundant
+    -t | --print-geotiff) printGeotiff="$2"    ; shift 2 ;; # required
+    -a | --stat)	  stats="$2"	       ; shift 2 ;; # optional
+    -q | --quantile)	  quantiles="$2"       ; shift 2 ;; # optional
+    -p | --prefix)	  prefix="$2"          ; shift 2 ;; # optional
+    -c | --cache)	  cache="$2"           ; shift 2 ;; # required
+
+    # -- means the end of the arguments; drop this, and break out of the while loop
+    --) shift; break ;;
+
+    # in case of invalid option
+    *)
+      echo "$(basename $0): ERROR! invalid option '$1'";
+      short_usage; exit 1 ;;
+  esac
+done
+
+# check if $ensemble is provided
+if [[ -n "$startDate" ]] || [[ -n "$endDate" ]]; then
+  echo "$(basename $0): ERROR! redundant argument (time extents) provided";
+  exit 1;
+fi
+
+# check the prefix if not set
+if [[ -z $prefix ]]; then
+  prefix="depth_to_bedrock_"
+fi
+
+# parse comma-delimited variables
+IFS=',' read -ra variables <<< "${variables}"
+
+
+# =====================
+# Necessary Assumptions
+# =====================
+# TZ to be set to UTC to avoid invalid dates due to Daylight Saving
+alias date='TZ=UTC date'
+# expand aliases for the one stated above
+shopt -s expand_aliases
+# hard-coded paths
+renvCache="/project/rpp-kshook/Climate_Forcing_Data/assets/r-envs/" # general renv cache path
+exactextractrCache="${renvCache}/exact-extract-env" # exactextractr renv cache path
+renvPackagePath="${renvCache}/renv_0.15.5.tar.gz" # renv_0.15.5 source path
+
+
+# ==========================
+# Necessary Global Variables
+# ==========================
+# the structure of the original file names are one of the following:
+#     * %{var}_M_sl%{soil_layer_depth_number}_250m_ll.tif
+#     * %{var}_250m_ll.tif
+# but the user is expected to include complete variable (file) name in the 
+# `--variable` input argument comma-separated list.
+
+
+# ===================
+# Necessary Functions
+# ===================
+# Modules below available on Compute Canada (CC) Graham Cluster Server
+load_core_modules () {
+    module -q purge
+    module -q load gcc/9.3.0
+    module -q load r/4.1.2
+    module -q load gdal/3.0.4
+    module -q load udunits/2.2.28
+    module -q load geos/3.10.2
+    module -q load proj/9.0.0
+}
+load_core_modules
+
+
+# =================
+# Useful One-liners
+# =================
+# sorting a comma-delimited string of real numbers
+sort_comma_delimited () { IFS=',' read -ra arr <<< "$*"; echo ${arr[*]} | tr " " "\n" | sort -n | tr "\n" " "; }
+
+# log date format
+logDate () { echo "($(date +"%Y-%m-%d %H:%M:%S")) "; }
+
+
+#######################################
+# subset GeoTIFFs
+#
+# Globals:
+#   latLims: comma-delimited latitude
+#            limits
+#   lonLims: comma-delimited longitude
+#            limits
+#
+# Arguments:
+#   sourceVrt: source vrt file
+#   destPath: destionation path (inclu-
+#	      ding file name)
+#
+# Outputs:
+#   one mosaiced (merged) GeoTIFF under
+#   the $destDir
+#######################################
+subset_geotiff () {
+  # local variables
+  local latMin
+  local latMax
+  local lonMin
+  local lonMax
+  local sortedLats
+  local sortedLons
+  # reading arguments
+  local sourceVrt="$1"
+  local destPath="$2"
+
+  # extracting minimum and maximum of latitude and longitude respectively
+  ## latitude
+  sortedLats=($(sort_comma_delimited "$latLims"))
+  latMin="${sortedLats[0]}"
+  latMax="${sortedLats[1]}"
+  ## longitude
+  sortedLons=($(sort_comma_delimited "$lonLims"))
+  lonMin="${sortedLons[0]}"
+  lonMax="${sortedLons[1]}"
+
+  # subset based on lat/lon - flush to disk at 500MB
+  GDAL_CACHEMAX=500
+  gdal_translate --config GDAL_CACHEMAX 500 \
+  		 -co COMPRESS="DEFLATE" \
+  		 -co BIGTIFF="YES" \
+		 -projwin $lonMin $latMax $lonMax $latMin "${sourceVrt}" "${destPath}" \
+		 > /dev/null;
+}
+
+
+# ===============
+# Data Processing
+# ===============
+# display info
+echo "$(logDate)$(basename $0): processing Depth-to-Bedrock GeoTIFF(s)..."
+
+# make the output directory
+echo "$(logDate)$(basename $0): creating output directory under $outputDir"
+mkdir -p "$outputDir" # making the output directory
+mkdir -p "$cache" # making the cache directory
+
+# if shapefile is provided extract the extents from it
+if [[ -n $shapefile ]]; then
+  # extract the shapefile extent
+  IFS=' ' read -ra shapefileExtents <<< "$(ogrinfo -so -al "$shapefile" | sed 's/[),(]//g' | grep Extent)"
+  # transform the extents in case they are not in EPSG:4326
+  IFS=':' read -ra sourceProj4 <<< "$(gdalsrsinfo $shapefile | grep -e "PROJ.4")" # source Proj4 value
+  if [[ -n $sourceProj4 ]]; then
+    :
+  else
+    echo "$(logDate)$(basename $0): WARNING! Assuming WSG84 CRS for the input ESRI shapefile"
+    sourceProj4=("PROJ.4" " +proj=longlat +datum=WGS84 +no_defs") # made an array for compatibility with the following statements
+  fi
+
+  # transform limits and assing to variables
+  IFS=' ' read -ra leftBottomLims <<< $(echo "${shapefileExtents[@]:1:2}" | gdaltransform -s_srs "${sourceProj4[1]}" -t_srs EPSG:4326 -output_xy)
+  IFS=' ' read -ra rightTopLims <<< $(echo "${shapefileExtents[@]:4:5}" | gdaltransform -s_srs "${sourceProj4[1]}" -t_srs EPSG:4326 -output_xy)
+  # define $latLims and $lonLims from $shapefileExtents
+  lonLims="${leftBottomLims[0]},${rightTopLims[0]}"
+  latLims="${leftBottomLims[1]},${rightTopLims[1]}"
+fi
+
+# subset and produce stats if needed
+if [[ "$printGeotiff" == "true" ]]; then
+  echo "$(logDate)$(basename $0): subsetting GeoTIFFs under $outputDir"
+  for var in "${variables[@]}"; do
+    # subset based on lat and lon values
+    subset_geotiff "${geotiffDir}/${var}.tif" "${outputDir}/${prefix}${var}.tif"
+  done
+fi
+
+## make R renv project directory
+if [[ -n "$shapefile" ]] && [[ -n $stats ]]; then
+  echo "$(logDate)$(basename $0): Extracting stats under $outputDir"
+  mkdir -p "$cache/r-virtual-env/"
+  ## make R renv in $cache
+  virtualEnvPath="$cache/r-virtual-env/"
+  cp "$(dirname $0)/../assets/renv.lock" "$virtualEnvPath"
+  ## load necessary modules - excessive, mainly for clarification
+  load_core_modules
+
+  ## make the temporary directory for installing r packages
+  tempInstallPath="$cache/r-packages"
+  mkdir -p "$tempInstallPath"
+  export R_LIBS_USER="$tempInstallPath"
+
+  # extract given stats for each variable
+  for var in "${variables[@]}"; do
+    ## build renv and create stats
+    Rscript "$(dirname $0)/../assets/stats.R" \
+    	    "$tempInstallPath" \
+  	    "$exactextractrCache" \
+  	    "$renvPackagePath" \
+	    "$virtualEnvPath" \
+	    "$virtualEnvPath" \
+	    "${virtualEnvPath}/renv.lock" \
+	    "${geotiffDir}/${var}.tif" \
+	    "$shapefile" \
+	    "$outputDir/${prefix}stats_${var}.csv" \
+	    "$stats" \
+	    "$quantiles" >> "${outputDir}/${prefix}stats_${var}.log" 2>&1;
+  done
+fi
+
+# remove unnecessary files 
+mkdir "$HOME/empty_dir" 
+echo "$(logDate)$(basename $0): deleting temporary files from $cache"
+rsync --quiet -aP --delete "$HOME/empty_dir/" "$cache"
+rm -r "$cache"
+echo "$(logDate)$(basename $0): temporary files from $cache are removed"
+echo "$(logDate)$(basename $0): results are produced under $outputDir"
+

example/depth-to-bedrock.sh

@@ -0,0 +1,44 @@
+#!/bin/bash
+
+# Geospatial Dataset Processing Workflow
+# Copyright (C) 2022, University of Saskatchewan
+#
+# This file is part of the Geospatial Dataset Processing Workflow
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+# This is a simple example to extract common statistics for the 
+# pfaf 71 - Saskatchewan-Nelson System Shapefiles of MERIT-Basins -
+# from the Depth-to-Bedrock rasters.
+
+# Always call the script in the root directory of this repository
+cd ..
+echo "The current directory is: $(pwd)"
+
+# first download a sample shapefile - could be any shapefile
+wget -m -nd -nv -q -A "cat_pfaf_71_MERIT_Hydro_v07_Basins_v01_bugfix1.*" \
+     "http://hydrology.princeton.edu/data/mpan/MERIT_Basins/MERIT_Hydro_v07_Basins_v01_bugfix1/pfaf_level_02/"; 
+
+# implement subsetting and zonal statistics
+./extract-gis.sh --dataset="depth-to-bedrock" \
+  --dataset-dir="/project/rpp-kshook/Model_Output/DTB/" \
+  --variable="BDRICM_M_10km_ll,BDRLOG_M_10km_ll" \
+  --shape-file="$(pwd)/cat_pfaf_71_MERIT_Hydro_v07_Basins_v01_bugfix1.shp" \
+  --print-geotiff=true \
+  --output-dir="$HOME/scratch/depth-to-bedrock-test/" \
+  --prefix="dtb_" \
+  --stat="mean,min,max" \
+  --email="you@company.ca" \ 
+  -j;
+

extract-gis.sh

@@ -325,7 +325,12 @@ case "${geotiff,,}" in
     call_processing_func "$(dirname $0)/modis/modis.sh"
     ;;
 
-  # Landsar
+  # Depth to Bedrock
+  "depth-to-bedrock" | "depth_to_bedrock" | "dtb" | "depth-to_bedrock" | "depth_to_bedrock")
+    call_processing_func "$(dirname $0)/depth_to_bedrock/depth_to_bedrock.sh"
+    ;;
+
+  # Landsat
   "landsat" | "landast" )
     call_processing_func "$(dirname $0)/landsat/landsat.sh"
     ;;