diff --git a/src/raster/r.surf.volcano/Makefile b/src/raster/r.surf.volcano/Makefile
new file mode 100644
index 0000000000..b0e30f319f
--- /dev/null
+++ b/src/raster/r.surf.volcano/Makefile
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.surf.volcano
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
diff --git a/src/raster/r.surf.volcano/r.surf.volcano.html b/src/raster/r.surf.volcano/r.surf.volcano.html
new file mode 100644
index 0000000000..a685ef21fe
--- /dev/null
+++ b/src/raster/r.surf.volcano/r.surf.volcano.html
@@ -0,0 +1,97 @@
+DESCRIPTION
+
+
+r.surf.volcano creates an artificial surface resembling a seamount
+or cone volcano. The user can alter the size and shape of the mountain and
+optionally roughen its surface.
+
+
+NOTES
+
+The friction of distance controls the shape of the mountain when using
+the default polynomial method. Higher values generate steeper slopes.
+
+The pseudo-kurtosis factor is used with all other methods
+to control the slope steepness. For the Gaussian method setting the value
+nearer to zero creates a flatter surface, while higher values generate
+steeper slopes. For Lorentzian, logarithmic, and exponential methods the
+opposite is true.
+
+The surface roughness factor controls the fixed standard deviation
+distance (sigma) used in the Gaussian random number generator.
+It is only used when the -r roughen surface flag is turned on.
+A value closer to zero makes a smoother surface, a higher value makes
+a rougher surface.
+
+It is possible to set a negative value for the peak in order
+to create a pit.
+
+
+
EXAMPLES
+
+Create a simple Gaussian bell:
+
+r.surf.volcano -r output=seamount method=gaussian
+
+# view in the display monitor
+r.colors seamount color=roygbiv
+d.rast seamount
+
+# render in 3D
+m.nviz.image elevation_map=seamount out=gaussian \
+ perspective=10 resolution_fine=1 height=3500
+pnmtopng gaussian.ppm > gaussian.png
+
+# export to Matlab
+r.out.mat in=seamount out=seamount.mat
+
+# integrate into existing DEM
+r.mapcalc "seamount_dem = if(seamount > dem, seamount, dem)"
+r.colors seamount_dem color=srtm
+
+
+
Gaussian bell
+
+r.surf.volcano -r output=volcano crater=250 --verbose
+r.relief in=volcano out=volcano.shade
+
+r.surf.volcano -r output=base_volcano peak=1000 crater=200
+r.surf.fractal output=base_fractal
+r.mapcalc "artificial_land = base_volcano*3.5 + base_fractal"
+
+m.nviz.image elevation_map=artificial_land out=volcano3D \
+ perspective=25 resolution_fine=1 height=30000
+pnmtopng volcano3D.ppm > volcano3D.png
+
+
+
Synthetic volcano from r.surf.volcano combined with fractal
+landscape from r.surf.fractal
+
SEE ALSO
+
+
+r.surf.fractal
+r.surf.gauss
+r.surf.random
+
+
+AUTHOR
+
+Hamish Bowman
+Dept. of Geology
+University of Otago
+Dunedin, New Zealand
diff --git a/src/raster/r.surf.volcano/r.surf.volcano.py b/src/raster/r.surf.volcano/r.surf.volcano.py
new file mode 100755
index 0000000000..3f7eae9eef
--- /dev/null
+++ b/src/raster/r.surf.volcano/r.surf.volcano.py
@@ -0,0 +1,385 @@
+#!/usr/bin/env python3
+############################################################################
+#
+# MODULE: r.surf.volcano
+#
+# AUTHOR: M. Hamish Bowman, Dept. of Geology, University of Otago
+# Dunedin, New Zealand
+# Ported to Python from GRASS 6 addons shell script 8/2024
+#
+# PURPOSE: Create an artificial surface resembling a seamount or cone volcano
+#
+# COPYRIGHT: (c) 2009-2024 Hamish Bowman, and the GRASS Development Team
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+# 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.
+#
+#############################################################################
+
+# %Module
+# % description: Creates an artificial surface resembling a seamount or cone volcano.
+# % keyword: raster
+# %End
+# %Option G_OPT_R_OUTPUT
+# %End
+# %Option
+# % key: peak
+# % type: double
+# % required: no
+# % description: Height of cone
+# % answer: 1000.0
+# %End
+# %Option
+# % key: crater
+# % type: double
+# % required: no
+# % label: Depth of crater below the cone
+# % description: A larger (deeper) value here also means a wider crater.
+# % answer: 0.0
+# %End
+# %Option
+# % key: method
+# % type: string
+# % required: no
+# % description: Mathematical function for creating the mountain
+# % answer: polynomial
+# % options: polynomial,gaussian,lorentzian,logarithmic,exponential
+# % descriptions: polynomial;1/distance^n;gaussian;Gaussian function;lorentzian;Cauchy-Lorentz distribution;logarithmic;Logarithmic decay;exponential;Exponential decay
+# %End
+# %Option
+# % key: friction
+# % type: integer
+# % required: no
+# % options: 1-25
+# % label: Polynomial friction of distance, (the 'n' in 1/d^n)
+# % description: Higher values generate steeper slopes. (only used with the polynomial method)
+# % answer: 6
+# %End
+# # FIXME: ok, it isn't really kurtosis but it's similar and I couldn't
+# # think of a better name.
+# %Option
+# % key: kurtosis
+# % type: double
+# % required: no
+# % label: Slope steepness (used with all methods except polynomial)
+# % description: For Gaussian: nearer to 0 is flatter, higher values generate steeper slopes. For Lorentzian, logarithmic, and exponential the opposite is true.
+# % answer: 1.0
+# %End
+# %Option
+# % key: sigma
+# % type: double
+# % required: no
+# % label: Surface roughness factor
+# % description: Nearer to 0 is smoother, higher values make a rougher surface.
+# % answer: 1.0
+# %End
+# %Flag
+# % key: r
+# % description: Roughen the surface
+# %End
+
+import os
+import sys
+from math import pi
+from grass.script import core as gc
+from grass.script import raster as gr
+
+# from grass.exceptions import CalledModuleError
+
+
+def remove_rast(maps):
+ """Remove raster maps"""
+ gc.run_command("g.remove", flags="f", quiet=True, type="rast", name=maps)
+
+
+def main():
+ outmap = options["output"]
+ method = options["method"]
+ friction = options["friction"]
+ peak = options["peak"]
+ crater = options["crater"]
+ sigma = options["sigma"]
+ kurtosis = options["kurtosis"]
+
+ tmp_base = "tmp__rsv_%d" % os.getpid()
+ map_dist_units = "%s_dist_units" % tmp_base
+ map_dist_norm = "%s_dist_norm" % tmp_base
+ map_peak = "%s_peak" % tmp_base
+ map_surf = "%s_surf" % tmp_base
+
+ gc.verbose(_("Finding cost from center of current region ..."))
+
+ region = gc.region(complete=True)
+
+ gr.mapcalc(
+ "$mdu = sqrt( (x() - $ce)^2 + (y() - $cn)^2 )",
+ mdu=map_dist_units,
+ ce=region["center_easting"],
+ cn=region["center_northing"],
+ )
+
+ gc.verbose(_("Normalizing cost map ..."))
+ rinfo = gr.raster_info(map_dist_units)
+
+ if method == "polynomial":
+ # Normalize with 1 in the center and 0 at outer edge
+ # r.mapcalc "volc.dist_norm.$$ = ($max - volc.dist_units.$$) / $max"
+ gr.mapcalc(
+ "$mdn = ($max - $mdu) / $max",
+ max=rinfo["max"],
+ mdn=map_dist_norm,
+ mdu=map_dist_units,
+ )
+ else:
+ # Normalize with 0 in the center and 1 at outer edge
+ # r.mapcalc "volc.dist_norm.$$ = volc.dist_units.$$ / $max"
+ gr.mapcalc(
+ "$mdn = $mdu / $max",
+ max=rinfo["max"],
+ mdn=map_dist_norm,
+ mdu=map_dist_units,
+ )
+
+ ## ##
+ # create peak map #
+ ## ##
+
+ if method == "polynomial":
+ gc.verbose(_("Creating IDW surface ..."))
+ # r.mapcalc "volc.peak.$$ = ($PEAK + abs($CRATER) ) \
+ # * pow( volc.dist_norm.$$, $FRICTION )"
+ gr.mapcalc(
+ "$mp = ($pk + abs($crat)) * pow($mdn, $fric)",
+ mp=map_peak,
+ pk=peak,
+ crat=crater,
+ mdn=map_dist_norm,
+ fric=friction,
+ )
+ elif method == "gaussian":
+ gc.verbose(_("Creating Gaussian surface ..."))
+ # % description: Use a Gaussian curve instead of 1/(d^n) for radial basis function
+ # normalized Gaussian curve: f(x) = a * e^( (x-b)^2 / 2*(c^2) )
+ # parameters: a = 1/(sigma*sqrt(2pi)), b = mu, c = sigma
+ # mu is mean value and sigma^2 is variance.
+ # so we only need b,c. and b can be locked at 0 here. so user only needs
+ # to give sigma (width)
+ # thus r.mapcalc expr could look like
+ # f(distance) = ( 1 / ($SIGMA*sqrt(2*$PI)) ) * exp( -1* $DIST^2 / (2 * $SIGMA^2) )
+
+ map_gauss = "%s_gauss" % tmp_base
+ SIGMA_C = 1.0
+
+ ## FIXME: the 10*kurtosis stuff is a completely bogus hack!
+ # r.mapcalc "volc.gauss.$$ = \
+ # ( 1 / ( $SIGMA_C * sqrt(2 * $PI) ) ) \
+ # * exp( -1* (10 * $KURTOSIS * volc.dist_norm.$$)^2 / (2 * $SIGMA_C^2) )"
+ gr.mapcalc(
+ "$mg = ( 1 / ( $sigC * sqrt(2 * $pi) ) ) \
+ * exp( -1* (10 * $kt * $mdn)^2 / (2 * $sigC^2) )",
+ mg=map_gauss,
+ sigC=SIGMA_C,
+ pi=pi,
+ kt=kurtosis,
+ mdn=map_dist_norm,
+ )
+
+ rinfo = gr.raster_info(map_gauss)
+ gc.verbose(_("Normalizing Gaussian surface ..."))
+ # r.mapcalc "volc.peak.$$ = \
+ # ( ($PEAK + abs($CRATER) ) / $max ) * volc.gauss.$$"
+ gr.mapcalc(
+ "$mp = ( ($pk + abs($crat) ) / $max ) * $mg",
+ mp=map_peak,
+ pk=peak,
+ crat=crater,
+ max=rinfo["max"],
+ mg=map_gauss,
+ )
+
+ remove_rast(map_gauss)
+ elif method == "lorentzian":
+ # Cauchy-Lorentz fn: f(distance, gamma, height) =
+ # height_of_peak * ( gamma^2 / ( distance^2 + gamma^2) )
+ # where gamma is the scale parameter giving half-width at half-maximum.
+ gc.verbose(_("Creating Lorentzian surface ..."))
+ # r.mapcalc "volc.peak.$$ = ($PEAK + abs($CRATER) ) \
+ # * ( ($KURTOSIS * 0.1)^2 / ( volc.dist_norm.$$ ^2 + ($KURTOSIS * 0.1)^2) )"
+ gr.mapcalc(
+ "$mp = ($pk + abs($crat) ) * ( ($kt * 0.1)^2 / ($mdn^2 + ($kt * 0.1)^2) )",
+ mp=map_peak,
+ pk=peak,
+ crat=crater,
+ kt=kurtosis,
+ mdn=map_dist_norm,
+ )
+ elif method == "exponential":
+ # exponential: 1 / ((e^distance) -1)
+ gc.verbose(_("Creating exponential decay surface ..."))
+ map_exp = "%s_exp" % tmp_base
+
+ # r.mapcalc "volc.exp.$$ = 1 / (exp(volc.dist_norm.$$ / $KURTOSIS) - 0.9)"
+ gr.mapcalc(
+ "$me = 1 / ( exp($mdn / $kt) ) - 0.9",
+ me=map_exp,
+ mdn=map_dist_norm,
+ kt=kurtosis,
+ )
+ rinfo = gr.raster_info(map_exp)
+
+ gc.verbose(_("Normalizing exponential surface ..."))
+ # r.mapcalc "volc.peak.$$ = \
+ # ( ($PEAK + abs($CRATER) ) / $max ) * volc.exp.$$"
+ gr.mapcalc(
+ "$mp = ( ($pk + abs($crat) ) / $max ) * $me",
+ mp=map_peak,
+ pk=peak,
+ crat=crater,
+ max=rinfo["max"],
+ me=map_exp,
+ )
+
+ remove_rast(map_exp)
+ elif method == "logarithmic":
+ # logarithmic: 1 / ( (d+1)^2 * log(d+1) )
+ gc.verbose(_("Creating logarithmic decay surface ..."))
+ map_log = "%s_log" % tmp_base
+
+ # r.mapcalc "volc.log.$$ = 1 / \
+ # ( (volc.dist_norm.$$ + pow(1.15,$KURTOSIS))^2 \
+ # * log((volc.dist_norm.$$) + pow(1.15,$KURTOSIS)) )"
+ gr.mapcalc(
+ "$ml = 1 / ( ($mdn + pow(1.15, $kt))^2 * log(($mdn) + pow(1.15, $kt)) )",
+ ml=map_log,
+ mdn=map_dist_norm,
+ kt=kurtosis,
+ )
+
+ rinfo = gr.raster_info(map_log)
+ gc.verbose(_("Normalizing logarithmic surface ..."))
+ # r.mapcalc "volc.peak.$$ = \
+ # ( ($PEAK + abs($CRATER) ) / $max ) * volc.log.$$"
+ gr.mapcalc(
+ "$mp = ( ($pk + abs($crat) ) / $max ) * $ml",
+ mp=map_peak,
+ pk=peak,
+ crat=crater,
+ max=rinfo["max"],
+ ml=map_log,
+ )
+
+ remove_rast(map_log)
+ else:
+ gc.error("Programmer error, method = <%s>" % method)
+ sys.exit(1)
+
+ if flags["r"]:
+ # roughen it up a bit
+ gc.verbose(_("Creating random Gaussian mottle ..."))
+ map_surf_gauss = "%s_surf_gauss" % tmp_base
+ map_peak_rough = "%s_peak_rough" % tmp_base
+
+ gc.run_command("r.surf.gauss", output=map_surf_gauss, sigma=sigma)
+
+ gc.verbose(_("Applying Gaussian mottle ..."))
+ # r.mapcalc "volc.peak_rough.$$ = \
+ # volc.peak.$$ + (volc.surf_gauss.$$ * $PEAK/400 )"
+ gr.mapcalc(
+ "$vpr = $vp + ($vsg * $pk/400)",
+ vpr=map_peak_rough,
+ vp=map_peak,
+ vsg=map_surf_gauss,
+ pk=peak,
+ )
+
+ gc.run_command("g.rename", raster=(map_peak_rough, map_surf), quiet=True)
+ remove_rast([map_surf_gauss, map_peak])
+
+ else: # no '-r'
+ gc.run_command("g.rename", raster=(map_peak, map_surf), quiet=True)
+
+ remove_rast([map_dist_units, map_dist_norm])
+
+ if float(crater) != 0:
+ gc.verbose(_("Creating crater ..."))
+ map_combo = "%s_full" % tmp_base
+ # r.mapcalc "volc.full.$$ = if( volc.surf.$$ > $PEAK, \
+ # 2*$PEAK - volc.surf.$$, volc.surf.$$ )"
+ gr.mapcalc(
+ "$vc = if($vs > $pk, 2*$pk - $vs, $vs)", vc=map_combo, vs=map_surf, pk=peak
+ )
+
+ gc.run_command("g.rename", raster=(map_combo, outmap), quiet=True)
+ remove_rast(map_surf)
+
+ else:
+ gc.run_command("g.rename", raster=(map_surf, outmap), quiet=True)
+
+ # DCELL is overkill here, convert what we made into FCELL
+ outmapD = "%s_DCELL" % tmp_base
+ gc.run_command("g.rename", raster=(outmap, outmapD), quiet=True)
+ gr.mapcalc("$out = float($outD)", out=outmap, outD=outmapD, quiet=True)
+ remove_rast(outmapD)
+
+ # test if it worked
+ if not gc.find_file(outmap, mapset=".")["name"]:
+ gc.error(_("Surface creation failed"))
+ sys.exit(1)
+
+ # write metadata
+ gc.run_command(
+ "r.support",
+ map=outmap,
+ description="generated by r.surf.volcano",
+ source1="Peak height = %s" % peak,
+ title="Artificial surface resembling a seamount or cone volcano",
+ )
+
+ if flags["r"]:
+ gc.run_command(
+ "r.support",
+ map=outmap,
+ history="Surface roughness used a Gaussian deviate with sigma of %s."
+ % sigma,
+ )
+
+ if float(crater) != 0:
+ gc.run_command("r.support", map=outmap, history="Crater depth %s." % crater)
+
+ if method == "polynomial":
+ s2msg = "Polynomial surface with friction of distance = %s" % friction
+ elif method == "gaussian":
+ s2msg = "Gaussian surface with pseudo-kurtosis factor = %s" % kurtosis
+ elif method == "lorentzian":
+ s2msg = "Lorentzian surface with pseudo-kurtosis factor = %s" % kurtosis
+ elif method == "exponential":
+ s2msg = "Exponential decay surface with pseudo-kurtosis factor = %s" % kurtosis
+ elif method == "logarithmic":
+ s2msg = "Logarithmic decay surface with pseudo-kurtosis factor = %s" % kurtosis
+ else:
+ gc.error("Programmer error, method = <%s>" % method)
+ sys.exit(1)
+
+ gc.run_command("r.support", map=outmap, source2=s2msg)
+
+ # perhaps,
+ # r.colors map=output color=srtm --quiet
+
+ # cleanup on abort/cancel
+ # FIXME: superquiet isn't
+ # TODO: make tmp_base global and move into cleanup_rast() error fn
+ # gc.run_command(
+ # "g.remove", flags="f", superquiet=True, type="rast", patt="%s_*" % tmp_base
+ # )
+
+ gc.verbose(_("Done."))
+
+
+if __name__ == "__main__":
+ options, flags = gc.parser()
+ main()
diff --git a/src/raster/r.surf.volcano/r_surf_volcano_gaussian.jpg b/src/raster/r.surf.volcano/r_surf_volcano_gaussian.jpg
new file mode 100644
index 0000000000..0a08394562
Binary files /dev/null and b/src/raster/r.surf.volcano/r_surf_volcano_gaussian.jpg differ
diff --git a/src/raster/r.surf.volcano/r_surf_volcano_volcano3D.jpg b/src/raster/r.surf.volcano/r_surf_volcano_volcano3D.jpg
new file mode 100644
index 0000000000..9e1c07d8f2
Binary files /dev/null and b/src/raster/r.surf.volcano/r_surf_volcano_volcano3D.jpg differ
diff --git a/src/raster/r.surf.volcano/testsuite/test_rvolc.sh b/src/raster/r.surf.volcano/testsuite/test_rvolc.sh
new file mode 100755
index 0000000000..14bb5e1bea
--- /dev/null
+++ b/src/raster/r.surf.volcano/testsuite/test_rvolc.sh
@@ -0,0 +1,69 @@
+#!/bin/sh
+# Test case for r.surf.volcano: check that resulting raster map
+# is within some arbitrary threshold of matching the map I get.
+# Exits with 0 if everything passed, 1 if something failed.
+# Hamish Bowman Sept. 2024
+#
+# See gunittest_testing.html's basic example for a pretty relevant
+# set of tests. See also r.surf.random's tests for inspiration.
+
+# Tested within the NC sample dataset (nc_spm_08).
+g.region n=5120 s=0 w=0 e=5120 res=10
+
+# touches r.mapcalc, g.region, r.info, g.rename, g.remove
+r.surf.volcano method=gaussian output="volcano.gauss.$$"
+
+if [ $? -ne 0 ] ; then
+ echo "ERROR r.surf.volcano test: did not exit cleanly" 1>&2
+ exit 1
+fi
+
+eval `r.univar -g "volcano.gauss.$$"`
+g.remove -f rast name="volcano.gauss.$$"
+
+
+# reference values
+ref_max=1000
+ref_mean=31.3
+ref_stddev=121.1
+
+# completely arbitrary cross platform variance threshold
+var=5
+
+
+result=`echo "$max $ref_max $var" | awk \
+ 'function abs(x){return ((x < 0) ? -x : x)}
+ {if( abs($1 - $2) < $3 ) {print "ok"}
+ else {print "bad"} }'`
+
+if [ "$result" = "bad" ] ; then
+ echo "ERROR r.surf.volcano test: out of bounds raster maximum" 1>&2
+ exit 1
+fi
+
+
+result=`echo "$mean $ref_mean $var" | awk \
+ 'function abs(x){return ((x < 0) ? -x : x)}
+ {if( abs($1 - $2) < $3 ) {print "ok"}
+ else {print "bad"} }'`
+
+if [ "$result" = "bad" ] ; then
+ echo "ERROR r.surf.volcano test: out of bounds raster mean" 1>&2
+ exit 1
+fi
+
+
+result=`echo "$stddev $ref_stddev $var" | awk \
+ 'function abs(x){return ((x < 0) ? -x : x)}
+ {if( abs($1 - $2) < $3 ) {print "ok"}
+ else {print "bad"} }'`
+
+if [ "$result" = "bad" ] ; then
+ echo "ERROR r.surf.volcano test: out of bounds raster std.dev." 1>&2
+ exit 1
+fi
+
+
+# all good
+exit 0
+