diff --git a/R/gdal_cli.R b/R/gdal_cli.R
index d9cf91e66..3ed0be609 100644
--- a/R/gdal_cli.R
+++ b/R/gdal_cli.R
@@ -401,7 +401,10 @@ gdal_run <- function(cmd, args, close = FALSE, quiet = FALSE,
             cli::cli_alert_warning("Error reported during algorithm finalize.")
     }
 
-    return(alg)
+    if (quiet)
+        return(invisible(alg))
+    else
+        return(alg)
 }
 
 #' @name gdal_cli
@@ -726,7 +729,10 @@ gdal_global_reg_names <- function() {
     }
 
     if (alginfo$name == "pipeline") {
-        x <- alg$usageAsJSON() |> yyjsonr::read_json_str()
+        # workaround for non-standard Infinity value in GDAL JSON
+        x <- gsub("Infinity", '"Infinity"', alg$usageAsJSON(), fixed = TRUE) |>
+            yyjsonr::read_json_str()
+
         if (!is.null(x$pipeline_algorithms) &&
             is.data.frame(x$pipeline_algorithms) &&
             nrow(x$pipeline_algorithms) > 0) {
diff --git a/inst/extdata/storml_elev_pal.txt b/inst/extdata/storml_elev_pal.txt
new file mode 100644
index 000000000..a391e13c7
--- /dev/null
+++ b/inst/extdata/storml_elev_pal.txt
@@ -0,0 +1,6 @@
+0% 0 166 14
+20% 129 205 0
+40% 232 214 37
+60% 235 182 97
+80% 239 191 167
+100% 242 242 242
diff --git a/src/gdalalg.cpp b/src/gdalalg.cpp
index 18f5b331b..2a3cf1efa 100644
--- a/src/gdalalg.cpp
+++ b/src/gdalalg.cpp
@@ -64,7 +64,8 @@ void append_subalg_names_desc_(const GDALAlgorithmH alg,
         GDALAlgorithmH subalg = nullptr;
         subalg = GDALAlgorithmInstantiateSubAlgorithm(alg, subnames[i]);
         if (!subalg) {
-            cli_alert_danger_("failed to instantiate alg name: "s + subnames[i]);
+            cli_alert_danger_(
+                "failed to instantiate alg name: "s + subnames[i]);
             continue;
         }
 
@@ -682,18 +683,19 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
                                arg_name_in.get_cstring());
     if (!hArg) {
         if (!quiet)
-            cli_alert_danger_("failed to set algorithm argument with the value "
-                              "given for {.arg "s + arg_name.get_cstring() +
-                              "}");
+            cli_alert_danger_("failed to instantiate algorithm argument named "
+                              "{.str "s + arg_name_in.get_cstring() + "}");
         return false;
     }
     if (!GDALAlgorithmArgIsInput(hArg)) {
         if (!quiet)
-            cli_alert_danger_("{.str "s + arg_name.get_cstring() + "} is not "
-                              "an input argument of the algorithm");
+            cli_alert_danger_("{.str "s + arg_name_in.get_cstring() + "} is "
+                              "not an input argument of the algorithm");
         return false;
     }
 
+    std::string true_arg_name = GDALAlgorithmArgGetName(hArg);
+
     bool ret = false;
 
     switch (GDALAlgorithmArgGetType(hArg)) {
@@ -702,7 +704,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
             if (!Rcpp::is<Rcpp::LogicalVector>(arg_value)) {
                 cli_alert_danger_("logical value required for a BOOLEAN type "
                                   "algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
 
@@ -710,7 +712,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
             if (v.size() != 1 || Rcpp::LogicalVector::is_na(v[0])) {
                 cli_alert_danger_("single logical value (TRUE|FALSE) required "
                                   "for BOOLEAN type algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
 
@@ -724,7 +726,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
             if (!Rcpp::is<Rcpp::CharacterVector>(arg_value)) {
                 cli_alert_danger_("character value required for a STRING type "
                                   "algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
 
@@ -732,7 +734,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
             if (v.size() != 1 || Rcpp::CharacterVector::is_na(v[0])) {
                 cli_alert_danger_("single character string is required for "
                                   "a STRING algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
 
@@ -749,7 +751,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
 
                 cli_alert_danger_("integer or numeric value required for an "
                                   "INTEGER type algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
 
@@ -757,7 +759,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
             if (v.size() != 1) {
                 cli_alert_danger_("single numeric value is required for an "
                                   "INTEGER type algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
             if (Rcpp::IntegerVector::is_na(v[0])) {
@@ -777,7 +779,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
 
                 cli_alert_danger_("numeric value required for a REAL type "
                                   "algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
 
@@ -785,7 +787,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
             if (v.size() != 1) {
                 cli_alert_danger_("single numeric value is required for a "
                                   "REAL type algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
             if (Rcpp::NumericVector::is_na(v[0])) {
@@ -802,7 +804,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
             if (!Rcpp::is<Rcpp::CharacterVector>(arg_value)) {
                 cli_alert_danger_("character vector required for a STRINGLIST "
                                   "algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
 
@@ -814,7 +816,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
             if (Rcpp::is_true(Rcpp::any(Rcpp::is_na(v)))) {
                 cli_alert_danger_("{.arg arg_value} cannot contain {.val NA} "
                                   "for a STRING_LIST algorithm argument"s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
 
@@ -837,7 +839,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
 
                 cli_alert_danger_("integer or numeric vector required for an "
                                   "INTEGER_LIST algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
 
@@ -864,7 +866,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
 
                 cli_alert_danger_("numeric vector required for a REAL_LIST "
                                   "algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
                 break;
             }
 
@@ -885,20 +887,36 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
 
         case GAAT_DATASET:
         {
+            if (m_map_input_hDS.count(arg_name_in.get_cstring()) > 0 ||
+                m_map_input_hDS.count(true_arg_name) > 0) {
+
+                cli_alert_danger_("this argument was set in {.cls GDALAlg} "
+                                  "object instantiation, cannot set manually");
+                break;
+            }
+
+            if (EQUAL(GDALAlgorithmArgGetName(hArg), "like") &&
+                m_map_input_hDS.count("like") > 0) {
+
+                cli_alert_danger_(".{str like} arg was set in {.cls GDALAlg} "
+                                  "object instantiation, cannot set manually");
+                break;
+            }
+
             int ds_input_flags = GDALAlgorithmArgGetDatasetInputFlags(hArg);
 
             if (Rcpp::is<Rcpp::CharacterVector>(arg_value)) {
                 if (!(ds_input_flags & GADV_NAME)) {
                     cli_alert_danger_("argument does not accept a dataset "
                                       "name as input (object required): "s +
-                                      arg_name.get_cstring());
+                                      arg_name_in.get_cstring());
                     break;
                 }
                 Rcpp::CharacterVector v(arg_value);
                 if (v.size() != 1 || Rcpp::CharacterVector::is_na(v[0])) {
                     cli_alert_danger_("string input must be a length-1 "
                                       "character vector for DATASET algorithm "
-                                      "argument"s + arg_name.get_cstring());
+                                      "argument"s + arg_name_in.get_cstring());
                     break;
                 }
                 Rcpp::String val(enc_to_utf8_(v));
@@ -910,14 +928,15 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
                     cli_alert_danger_("argument does not accept a dataset "
                                       "object as input (DSN required, may be "
                                       "created by algorithm)"s +
-                                      arg_name.get_cstring());
+                                      arg_name_in.get_cstring());
                     break;
                 }
                 const Rcpp::String cls = arg_value.attr("class");
                 if (cls == "Rcpp_GDALRaster") {
                     const GDALRaster &ds = Rcpp::as<GDALRaster &>(arg_value);
-                    ret = GDALAlgorithmArgSetDataset(hArg,
-                                                     ds.getGDALDatasetH_());
+                    ret = GDALAlgorithmArgSetDataset(
+                        hArg, ds.getGDALDatasetH_());
+
                     break;
                 }
                 else if (cls == "Rcpp_GDALVector") {
@@ -925,32 +944,74 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
                     ret = GDALAlgorithmArgSetDataset(
                         hArg, ds.getGDALDatasetH_());
 
+                    // FIXME: "input" or "dataset" should be okay as of
+                    // GDAL 3.13 for vector algorithms, but probably
+                    // should not rely on hard coding the names
+                    if (ret) {
+                        if (EQUAL(arg_name_in.get_cstring(), "input") ||
+                            EQUAL(arg_name_in.get_cstring(), "dataset")) {
+
+                            m_in_vector_props.arg_name = true_arg_name;
+                            m_in_vector_props.driver_short_name =
+                                ds.getDriverShortName();
+                            m_in_vector_props.layer_name = ds.getName();
+                            m_in_vector_props.is_sql = ds.m_is_sql;
+                            m_in_vector_props.layer_sql = ds.m_layer_name;
+                            m_in_vector_props.sql_dialect = ds.m_dialect;
+                            m_in_vector_props.is_set = true;
+                        }
+                        else if (EQUAL(arg_name_in.get_cstring(), "like")) {
+                            m_like_vector_props.driver_short_name =
+                                ds.getDriverShortName();
+                            m_like_vector_props.layer_name = ds.getName();
+                            m_like_vector_props.is_sql = ds.m_is_sql;
+                            m_like_vector_props.layer_sql = ds.m_layer_name;
+                            m_like_vector_props.sql_dialect = ds.m_dialect;
+                            m_like_vector_props.is_set = true;
+                        }
+                    }
                     break;
                 }
                 else {
                     cli_alert_danger_("unhandled object of class {.cls "s +
                                       cls.get_cstring() + "} given for: " +
-                                      arg_name.get_cstring());
+                                      arg_name_in.get_cstring());
                     break;
                 }
             }
             else {
                 cli_alert_danger_("unsupported input type for a DATASET "
                                   "algorithm argument: "s +
-                                  arg_name.get_cstring());
+                                  arg_name_in.get_cstring());
             }
         }
         break;
 
         case GAAT_DATASET_LIST:
         {
+            if (m_map_input_hDS.count(arg_name_in.get_cstring()) > 0 ||
+                m_map_input_hDS.count(true_arg_name) > 0) {
+
+                cli_alert_danger_("this argument was set in {.cls GDALAlg} "
+                                  "object instantiation, cannot set manually");
+                break;
+            }
+
+            if (EQUAL(GDALAlgorithmArgGetName(hArg), "like") &&
+                m_map_input_hDS.count("like") > 0) {
+
+                cli_alert_danger_(".{str like} arg was set in {.cls GDALAlg} "
+                                  "object instantiation, cannot set manually");
+                break;
+            }
+
             int ds_input_flags = GDALAlgorithmArgGetDatasetInputFlags(hArg);
 
             if (Rcpp::is<Rcpp::CharacterVector>(arg_value)) {
                 if (!(ds_input_flags & GADV_NAME)) {
                     cli_alert_danger_("argument does not accept dataset "
                                       "names as input: "s +
-                                      arg_name.get_cstring());
+                                      arg_name_in.get_cstring());
                     break;
                 }
                 Rcpp::CharacterVector v(arg_value);
@@ -980,7 +1041,7 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
                 if (!(ds_input_flags & GADV_OBJECT)) {
                     cli_alert_danger_("argument does not accept dataset "
                                       "objects as input: "s +
-                                      arg_name.get_cstring());
+                                      arg_name_in.get_cstring());
                     break;
                 }
 
@@ -991,6 +1052,8 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
                     ds_list = Rcpp::List(arg_value);
 
                 std::vector<GDALDatasetH> pahDS;
+                bool in_vector_props_were_set = false;
+                int nInputVectorDatasets = 0;
                 for (R_xlen_t i = 0; i < ds_list.size(); ++i) {
                     if (is_gdalraster_obj_(ds_list[i])) {
                         const Rcpp::RObject &x(ds_list[i]);
@@ -1002,25 +1065,50 @@ bool GDALAlg::setArg(const Rcpp::String &arg_name,
                         else if (cls == "Rcpp_GDALVector") {
                             const GDALVector &ds = Rcpp::as<GDALVector &>(x);
                             pahDS.push_back(ds.getGDALDatasetH_());
+                            ++nInputVectorDatasets;
+
+                            // FIXME: "input" or "dataset" should be okay as of
+                            // GDAL 3.13 for vector algorithms, but probably
+                            // should not rely on hard coding the names
+                            if ((EQUAL(arg_name_in.get_cstring(), "input") ||
+                                 EQUAL(arg_name_in.get_cstring(), "dataset")) &&
+                                !in_vector_props_were_set) {
+
+                                m_in_vector_props.arg_name = true_arg_name;
+                                m_in_vector_props.driver_short_name =
+                                    ds.getDriverShortName();
+                                m_in_vector_props.layer_name = ds.getName();
+                                m_in_vector_props.is_sql = ds.m_is_sql;
+                                m_in_vector_props.layer_sql = ds.m_layer_name;
+                                m_in_vector_props.sql_dialect = ds.m_dialect;
+                                m_in_vector_props.is_set = true;
+                                in_vector_props_were_set = true;
+                            }
                         }
                         else {
                             cli_alert_danger_("unhandled object of class "
                                               "{.cls "s + cls.get_cstring() +
                                               "} given for: " +
-                                              arg_name.get_cstring());
+                                              arg_name_in.get_cstring());
                             break;
                         }
                     }
                     else {
                         cli_alert_danger_("an element contains invalid input, "
                                           "in the list given for: "s +
-                                          arg_name.get_cstring());;
+                                          arg_name_in.get_cstring());;
                         break;
                     }
                 }
                 if (!pahDS.empty()) {
                     ret = GDALAlgorithmArgSetDatasets(hArg, pahDS.size(),
                                                       pahDS.data());
+
+                    if (nInputVectorDatasets > 1) {
+                        // since currently, auto setting args from GDALVector
+                        // input is not supported for multiple input datasets
+                        m_in_vector_props = {};
+                    }
                     break;
                 }
             }
@@ -1110,7 +1198,9 @@ bool GDALAlg::parseCommandLineArgs() {
         // directly set certain arguments from input GDALVector object(s)
 
         if (m_in_vector_props.is_set) {
-            GDALAlgorithmArgH hArg = GDALAlgorithmGetArg(m_hAlg, "input");
+            GDALAlgorithmArgH hArg = GDALAlgorithmGetArg(
+                m_hAlg, m_in_vector_props.arg_name.c_str());
+
             if (GDALAlgorithmArgGetType(hArg) == GAAT_DATASET_LIST) {
                 if (m_num_input_datasets > 1) {
                     GDALAlgorithmArgRelease(hArg);
@@ -1121,7 +1211,7 @@ bool GDALAlg::parseCommandLineArgs() {
             }
             else if (GDALAlgorithmArgGetType(hArg) != GAAT_DATASET) {
                 GDALAlgorithmArgRelease(hArg);
-                Rcpp::stop("setting args from GDALVector incompatible with "
+                Rcpp::stop("setting args from `GDALVector` incompatible with "
                            "the algorithm argument type of \"input\" arg");
             }
             GDALAlgorithmArgRelease(hArg);
@@ -1244,6 +1334,8 @@ bool GDALAlg::parseCommandLineArgs() {
         m_haveParsedCmdLineArgs = true;
         if (m_hActualAlg == nullptr) {
             m_hActualAlg = GDALAlgorithmGetActualAlgorithm(m_hAlg);
+            if (!m_hActualAlg)
+                res = false;
         }
     }
 
@@ -1522,7 +1614,7 @@ void GDALAlg::show() const {
 // ****************************************************************************
 
 Rcpp::CharacterVector GDALAlg::parseListArgs_(
-                                const Rcpp::List &list_args) {
+    const Rcpp::List &list_args) {
 
     // convert arguments in a named list to character vector form and return it
     // arguments in list form must use arg long names
@@ -1558,7 +1650,7 @@ Rcpp::CharacterVector GDALAlg::parseListArgs_(
             cli_alert_danger_("argument: "s + arg_names[i]);
             Rcpp::stop("arguments in list format must use \"long\" names");
         }
-        std::string nm_no_lead_dashes(nm.get_cstring());
+        std::string nm_no_lead_dashes(nm);
         nm.push_front("--");
 
         // boolean
@@ -1608,10 +1700,15 @@ Rcpp::CharacterVector GDALAlg::parseListArgs_(
                         const GDALVector &ds =
                             Rcpp::as<GDALVector &>(list_tmp[j]);
 
+                        // FIXME: "input" or "dataset" should be okay as of
+                        // GDAL 3.13 for vector algorithms, but probably
+                        // should not rely on hard coding the names
                         ds_list.push_back(ds.getGDALDatasetH_());
-                        if (EQUAL(nm.get_cstring(), "--input") &&
+                        if ((EQUAL(nm.get_cstring(), "--input") ||
+                             EQUAL(nm.get_cstring(), "--dataset")) &&
                             !m_in_vector_props.is_set) {
 
+                            m_in_vector_props.arg_name = nm_no_lead_dashes;
                             m_in_vector_props.driver_short_name =
                                 ds.getDriverShortName();
                             m_in_vector_props.layer_name = ds.getName();
@@ -1654,7 +1751,14 @@ Rcpp::CharacterVector GDALAlg::parseListArgs_(
                 std::vector<GDALDatasetH> ds_list = {};
                 ds_list.push_back(ds.getGDALDatasetH_());
                 m_map_input_hDS[nm_no_lead_dashes] = ds_list;
-                if (EQUAL(nm.get_cstring(), "--input")) {
+
+                // FIXME: "input" or "dataset" should be okay as of
+                // GDAL 3.13 for vector algorithms, but probably
+                // should not rely on hard coding the names
+                if (EQUAL(nm.get_cstring(), "--input") ||
+                    EQUAL(nm.get_cstring(), "--dataset")) {
+
+                    m_in_vector_props.arg_name = nm_no_lead_dashes;
                     m_in_vector_props.driver_short_name =
                         ds.getDriverShortName();
                     m_in_vector_props.layer_name = ds.getName();
@@ -1680,10 +1784,19 @@ Rcpp::CharacterVector GDALAlg::parseListArgs_(
                           nm_no_lead_dashes.c_str());
     }
 
-    if (m_map_input_hDS.count("input") > 0) {
-        m_input_is_object = true;
-        std::vector<GDALDatasetH> ds_list = m_map_input_hDS["input"];
-        m_num_input_datasets = ds_list.size();
+    if (m_map_input_hDS.count("input") > 0 ||
+        m_map_input_hDS.count("dataset") > 0) {
+
+        if (m_map_input_hDS.count("input") > 0) {
+            std::vector<GDALDatasetH> ds_list = m_map_input_hDS["input"];
+            m_num_input_datasets += ds_list.size();
+        }
+
+        if (m_map_input_hDS.count("dataset") > 0) {
+            std::vector<GDALDatasetH> ds_list = m_map_input_hDS["dataset"];
+            m_num_input_datasets += ds_list.size();
+        }
+
         if (m_num_input_datasets > 1 && m_in_vector_props.is_set) {
             // since currently, auto setting args from GDALVector input is
             // not supported for multiple input datasets
@@ -1914,7 +2027,7 @@ SEXP GDALAlg::getArgValue_(const GDALAlgorithmArgH &hArg) const {
             hArgDSValue = GDALAlgorithmArgGetAsDatasetValue(hArg);
             if (!hArgDSValue) {
                 if (!quiet)
-                    cli_alert_warning_( "output dataset value is NULL");
+                    cli_alert_warning_("output dataset value is NULL");
                 return R_NilValue;
             }
 
diff --git a/src/gdalalg.h b/src/gdalalg.h
index 56ca2058e..77680ecf6 100644
--- a/src/gdalalg.h
+++ b/src/gdalalg.h
@@ -26,6 +26,7 @@ Rcpp::DataFrame gdal_commands(const std::string &contains, bool recurse,
 Rcpp::CharacterVector gdal_global_reg_names();
 
 struct VectorObjectProperties {
+    std::string arg_name {};
     bool is_set {false};
     std::string driver_short_name {};
     std::string layer_name {};
@@ -86,7 +87,6 @@ class GDALAlg {
     GDALAlgorithmH m_hAlg {nullptr};
     GDALAlgorithmH m_hActualAlg {nullptr};
 #endif
-    bool m_input_is_object {false};
     std::map<std::string, std::vector<GDALDatasetH>> m_map_input_hDS {};
     size_t m_num_input_datasets {0};
     VectorObjectProperties m_in_vector_props;
diff --git a/tests/testthat/test-GDALAlg-class.R b/tests/testthat/test-GDALAlg-class.R
index d30c638e2..f1236a742 100644
--- a/tests/testthat/test-GDALAlg-class.R
+++ b/tests/testthat/test-GDALAlg-class.R
@@ -750,3 +750,31 @@ test_that("setArg works", {
     out_ds$close()
     ds$close()
 })
+
+test_that("raster compare works", {
+    # test an algorithm with an atypical argument name for an input dataset
+    skip_if(gdal_version_num() < gdal_compute_version(3, 12, 0))
+
+    f_ref <- system.file("extdata/storml_elev.tif", package="gdalraster")
+    f2 <- system.file("extdata/storml_elev_orig.tif", package="gdalraster")
+    ds_ref <- new(GDALRaster, f_ref)
+    on.exit(ds_ref$close(), add = TRUE)
+    ds2 <- new(GDALRaster, f2)
+    on.exit(ds2$close(), add = TRUE)
+    args <- list(reference = ds_ref, input = ds2)
+    alg <- new(GDALAlg, "raster compare", args)
+    expect_true(alg$run())
+    expect_true(
+        grepl("pixels differing", alg$outputs()$output_string, fixed = TRUE))
+
+    alg$release()
+
+    ds_mem <- createCopy("MEM", "", f_ref, return_obj = TRUE)
+    on.exit(ds_mem$close(), add = TRUE)
+    args <- list(reference = ds_ref, input = ds_mem, skip_binary = TRUE)
+    alg <- new(GDALAlg, "raster compare", args)
+    expect_true(alg$run())
+    expect_equal(alg$outputs()$output_string, "")
+
+    alg$release()
+})
diff --git a/tests/testthat/test-gdal_cli.R b/tests/testthat/test-gdal_cli.R
index 3939b9737..bc4259b9a 100644
--- a/tests/testthat/test-gdal_cli.R
+++ b/tests/testthat/test-gdal_cli.R
@@ -166,8 +166,42 @@ test_that("gdal_usage works", {
     expect_output(gdal_usage(cmd), "Options:")
     expect_output(gdal_usage(cmd), "Advanced options:")
     expect_output(gdal_usage(cmd), "For more details:")
+
+    expect_no_error(gdal_usage("pipeline"))
 })
 
 test_that("gdal_global_reg_names returns a character vector", {
     expect_vector(gdal_global_reg_names(), character())
 })
+
+test_that("raster pipeline works", {
+    ## test raster pipeline algorithms
+    skip_if(gdal_version_num() < gdal_compute_version(3, 12, 0))
+
+    ## with a nested input pipeline
+    f <- system.file("extdata/storml_elev.tif", package="gdalraster")
+    f_elev <- tempfile(fileext = ".tif")
+    file.copy(f, f_elev)
+    on.exit(deleteDataset(f_elev), add = TRUE)
+    f_pal <- system.file("extdata/storml_elev_pal.txt", package="gdalraster")
+    f_out <- file.path(tempdir(), "storml_col_relief.tif")
+    on.exit(deleteDataset(f_out), add = TRUE)
+
+    args <- paste(
+        "read --input", f_elev,
+        "! color-map --color-map", f_pal,
+        "! blend --overlay [ read --input", f_elev, "! hillshade -z 1.5 ]",
+            "--operator=hsv-value",
+        "! write --output", f_out, "--overwrite")
+
+    expect_no_error(alg <- gdal_run("raster pipeline", args))
+    expect_true(is.list(alg$outputs()))
+    ds <- alg$outputs()$output
+    expect_true(is(ds, "Rcpp_GDALRaster"))
+    expect_equal(ds$res(), c(30, 30))
+    expect_equal(ds$dim(), c(143, 107, 3))
+
+    ds$close()
+    expect_true(alg$close())
+    alg$release()
+})
\ No newline at end of file
diff --git a/vignettes/articles/img/storml_col_shaded_relief.png b/vignettes/articles/img/storml_col_shaded_relief.png
new file mode 100644
index 000000000..fcd2b46c3
Binary files /dev/null and b/vignettes/articles/img/storml_col_shaded_relief.png differ
diff --git a/vignettes/articles/use-gdal-cli-from-r.Rmd b/vignettes/articles/use-gdal-cli-from-r.Rmd
index 07084cfb9..a0f348b96 100644
--- a/vignettes/articles/use-gdal-cli-from-r.Rmd
+++ b/vignettes/articles/use-gdal-cli-from-r.Rmd
@@ -34,40 +34,44 @@ The function `gdal_commands()` prints a list of available commands to the consol
 
 ``` r
 library(gdalraster)
-#> GDAL 3.11.3 (released 2025-07-12), GEOS 3.12.2, PROJ 9.4.1
+#> GDAL 3.12.1 (released 2025-12-12), GEOS 3.12.2, PROJ 9.4.1
 
 ## top-level commands
 gdal_commands(recurse = FALSE)
-#> convert:
-#> Convert a dataset (shortcut for 'gdal raster convert' or 'gdal vector convert').
-#> https://gdal.org/programs/gdal_convert.html
+#> → "convert"
+#> ℹ Convert a dataset (shortcut for 'gdal raster convert' or 'gdal vector convert').
+#> ℹ <https://gdal.org/programs/gdal_convert.html>
 #> 
-#> driver:
-#> Command for driver specific operations.
+#> → "dataset"
+#> ℹ Commands to manage datasets.
+#> ℹ <https://gdal.org/programs/gdal_dataset.html>
 #> 
-#> info:
-#> Return information on a dataset (shortcut for 'gdal raster info' or 'gdal vector info').
-#> https://gdal.org/programs/gdal_info.html
+#> → "driver"
+#> ℹ Command for driver specific operations.
 #> 
-#> mdim:
-#> Multidimensional commands.
-#> https://gdal.org/programs/gdal_mdim.html
+#> → "info"
+#> ℹ Return information on a dataset (shortcut for 'gdal raster info' or 'gdal vector info').
+#> ℹ <https://gdal.org/programs/gdal_info.html>
 #> 
-#> pipeline:
-#> Execute a pipeline (shortcut for 'gdal raster pipeline' or 'gdal vector pipeline').
-#> https://gdal.org/programs/gdal_pipeline.html
+#> → "mdim"
+#> ℹ Multidimensional commands.
+#> ℹ <https://gdal.org/programs/gdal_mdim.html>
 #> 
-#> raster:
-#> Raster commands.
-#> https://gdal.org/programs/gdal_raster.html
+#> → "pipeline"
+#> ℹ Process a dataset applying several steps.
+#> ℹ <https://gdal.org/programs/gdal_pipeline.html>
 #> 
-#> vector:
-#> Vector commands.
-#> https://gdal.org/programs/gdal_vector.html
+#> → "raster"
+#> ℹ Raster commands.
+#> ℹ <https://gdal.org/programs/gdal_raster.html>
 #> 
-#> vsi:
-#> GDAL Virtual System Interface (VSI) commands.
-#> https://gdal.org/programs/gdal_vsi.html
+#> → "vector"
+#> ℹ Vector commands.
+#> ℹ <https://gdal.org/programs/gdal_vector.html>
+#> 
+#> → "vsi"
+#> ℹ GDAL Virtual System Interface (VSI) commands.
+#> ℹ <https://gdal.org/programs/gdal_vsi.html>
 ```
 
 A character string can also be given to filter for commands containing specific text.
@@ -75,157 +79,201 @@ A character string can also be given to filter for commands containing specific
 ``` r
 ## list commands relevant to raster data
 gdal_commands("raster")
-#> raster:
-#> Raster commands.
-#> https://gdal.org/programs/gdal_raster.html
+#> → "raster"
+#> ℹ Raster commands.
+#> ℹ <https://gdal.org/programs/gdal_raster.html>
+#> 
+#> → "raster as-features"
+#> ℹ Create features from pixels of a raster dataset
+#> ℹ <https://gdal.org/programs/gdal_raster_as_features.html>
+#> 
+#> → "raster aspect"
+#> ℹ Generate an aspect map
+#> ℹ <https://gdal.org/programs/gdal_raster_aspect.html>
+#> 
+#> → "raster blend"
+#> ℹ Blend/compose two raster datasets
+#> ℹ <https://gdal.org/programs/gdal_raster_blend.html>
+#> 
+#> → "raster calc"
+#> ℹ Perform raster algebra
+#> ℹ <https://gdal.org/programs/gdal_raster_calc.html>
+#> 
+#> → "raster clean-collar"
+#> ℹ Clean the collar of a raster dataset, removing noise.
+#> ℹ <https://gdal.org/programs/gdal_raster_clean_collar.html>
+#> 
+#> → "raster clip"
+#> ℹ Clip a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_clip.html>
+#> 
+#> → "raster color-map"
+#> ℹ Generate a RGB or RGBA dataset from a single band, using a color map
+#> ℹ <https://gdal.org/programs/gdal_raster_color_map.html>
+#> 
+#> → "raster compare"
+#> ℹ Compare two raster datasets.
+#> ℹ <https://gdal.org/programs/gdal_raster_compare.html>
+#> 
+#> → "raster contour"
+#> ℹ Creates a vector contour from a raster elevation model (DEM).
+#> ℹ <https://gdal.org/programs/gdal_raster_contour.html>
 #> 
-#> raster aspect:
-#> Generate an aspect map
-#> https://gdal.org/programs/gdal_raster_aspect.html
+#> → "raster convert"
+#> ℹ Convert a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_convert.html>
 #> 
-#> raster calc:
-#> Perform raster algebra
-#> https://gdal.org/programs/gdal_raster_calc.html
+#> → "raster create"
+#> ℹ Create a new raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_create.html>
 #> 
-#> raster clean-collar:
-#> Clean the collar of a raster dataset, removing noise.
-#> https://gdal.org/programs/gdal_raster_clean_collar.html
+#> → "raster edit"
+#> ℹ Edit a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_edit.html>
 #> 
-#> raster clip:
-#> Clip a raster dataset.
-#> https://gdal.org/programs/gdal_raster_clip.html
+#> → "raster fill-nodata"
+#> ℹ Fill nodata raster regions by interpolation from edges.
+#> ℹ <https://gdal.org/programs/gdal_raster_fill_nodata.html>
 #> 
-#> raster color-map:
-#> Generate a RGB or RGBA dataset from a single band, using a color map
-#> https://gdal.org/programs/gdal_raster_color_map.html
+#> → "raster footprint"
+#> ℹ Compute the footprint of a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_footprint.html>
 #> 
-#> raster contour:
-#> Creates a vector contour from a raster elevation model (DEM).
-#> https://gdal.org/programs/gdal_raster_contour.html
+#> → "raster hillshade"
+#> ℹ Generate a shaded relief map
+#> ℹ <https://gdal.org/programs/gdal_raster_hillshade.html>
 #> 
-#> raster convert:
-#> Convert a raster dataset.
-#> https://gdal.org/programs/gdal_raster_convert.html
+#> → "raster index"
+#> ℹ Create a vector index of raster datasets.
+#> ℹ <https://gdal.org/programs/gdal_raster_index.html>
 #> 
-#> raster create:
-#> Create a new raster dataset.
-#> https://gdal.org/programs/gdal_raster_create.html
+#> → "raster info"
+#> ℹ Return information on a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_info.html>
 #> 
-#> raster edit:
-#> Edit a raster dataset.
-#> https://gdal.org/programs/gdal_raster_edit.html
+#> → "raster mosaic"
+#> ℹ Build a mosaic, either virtual (VRT) or materialized.
+#> ℹ <https://gdal.org/programs/gdal_raster_mosaic.html>
 #> 
-#> raster fill-nodata:
-#> Fill nodata raster regions by interpolation from edges.
-#> https://gdal.org/programs/gdal_raster_fill_nodata.html
+#> → "raster neighbors"
+#> ℹ Compute the value of each pixel from its neighbors (focal statistics)
+#> ℹ <https://gdal.org/programs/gdal_raster_neighbors.html>
 #> 
-#> raster footprint:
-#> Compute the footprint of a raster dataset.
-#> https://gdal.org/programs/gdal_raster_footprint.html
+#> → "raster nodata-to-alpha"
+#> ℹ Replace nodata value(s) with an alpha band.
+#> ℹ <https://gdal.org/programs/gdal_raster_nodata_to_alpha.html>
 #> 
-#> raster hillshade:
-#> Generate a shaded relief map
-#> https://gdal.org/programs/gdal_raster_hillshade.html
+#> → "raster overview"
+#> ℹ Manage overviews of a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_overview.html>
 #> 
-#> raster index:
-#> Create a vector index of raster datasets.
-#> https://gdal.org/programs/gdal_raster_index.html
+#> → "raster overview add"
+#> ℹ Adding overviews.
+#> ℹ <https://gdal.org/programs/gdal_raster_overview_add.html>
 #> 
-#> raster info:
-#> Return information on a raster dataset.
-#> https://gdal.org/programs/gdal_raster_info.html
+#> → "raster overview delete"
+#> ℹ Deleting overviews.
+#> ℹ <https://gdal.org/programs/gdal_raster_overview_delete.html>
 #> 
-#> raster mosaic:
-#> Build a mosaic, either virtual (VRT) or materialized.
-#> https://gdal.org/programs/gdal_raster_mosaic.html
+#> → "raster overview refresh"
+#> ℹ Refresh overviews.
+#> ℹ <https://gdal.org/programs/gdal_raster_overview_refresh.html>
 #> 
-#> raster overview:
-#> Manage overviews of a raster dataset.
-#> https://gdal.org/programs/gdal_raster_overview.html
+#> → "raster pansharpen"
+#> ℹ Perform a pansharpen operation.
+#> ℹ <https://gdal.org/programs/gdal_raster_pansharpen.html>
 #> 
-#> raster overview add:
-#> Adding overviews.
-#> https://gdal.org/programs/gdal_raster_overview_add.html
+#> → "raster pipeline"
+#> ℹ Process a raster dataset applying several steps.
+#> ℹ <https://gdal.org/programs/gdal_raster_pipeline.html>
 #> 
-#> raster overview delete:
-#> Deleting overviews.
-#> https://gdal.org/programs/gdal_raster_overview_delete.html
+#> → "raster pixel-info"
+#> ℹ Return information on a pixel of a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_pixel_info.html>
 #> 
-#> raster pipeline:
-#> Process a raster dataset.
-#> https://gdal.org/programs/gdal_raster_pipeline.html
+#> → "raster polygonize"
+#> ℹ Create a polygon feature dataset from a raster band.
+#> ℹ <https://gdal.org/programs/gdal_raster_polygonize.html>
 #> 
-#> raster pixel-info:
-#> Return information on a pixel of a raster dataset.
-#> https://gdal.org/programs/gdal_raster_pixel_info.html
+#> → "raster proximity"
+#> ℹ Produces a raster proximity map.
+#> ℹ <https://gdal.org/programs/gdal_raster_proximity.html>
 #> 
-#> raster polygonize:
-#> Create a polygon feature dataset from a raster band.
-#> https://gdal.org/programs/gdal_raster_polygonize.html
+#> → "raster reclassify"
+#> ℹ Reclassify values in a raster dataset
+#> ℹ <https://gdal.org/programs/gdal_raster_reclassify.html>
 #> 
-#> raster reclassify:
-#> Reclassify values in a raster dataset
-#> https://gdal.org/programs/gdal_raster_reclassify.html
+#> → "raster reproject"
+#> ℹ Reproject a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_reproject.html>
 #> 
-#> raster reproject:
-#> Reproject a raster dataset.
-#> https://gdal.org/programs/gdal_raster_reproject.html
+#> → "raster resize"
+#> ℹ Resize a raster dataset without changing the georeferenced extents.
+#> ℹ <https://gdal.org/programs/gdal_raster_resize.html>
 #> 
-#> raster resize:
-#> Resize a raster dataset without changing the georeferenced extents.
-#> https://gdal.org/programs/gdal_raster_resize.html
+#> → "raster rgb-to-palette"
+#> ℹ Convert a RGB image into a pseudo-color / paletted image.
+#> ℹ <https://gdal.org/programs/gdal_raster_rgb_to_palette.html>
 #> 
-#> raster roughness:
-#> Generate a roughness map
-#> https://gdal.org/programs/gdal_raster_roughness.html
+#> → "raster roughness"
+#> ℹ Generate a roughness map
+#> ℹ <https://gdal.org/programs/gdal_raster_roughness.html>
 #> 
-#> raster scale:
-#> Scale the values of the bands of a raster dataset.
-#> https://gdal.org/programs/gdal_raster_scale.html
+#> → "raster scale"
+#> ℹ Scale the values of the bands of a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_scale.html>
 #> 
-#> raster select:
-#> Select a subset of bands from a raster dataset.
-#> https://gdal.org/programs/gdal_raster_select.html
+#> → "raster select"
+#> ℹ Select a subset of bands from a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_select.html>
 #> 
-#> raster set-type:
-#> Modify the data type of bands of a raster dataset.
-#> https://gdal.org/programs/gdal_raster_set_type.html
+#> → "raster set-type"
+#> ℹ Modify the data type of bands of a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_set_type.html>
 #> 
-#> raster sieve:
-#> Remove small polygons from a raster dataset.
-#> https://gdal.org/programs/gdal_raster_sieve.html
+#> → "raster sieve"
+#> ℹ Remove small polygons from a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_sieve.html>
 #> 
-#> raster slope:
-#> Generate a slope map
-#> https://gdal.org/programs/gdal_raster_slope.html
+#> → "raster slope"
+#> ℹ Generate a slope map
+#> ℹ <https://gdal.org/programs/gdal_raster_slope.html>
 #> 
-#> raster stack:
-#> Combine together input bands into a multi-band output, either virtual (VRT) or materialized.
-#> https://gdal.org/programs/gdal_raster_stack.html
+#> → "raster stack"
+#> ℹ Combine together input bands into a multi-band output, either virtual (VRT) or materialized.
+#> ℹ <https://gdal.org/programs/gdal_raster_stack.html>
 #> 
-#> raster tile:
-#> Generate tiles in separate files from a raster dataset.
-#> https://gdal.org/programs/gdal_raster_tile.html
+#> → "raster tile"
+#> ℹ Generate tiles in separate files from a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_tile.html>
 #> 
-#> raster tpi:
-#> Generate a Topographic Position Index (TPI) map
-#> https://gdal.org/programs/gdal_raster_tpi.html
+#> → "raster tpi"
+#> ℹ Generate a Topographic Position Index (TPI) map
+#> ℹ <https://gdal.org/programs/gdal_raster_tpi.html>
 #> 
-#> raster tri:
-#> Generate a Terrain Ruggedness Index (TRI) map
-#> https://gdal.org/programs/gdal_raster_tri.html
+#> → "raster tri"
+#> ℹ Generate a Terrain Ruggedness Index (TRI) map
+#> ℹ <https://gdal.org/programs/gdal_raster_tri.html>
 #> 
-#> raster unscale:
-#> Convert scaled values of a raster dataset into unscaled values.
-#> https://gdal.org/programs/gdal_raster_unscale.html
+#> → "raster unscale"
+#> ℹ Convert scaled values of a raster dataset into unscaled values.
+#> ℹ <https://gdal.org/programs/gdal_raster_unscale.html>
 #> 
-#> raster viewshed:
-#> Compute the viewshed of a raster dataset.
-#> https://gdal.org/programs/gdal_raster_viewshed.html
+#> → "raster update"
+#> ℹ Update the destination raster with the content of the input one.
+#> ℹ <https://gdal.org/programs/gdal_raster_update.html>
 #> 
-#> vector rasterize:
-#> Burns vector geometries into a raster.
-#> https://gdal.org/programs/gdal_vector_rasterize.html
+#> → "raster viewshed"
+#> ℹ Compute the viewshed of a raster dataset.
+#> ℹ <https://gdal.org/programs/gdal_raster_viewshed.html>
+#> 
+#> → "raster zonal-stats"
+#> ℹ Calculate raster zonal statistics
+#> ℹ <https://gdal.org/programs/gdal_raster_zonal_stats.html>
+#> 
+#> → "vector rasterize"
+#> ℹ Burns vector geometries into a raster.
+#> ℹ <https://gdal.org/programs/gdal_vector_rasterize.html>
 ```
 
 The function `gdal_usage()` prints a help message to the console for a specific command.
@@ -414,21 +462,38 @@ unlink(f_gpkg)
 
 ``` r
 gdal_usage("vector")
+#> 
 #> Usage: vector <SUBCOMMAND> [OPTIONS]
 #> where <SUBCOMMAND> is one of:
-#>   - clip      : Clip a vector dataset.
-#>   - concat    : Concatenate vector datasets.
-#>   - convert   : Convert a vector dataset.
-#>   - edit      : Edit metadata of a vector dataset.
-#>   - filter    : Filter a vector dataset.
-#>   - geom      : Geometry operations on a vector dataset.
-#>   - grid      : Create a regular grid from scattered points.
-#>   - info      : Return information on a vector dataset.
-#>   - pipeline  : Process a vector dataset.
-#>   - rasterize : Burns vector geometries into a raster.
-#>   - reproject : Reproject a vector dataset.
-#>   - select    : Select a subset of fields from a vector dataset.
-#>   - sql       : Apply SQL statement(s) to a dataset.
+#>   - buffer              : Compute a buffer around geometries of a vector dataset.
+#>   - check-coverage      : Check a polygon coverage for validity
+#>   - check-geometry      : Check a dataset for invalid geometries
+#>   - clean-coverage      : Alter polygon boundaries to make shared edges identical, removing gaps and overlaps
+#>   - clip                : Clip a vector dataset.
+#>   - concat              : Concatenate vector datasets.
+#>   - convert             : Convert a vector dataset.
+#>   - edit                : Edit metadata of a vector dataset.
+#>   - explode-collections : Explode geometries of type collection of a vector dataset.
+#>   - filter              : Filter a vector dataset.
+#>   - geom                : Geometry operations on a vector dataset.
+#>   - grid                : Create a regular grid from scattered points.
+#>   - index               : Create a vector index of vector datasets.
+#>   - info                : Return information on a vector dataset.
+#>   - layer-algebra       : Perform algebraic operation between 2 layers.
+#>   - make-point          : Create point geometries from attribute fields
+#>   - make-valid          : Fix validity of geometries of a vector dataset.
+#>   - partition           : Partition a vector dataset into multiple files.
+#>   - pipeline            : Process a vector dataset applying several steps.
+#>   - rasterize           : Burns vector geometries into a raster.
+#>   - reproject           : Reproject a vector dataset.
+#>   - segmentize          : Segmentize geometries of a vector dataset.
+#>   - select              : Select a subset of fields from a vector dataset.
+#>   - set-field-type      : Modify the type of a field of a vector dataset.
+#>   - set-geom-type       : Modify the geometry type of a vector dataset.
+#>   - simplify            : Simplify geometries of a vector dataset.
+#>   - simplify-coverage   : Simplify shared boundaries of a polygonal vector dataset.
+#>   - sql                 : Apply SQL statement(s) to a dataset.
+#>   - swap-xy             : Swap X and Y coordinates of geometries of a vector dataset.
 #> 
 #> Options:
 #>   --drivers
@@ -436,8 +501,7 @@ gdal_usage("vector")
 #>   --output-string <OUTPUT-STRING>
 #>     Output string, in which the result is placed
 #> 
-#> 
-#> For more details: https://gdal.org/programs/gdal_vector.html
+#> For more details: <https://gdal.org/programs/gdal_vector.html>
 ```
 
 #### vector clip
@@ -710,6 +774,91 @@ deleteDataset(f_out)
 #> [1] TRUE
 ```
 
+### Pipeline examples
+
+#### Colorized shaded relief
+
+A raster pipline to create a color shaded relief file, based on Example 3 for [GDAL nested pipeline](https://gdal.org/en/latest/programs/gdal_pipeline.html#nested-pipeline).
+
+``` r
+# requires GDAL >= 3.12 for nested pipelines and "raster blend"
+library(gdalraster)
+#> GDAL 3.12.1 (released 2025-12-12), GEOS 3.12.2, PROJ 9.4.1
+
+# gdal_usage("raster pipeline")
+# gdal_usage("raster color-map")
+# gdal_usage("raster hillshade")
+# gdal_usage("raster blend")
+
+# make a working copy the Storm Lake DEM since GDAL will compute stats on it
+# (which creates a .aux.xml sidecar file in the read-only case)
+f <- system.file("extdata/storml_elev.tif", package="gdalraster")
+f_elev <- tempfile(fileext = ".tif")
+file.copy(f, f_elev)
+#> [1] TRUE
+
+# this is a color map file for passing to "raster color-map"
+# see https://gdal.org/programs/gdal_raster_color_map.html
+f_pal <- system.file("extdata/storml_elev_pal.txt", package="gdalraster")
+# writeLines(readLines(f_pal))
+
+# output will be an RGB shaded relief file
+f_out <- file.path(tempdir(), "storml_col_relief.tif")
+
+# arguments for "raster pipeline", using a nested input pipeline for "blend"
+args <- paste(
+    "read --input", f_elev,
+    "! color-map --color-map", f_pal,
+    "! blend --overlay [ read --input", f_elev, "! hillshade -z 1.5 ] --operator=hsv-value",
+    "! write --output", f_out, "--overwrite")
+
+(alg <- gdal_run("raster pipeline", args))
+#> ✔ Done (12ms)
+#> 
+#> C++ object of class <GDALAlg>
+#> • Command: "raster pipeline"
+#> • Description: Process a raster dataset applying several steps.
+#> • Help URL: <https://gdal.org/programs/gdal_raster_pipeline.html>
+
+# note that the "pipeline" algorithm has two outputs: the `output` dataset, and
+# `output_string` (which is empty in this case)
+alg$outputs()
+#> $output
+#> C++ object of class <GDALRaster>
+#>   • Driver: GeoTIFF (GTiff)
+#>   • DSN: "/tmp/Rtmpd9pug9/storml_col_relief.tif"
+#>   • Dimensions: 143, 107, 3
+#>   • CRS: NAD83 / UTM zone 12N (EPSG:26912)
+#>   • Pixel resolution: 30.000000, 30.000000
+#>   • Bbox: 323476.071971, 5101871.983031, 327766.071971, 5105081.983031
+#> 
+#> $output_string
+#> [1] ""
+
+ds <- alg$outputs()$output
+
+alg$release()
+
+plot_raster(ds, main = "Storm Lake AOI colorized shaded relief")
+```
+
+```{r out.width = '75%', echo = FALSE}
+#| fig.alt: >
+#|   A plot of the colorized shaded relief map for the Storm Lake area of
+#|   interest. Lower elevation areas are green, transitioning to yellow then
+#|   light brown, pinkish and almost white at the highest elevations.
+knitr::include_graphics("img/storml_col_shaded_relief.png")
+```
+
+``` r
+# cleanup
+ds$close()
+deleteDataset(f_elev)
+#> [1] TRUE
+deleteDataset(f_out)
+#> [1] TRUE
+```
+
 ## See also
 
 **gdalraster** documentation: