SMESHR - Improve docs

.gitignore

@@ -8,3 +8,7 @@
 
 # Dask
 dask-worker-space
+
+*.ipynb_checkpoints
+*__pycache__
+*archive*

scripts/smeshr/HRRR_linker.sh

@@ -11,7 +11,8 @@
 
 for DAY in $1; do
   pushd $DAY || exit 1
-  DAY=$(echo $DAY | cut -d '/' -f11,12 | sed 's/\//./g' | sed 's/run//g' );
-  ln -fs $2.$DAY.nc $3
+  # get date in YYYYMMDD format
+  DAY=$(echo $DAY | cut -d '/' -f12 | sed 's/\//./g' | sed 's/run//g' );
+  ln -fsv ${2}/*.MST.$DAY.nc $3
   popd
 done

scripts/smeshr/MODIS_albedo_ERW.sh → scripts/smeshr/MODIS-STC/MODIS_albedo_ERW.sh

@@ -9,17 +9,19 @@
 # ${1}: /path/to/MODIS/inputs
 # ${2}: water year
 
-# command from UofU snow-rs package
-# https://github.com/UofU-Cryosphere/snow-rs
-# Converts the MODIS matlab to GeoTiff files
-variable_from_modis --source-dir ${1} \
-                    --water-year ${2} \
-                    --variable albedo_observed_muZ
-
-if [ $? != 0 ]; then
-  echo "Error extracting variable from MODIS inputs"
-  exit 1
-fi
+extract_modis() {
+  # command from UofU snow-rs package
+  # https://github.com/UofU-Cryosphere/snow-rs
+  # Converts the MODIS matlab to GeoTiff files
+  variable_from_modis --source-dir ${1} \
+                      --water-year ${2} \
+                      --variable albedo_observed_muZ
+
+  if [ $? != 0 ]; then
+    echo "Error extracting variable from MODIS inputs"
+    exit 1
+  fi
+}
 
 export ERW_NC="_ERW.nc"
 export ERW_ONE_DAY_SUFFIX="_24.nc"
@@ -35,7 +37,7 @@ modis_erw() {
   ERW_TMP_NC=${1/\.tif/_tmp.nc}
 
   ERW_DOMAIN="\
--t_srs EPSG:32613 -tr 50 50 -dstnodata 65535
+-t_srs EPSG:32613 -tr 100 100 -dstnodata 65535
 -te 315900.0 4280850.0 348700.0 4322700.0
 "
   FILTER_MATH="A*(A<=numpy.max(B)) + numpy.max(B)*(A>numpy.max(B))" 
@@ -67,7 +69,7 @@ modis_erw() {
     return 1
   fi
 
-  date=$(date -d $(basename $1 | cut -d '_' -f2) +%Y-%m-%d)
+  date=$(date -d $(basename $1 | cut -d '_' -f1) +%Y-%m-%d)
   # Set time and variable name in NetCDF
   ${CDO_call} \
     -chname,Band1,albedo \

scripts/smeshr/MODIS-STC/MODIS_albedo_basin.sh

@@ -0,0 +1,154 @@
+#!/usr/bin/env bash
+
+# Script to extract albedo values for a model domain from MODIS.
+# This prepares a NetCDF file with per grid cell albedo values for every
+# hour of one day. It can be directly used to calculate net solar 
+# with incoming shortwave radiation.
+
+# Arguments:
+# ${1}: /path/to/MODIS/inputs
+# ${2}: water year
+
+set -e
+
+extract_modis() {
+  # command from UofU snow-rs package
+  # https://github.com/UofU-Cryosphere/snow-rs
+  # Converts the MODIS matlab to GeoTiff files
+  variable_from_modis --source-dir ${1} \
+                      --water-year ${2} \
+                      --variable albedo_observed_muZ
+
+  if [ $? != 0 ]; then
+    echo "Error extracting variable from MODIS inputs"
+    exit 1
+  fi
+}
+
+# Basin name and processing parameters
+#export BASIN='yampa'
+#BASIN_EXTENT='156726.4 4426192.0 365326.4 4530492.0' #yampa
+export BASIN='littleanimas_jmh'
+BASIN_EXTENT='251200 4170400 256400 4177500'
+RES=100
+EPSG='EPSG:32613'
+BASIN_DOMAIN=" -t_srs $EPSG -tr $RES $RES -dstnodata 65535 -te $BASIN_EXTENT"
+
+SLURM_NTASKS=4
+
+export BASIN_NC="_${BASIN}.nc"
+export BASIN_ONE_DAY_SUFFIX="_24.nc"
+
+export CDO_call="cdo -O -z zip_4 -f nc4 -s"
+export PARALLEL_call="parallel --tagstring {#} --tag --line-buffer --jobs ${SLURM_NTASKS}"
+
+# Extracts the BASIN domain from the Western US MODIS albedo GeoTiffs
+modis_basin() {
+  BASIN_TMP_CUBIC=${1/\.tif/_cubic_tmp.vrt}
+  BASIN_TMP_NN=${1/\.tif/_nn_tmp.vrt}
+  BASIN_TMP=${1/\.tif/_tmp.tif}
+  BASIN_TMP_NC=${1/\.tif/_tmp.nc}
+  BASIN_DOMAIN="${@:2}"
+  # add catch for existing 24.nc files (final output), remove later as needed
+  if [[ -f "${1%.*}_${BASIN}${BASIN_ONE_DAY_SUFFIX}" ]]; then
+    return
+  else
+    echo "${1%.*}_${BASIN}${BASIN_ONE_DAY_SUFFIX}" DNE
+    # echo "BASIN_DOMAIN: ${BASIN_DOMAIN}"
+    FILTER_MATH="A*(A<=numpy.max(B)) + numpy.max(B)*(A>numpy.max(B))" 
+    echo "gdalwarp -q -overwrite -multi \
+      -r cubic ${BASIN_DOMAIN} \
+      ${1} ${BASIN_TMP_CUBIC}"
+
+    # Use cubic resampling to smooth the image
+    gdalwarp -q -overwrite -multi \
+      -r cubic ${BASIN_DOMAIN} \
+      ${1} ${BASIN_TMP_CUBIC}
+
+    # Generate a nearest neighbor to get the max basin value as filter
+    gdalwarp -q -overwrite -multi \
+      ${BASIN_DOMAIN} \
+      ${1} ${BASIN_TMP_NN}
+
+    # Combine products, filter any values from cubic resampling higher
+    # than the max of nearest neighbor and set to that value
+    gdal_calc.py --overwrite --quiet --co COMPRESS=LZW \
+      -A ${BASIN_TMP_CUBIC} -B ${BASIN_TMP_NN} \
+      --calc="${FILTER_MATH}" \
+      --outfile=${BASIN_TMP}
+
+    # Add threshold enforcement of the values lower than the
+    # minimum of the nearest neighbor resampling and set to that min
+    gdal_calc.py --overwrite --quiet --co COMPRESS=LZW \
+      -A ${BASIN_TMP_CUBIC} -B ${BASIN_TMP_NN} \
+      --calc="A*(A>=numpy.min(B)) + numpy.min(B)*(A<numpy.min(B))" \
+      --outfile=${BASIN_TMP}
+
+    # Convert to NetCDF
+    gdalwarp -overwrite -multi -q \
+      -co FORMAT=NC4C -co WRITE_BOTTOMUP=NO \
+      ${BASIN_TMP} ${BASIN_TMP_NC}
+
+    if [ $? != 0 ]; then
+      printf "Error extracting domain from: \n   ${1}\n"
+      return 1
+    fi
+
+    date=$(date -d $(basename $1 | cut -d '_' -f3) +%Y-%m-%d)
+    # Set time and variable name in NetCDF
+    ${CDO_call} \
+      -chname,Band1,albedo \
+      -setdate,"${date}" \
+      -settunits,1hour \
+      ${BASIN_TMP_NC} ${1/\.tif/${BASIN_NC}}
+
+    if [ $? != 0 ]; then
+      echo "Error setting NetCDF variable name and time"
+      return 1
+    fi
+
+    printf "DONE: ${date}\n"
+
+    rm ${BASIN_TMP_CUBIC}
+    rm ${BASIN_TMP_NN}
+    rm ${BASIN_TMP}
+    rm ${BASIN_TMP_NC}
+  fi
+  }
+export -f modis_basin
+echo "Extracting MODIS albedo"
+${PARALLEL_call} modis_basin {} ${BASIN_DOMAIN} ::: ${1}/wy${2}/*.tif 
+
+# Creates a NetCDF using the extracted model domain values and
+# populates every hour of one day with those.
+albedo_day() {
+  HOUR_FILE=${1/\.nc$//}
+
+  if [[ ! -f "${1}" ]]; then
+    printf "Error: Source NetCDF:  \n  ${1}\n  does not exist\n"
+    return 1
+  fi
+
+  # Duplicate the single day value for every hour
+  for hour in {1..23}; do
+    ${CDO_call} -shifttime,${hour}hour \
+      ${1} ${HOUR_FILE}_${hour}.nc
+  done
+
+  # Make on file with values for every full hour and given file
+  ${CDO_call} mergetime ${1} ${HOUR_FILE}*.nc ${1/\.nc/${BASIN_ONE_DAY_SUFFIX}}
+
+  if [ $? != 0 ]; then
+    printf "Error: Could not merge hours for:\n  ${1}\n"
+    return 1
+  fi
+
+  printf "*\n"
+
+  rm ${HOUR_FILE}*.nc
+  rm ${1}
+}
+
+export -f albedo_day
+echo "Creating NetCDFs"
+${PARALLEL_call} albedo_day ::: ${1}/wy${2}/*${BASIN_NC}

scripts/smeshr/net_dswrf_for_topo_MODIS.sh → scripts/smeshr/MODIS-STC/net_dswrf_for_topo_MODIS.sh

@@ -47,7 +47,7 @@ function net_hrrr_for_month() {
 
   echo "  Merge HRRR month"
   MONTH_FILE="${DSWRF_OUT}/dswrf.${MONTH_SELECTOR}.nc"
-  ${CDO_COMMAND} mergetime ${HRRR_SELECT} ${DSWRF_IN}/*${LAST_DAY}* ${DSWRF_IN}/*${MONTH_SELECTOR}* ${MONTH_FILE}
+  ${CDO_COMMAND} ${HRRR_SELECT} [ -mergetime ${DSWRF_IN}/*${LAST_DAY}* ${DSWRF_IN}/*${MONTH_SELECTOR}* ] ${MONTH_FILE}
 
   if [[ $? != 0 ]]; then
     echo "  ** Error merging HRRR month **"

scripts/smeshr/README.md

@@ -5,15 +5,58 @@ Snow Melt Energy from Shortwave Radiation (SMESHR)
 Utilities to prepare using HRRR cloud cover (TCDC) and downward short-wave 
 radiation flux (DSWRF) in iSnobal as forcing inputs.
 
+## Required Tools
+Most scripts rely on these commands to be available/installed:
+* [Climate Data Operators (CDO)](https://code.mpimet.mpg.de/projects/cdo)
+* [GNU parallel](https://www.gnu.org/software/parallel/parallel_tutorial.html)
+
+Both are available via [conda forge](https://conda-forge.org/)
+
+All below scripts also require ths `snobedo` package of this repository to be 
+installed. See the [installation instructions](../../package/README.md) for 
+getting this step complete.
+
+### Conda
+The provided [environment.yaml](/environment.yaml) in this directory
+can be used to setup a execution environment. It also includes the above
+listed packages and the libraries needed for the `snobedo` package.
+
 ## Steps
+Below the order execution of scripts to create forcing input files for iSnobal
+in the expected format. Each script has it's functionality described in the
+header.
+
 ### For cloud cover
-Run `tcdc_for_topo.sh`
-Run `HRRR_linker.sh`
+1. `tcdc_for_topo.sh`
+1. `HRRR_linker.sh`
 
-### For short-wave radiation
+### For shortwave radiation
+There are two options to use the HRRR variable and produce a `net_solar.nc` 
+required to run iSnobal:
+* Combine with Time-Decay albedo (SMRF default)
+* Use the MODIS STC product
 
-Execution steps when using SMRF albedo:
+#### Common steps in both cases
+Extract the HRRR shortwave information for model domain
 1. `dswrf_for_topo.sh`
-2. Run SMRF and output variables `albedo_vis` and `albedo_ir`
-3. `net_dswrf_for_topo.sh`
-4. `HRRR_linker.sh`
+
+#### MODIS Albedo
+__REQUIRED SOFTWARE__:  
+The UofU [snow-rs](https://github.com/UofU-Cryosphere/snow-rs) package is 
+required to extract the relevant albedo variable from the MODIS STC product.
+The supplied conda environment of that package is not needed and already
+included with the provided `environment.yaml` of this directory.
+
+__IMPORTANT__:  
+The scripts are currently tailored to the East-River domain.
+Please adjust the domain extent as needed in `MODIS_albedo_ERW.sh` line 37
+
+Steps:  
+1. `MODIS_albedo_ERW.sh`
+1. `net_dswrf_for_topo_MODIS.sh`
+1. `HRRR_linker.sh`
+
+#### SMRF Albedo
+1. Run SMRF and output variables `albedo_vis` and `albedo_ir`
+1. `net_dswrf_for_topo_SMRF.sh`
+1. `HRRR_linker.sh`

scripts/smeshr/net_dswrf_for_topo_SMRF.sh → scripts/smeshr/SMRF/net_dswrf_for_topo_SMRF.sh

@@ -22,7 +22,7 @@ export WATER_START_MONTH=10
 export DSWRF_IN=${3}
 export SMRF_IN=${4}
 export DSWRF_OUT=${5}
-export ERW_DAY_MST="${DSWRF_OUT}/net_dswrf.MST"
+export BASIN_DAY_MST="${DSWRF_OUT}/net_dswrf.MST"
 
 export CDO_COMMAND='cdo -z zip4 -O'
 export HRRR_SELECT="-select,name=illumination_angle,DSWRF"
@@ -49,7 +49,8 @@ function net_hrrr_for_month() {
 
   echo "  Merge HRRR month"
   MONTH_FILE="${DSWRF_OUT}/dswrf.${MONTH_SELECTOR}.nc"
-  ${CDO_COMMAND} mergetime ${HRRR_SELECT} ${DSWRF_IN}.${LAST_DAY}* ${DSWRF_IN}.${MONTH_SELECTOR}* ${MONTH_FILE}
+  ${CDO_COMMAND} ${HRRR_SELECT} [ -mergetime ${DSWRF_IN}/*${LAST_DAY}* ${DSWRF_IN}/*${MONTH_SELECTOR}* ] ${MONTH_FILE}
+
 
   if [[ $? != 0 ]]; then
     echo "  ** Error merging HRRR month **"
@@ -80,7 +81,7 @@ function net_hrrr_for_month() {
   ${CDO_COMMAND} -aexpr,"${NET_MATH}" ${MERGE_FILE} ${MONTH_CALC_FILE}
 
   echo "  Split by day MST"
-  ${CDO_COMMAND} splitday ${MONTH_CALC_FILE} ${ERW_DAY_MST}.${MONTH_SELECTOR}
+  ${CDO_COMMAND} splitday ${MONTH_CALC_FILE} ${BASIN_DAY_MST}.${MONTH_SELECTOR}
 
   if [[ $? != 0 ]]; then
     echo "  ** Error processing ${MONTH_SELECTOR} **"

scripts/smeshr/dswrf_for_topo.sh

@@ -1,7 +1,7 @@
 #!/usr/bin/env bash
 # Extract DSWRF from given HRRR files for given topo extent
 #
-# NOTE: First argument needs to bin quotes to prevent shell expansion
+# NOTE: First argument needs to be in quotes to prevent shell expansion
 #
 # Arguments:
 #   ./dswrf_for_topo.sh <TOPO> <FOLDER_PATTERN> <HRRR_FILE_PATTERN> <OUTPUT_PATH_WITH_PREFIX>

scripts/smeshr/environment.yaml

@@ -0,0 +1,20 @@
+name: smeshr
+channels:
+- conda-forge
+dependencies:
+- cdo=2.2
+- cython
+- dask
+- dask-jobqueue
+- gdal=3.5
+- h5py
+- parallel
+- pip
+- python<3.11
+- skyfield
+- utm
+- wgrib2
+- netCDF4
+- numpy<2
+- pip:
+  - topocalc