diff --git a/worldeditadditions/init.lua b/worldeditadditions/init.lua
index d391094..ab628ec 100644
--- a/worldeditadditions/init.lua
+++ b/worldeditadditions/init.lua
@@ -23,6 +23,7 @@ wea.io = dofile(wea.modpath.."/utils/io.lua")
 
 wea.bit = dofile(wea.modpath.."/utils/bit.lua")
 
+wea.terrain = dofile(wea.modpath.."/utils/terrain.lua")
 
 dofile(wea.modpath.."/utils/vector.lua")
 dofile(wea.modpath.."/utils/strings/init.lua")
@@ -33,7 +34,7 @@ dofile(wea.modpath.."/utils/tables/init.lua")
 dofile(wea.modpath.."/utils/numbers.lua")
 dofile(wea.modpath.."/utils/nodes.lua")
 dofile(wea.modpath.."/utils/node_identification.lua")
-dofile(wea.modpath.."/utils/terrain.lua")
+
 dofile(wea.modpath.."/utils/raycast_adv.lua") -- For the farwand
 dofile(wea.modpath.."/utils/player.lua") -- Player info functions
 
diff --git a/worldeditadditions/lib/conv/conv.lua b/worldeditadditions/lib/conv/conv.lua
index b374fea..eec69b3 100644
--- a/worldeditadditions/lib/conv/conv.lua
+++ b/worldeditadditions/lib/conv/conv.lua
@@ -1,11 +1,12 @@
-local Vector3 = worldeditadditions.Vector3
+local wea = worldeditadditions
+local Vector3 = wea.Vector3
 
-worldeditadditions.conv = {}
+wea.conv = {}
 
-dofile(worldeditadditions.modpath.."/lib/conv/kernels.lua")
-dofile(worldeditadditions.modpath.."/lib/conv/kernel_gaussian.lua")
+dofile(wea.modpath.."/lib/conv/kernels.lua")
+dofile(wea.modpath.."/lib/conv/kernel_gaussian.lua")
 
-dofile(worldeditadditions.modpath.."/lib/conv/convolve.lua")
+dofile(wea.modpath.."/lib/conv/convolve.lua")
 
 --- Creates a new kernel.
 -- Note that the gaussian kernel only allows for creating a square kernel.
@@ -14,7 +15,7 @@ dofile(worldeditadditions.modpath.."/lib/conv/convolve.lua")
 -- @param	width	number	The width of the kernel to create (must be an odd integer).
 -- @param	height	number	The height of the kernel to create (must be an odd integer).
 -- @param	arg		number	The argument to pass when creating the kernel. Currently only used by gaussian kernel as the sigma value.
-function worldeditadditions.get_conv_kernel(name, width, height, arg)
+function wea.get_conv_kernel(name, width, height, arg)
 	if width % 2 ~= 1 then
 		return false, "Error: The width must be an odd integer.";
 	end
@@ -23,16 +24,16 @@ function worldeditadditions.get_conv_kernel(name, width, height, arg)
 	end
 	
 	if name == "box" then
-		return true, worldeditadditions.conv.kernel_box(width, height)
+		return true, wea.conv.kernel_box(width, height)
 	elseif name == "pascal" then
-		return true, worldeditadditions.conv.kernel_pascal(width, height, true)
+		return true, wea.conv.kernel_pascal(width, height, true)
 	elseif name == "gaussian" then
 		if width ~= height then
 			return false, "Error: When using a gaussian kernel the width and height must be identical."
 		end
 		-- Default to sigma = 2
 		if arg == nil then arg = 2 end
-		local success, result = worldeditadditions.conv.kernel_gaussian(width, arg)
+		local success, result = wea.conv.kernel_gaussian(width, arg)
 		return success, result
 	end
 	
@@ -40,7 +41,7 @@ function worldeditadditions.get_conv_kernel(name, width, height, arg)
 end
 
 
-function worldeditadditions.convolve(pos1, pos2, kernel, kernel_size)
+function wea.convolve(pos1, pos2, kernel, kernel_size)
 	pos1, pos2 = worldedit.sort_pos(pos1, pos2)
 	
 	local border_size = Vector3.new(
@@ -61,10 +62,10 @@ function worldeditadditions.convolve(pos1, pos2, kernel, kernel_size)
 	
 	local node_id_air = minetest.get_content_id("air")
 	
-	local heightmap, heightmap_size = worldeditadditions.make_heightmap(pos1, pos2, manip, area, data)
-	local heightmap_conv = worldeditadditions.table.shallowcopy(heightmap)
+	local heightmap, heightmap_size = wea.terrain.make_heightmap(pos1, pos2, manip, area, data)
+	local heightmap_conv = wea.table.shallowcopy(heightmap)
 	
-	worldeditadditions.conv.convolve(
+	wea.conv.convolve(
 		heightmap_conv,
 		heightmap_size,
 		kernel,
@@ -72,11 +73,11 @@ function worldeditadditions.convolve(pos1, pos2, kernel, kernel_size)
 	)
 	
 	-- print("original")
-	-- worldeditadditions.format.array_2d(heightmap, (pos2.z - pos1.z) + 1)
+	-- wea.format.array_2d(heightmap, (pos2.z - pos1.z) + 1)
 	-- print("transformed")
-	-- worldeditadditions.format.array_2d(heightmap_conv, (pos2.z - pos1.z) + 1)
+	-- wea.format.array_2d(heightmap_conv, (pos2.z - pos1.z) + 1)
 	
-	worldeditadditions.apply_heightmap_changes(
+	wea.terrain.apply_heightmap_changes(
 		pos1, pos2, area, data,
 		heightmap, heightmap_conv, heightmap_size
 	)
diff --git a/worldeditadditions/lib/erode/erode.lua b/worldeditadditions/lib/erode/erode.lua
index 7c263b6..be16a9a 100644
--- a/worldeditadditions/lib/erode/erode.lua
+++ b/worldeditadditions/lib/erode/erode.lua
@@ -1,10 +1,11 @@
-worldeditadditions.erode = {}
+local wea = worldeditadditions
+wea.erode = {}
 
-dofile(worldeditadditions.modpath.."/lib/erode/snowballs.lua")
-dofile(worldeditadditions.modpath.."/lib/erode/river.lua")
+dofile(wea.modpath.."/lib/erode/snowballs.lua")
+dofile(wea.modpath.."/lib/erode/river.lua")
 
 
-function worldeditadditions.erode.run(pos1, pos2, algorithm, params)
+function wea.erode.run(pos1, pos2, algorithm, params)
 	pos1, pos2 = worldedit.sort_pos(pos1, pos2)
 	
 	local manip, area = worldedit.manip_helpers.init(pos1, pos2)
@@ -17,15 +18,15 @@ function worldeditadditions.erode.run(pos1, pos2, algorithm, params)
 	
 	local region_height = (pos2.y - pos1.y) + 1
 	
-	local heightmap = worldeditadditions.make_heightmap(pos1, pos2, manip, area, data)
-	local heightmap_eroded = worldeditadditions.table.shallowcopy(heightmap)
+	local heightmap = wea.terrain.make_heightmap(pos1, pos2, manip, area, data)
+	local heightmap_eroded = wea.table.shallowcopy(heightmap)
 	
 	-- print("[erode.run] algorithm: "..algorithm..", params:");
-	-- print(worldeditadditions.format.map(params))
-	-- worldeditadditions.format.array_2d(heightmap, heightmap_size.x)
+	-- print(wea.format.map(params))
+	-- wea.format.array_2d(heightmap, heightmap_size.x)
 	local success, msg, stats
 	if algorithm == "snowballs" then
-		success, msg = worldeditadditions.erode.snowballs(
+		success, msg = wea.erode.snowballs(
 			heightmap, heightmap_eroded,
 			heightmap_size,
 			region_height,
@@ -33,7 +34,7 @@ function worldeditadditions.erode.run(pos1, pos2, algorithm, params)
 		)
 		if not success then return success, msg end
 	elseif algorithm == "river" then
-		success, msg = worldeditadditions.erode.river(
+		success, msg = wea.erode.river(
 			heightmap, heightmap_eroded,
 			heightmap_size,
 			region_height,
@@ -48,7 +49,7 @@ function worldeditadditions.erode.run(pos1, pos2, algorithm, params)
 		return false, "Error: Unknown algorithm '"..algorithm.."'. Currently implemented algorithms: snowballs (2d; hydraulic-like), river (2d; cellular automata-like; fills potholes and lowers towers). Ideas for algorithms to implement are welcome!"
 	end
 	
-	success, stats = worldeditadditions.apply_heightmap_changes(
+	success, stats = wea.terrain.apply_heightmap_changes(
 		pos1, pos2, area, data,
 		heightmap, heightmap_eroded, heightmap_size
 	)
diff --git a/worldeditadditions/lib/erode/snowballs.lua b/worldeditadditions/lib/erode/snowballs.lua
index 29853a5..cc10ba5 100644
--- a/worldeditadditions/lib/erode/snowballs.lua
+++ b/worldeditadditions/lib/erode/snowballs.lua
@@ -1,4 +1,5 @@
-local Vector3 = worldeditadditions.Vector3
+local wea = worldeditadditions
+local Vector3 = wea.Vector3
 
 -- Test command: //multi //fp set1 1313 6 5540 //fp set2 1338 17 5521 //erode snowballs
 
@@ -28,7 +29,7 @@ local function snowball(heightmap, normalmap, heightmap_size, startpos, params)
 		end
 		
 		if #hist_velocity > 0 and i > 5
-			and worldeditadditions.average(hist_velocity) < 0.03 then
+			and wea.average(hist_velocity) < 0.03 then
 			-- print("[snowball] It looks like we've stopped")
 			return true, i
 		end
@@ -51,11 +52,11 @@ local function snowball(heightmap, normalmap, heightmap_size, startpos, params)
 		velocity.x = params.friction * velocity.x + normalmap[hi].x * params.speed
 		velocity.z = params.friction * velocity.z + normalmap[hi].y * params.speed
 		
-		-- print("[snowball] now at ("..x..", "..z..") velocity "..worldeditadditions.vector.lengthsquared(velocity)..", sediment "..sediment)
-		local new_vel_sq = worldeditadditions.vector.lengthsquared(velocity)
+		-- print("[snowball] now at ("..x..", "..z..") velocity "..wea.vector.lengthsquared(velocity)..", sediment "..sediment)
+		local new_vel_sq = wea.vector.lengthsquared(velocity)
 		if new_vel_sq > 1 then
 			-- print("[snowball] velocity squared over 1, normalising")
-			velocity = worldeditadditions.vector.normalize(velocity)
+			velocity = wea.vector.normalize(velocity)
 		end
 		table.insert(hist_velocity, new_vel_sq)
 		if #hist_velocity > params.velocity_hist_count then table.remove(hist_velocity, 1) end
@@ -74,7 +75,7 @@ Note that this *mutates* the given heightmap.
 @source https://jobtalle.com/simulating_hydraulic_erosion.html
 
 ]]--
-function worldeditadditions.erode.snowballs(heightmap_initial, heightmap, heightmap_size, region_height, params_custom)
+function wea.erode.snowballs(heightmap_initial, heightmap, heightmap_size, region_height, params_custom)
 	local params = {
 		rate_deposit = 0.03, -- 0.03
 		rate_erosion = 0.04, -- 0.04
@@ -88,12 +89,12 @@ function worldeditadditions.erode.snowballs(heightmap_initial, heightmap, height
 		count = 25000
 	}
 	-- Apply the custom settings
-	worldeditadditions.table.apply(params_custom, params)
+	wea.table.apply(params_custom, params)
 	
 	-- print("[erode/snowballs] params: ")
-	-- print(worldeditadditions.format.map(params))
+	-- print(wea.format.map(params))
 	
-	local normals = worldeditadditions.calculate_normals(heightmap, heightmap_size)
+	local normals = wea.terrain.calculate_normals(heightmap, heightmap_size)
 	
 	local stats_steps = {}
 	for i = 1, params.count do
@@ -111,7 +112,7 @@ function worldeditadditions.erode.snowballs(heightmap_initial, heightmap, height
 		if not success then return false, "Error: Failed at snowball "..i..":"..steps end
 	end
 	
-	-- print("[snowballs] "..#stats_steps.." snowballs simulated, max "..params.max_steps.." steps, averaged ~"..worldeditadditions.average(stats_steps).."")
+	-- print("[snowballs] "..#stats_steps.." snowballs simulated, max "..params.max_steps.." steps, averaged ~"..wea.average(stats_steps).."")
 	
 	-- Round everything to the nearest int, since you can't really have
 	-- something like .141592671 of a node
@@ -131,15 +132,15 @@ function worldeditadditions.erode.snowballs(heightmap_initial, heightmap, height
 	end
 	
 	if not params.noconv then
-		local success, matrix = worldeditadditions.get_conv_kernel("gaussian", 3, 3)
+		local success, matrix = wea.get_conv_kernel("gaussian", 3, 3)
 		if not success then return success, matrix end
 		local matrix_size = Vector3.new(3, 0, 3)
-		worldeditadditions.conv.convolve(
+		wea.conv.convolve(
 			heightmap, heightmap_size,
 			matrix,
 			matrix_size
 		)
 	end
 	
-	return true, ""..#stats_steps.." snowballs simulated, max "..params.max_steps.." steps (averaged ~"..worldeditadditions.average(stats_steps).." steps)"
+	return true, ""..#stats_steps.." snowballs simulated, max "..params.max_steps.." steps (averaged ~"..wea.average(stats_steps).." steps)"
 end
diff --git a/worldeditadditions/lib/layers.lua b/worldeditadditions/lib/layers.lua
index 80fef4e..8ace75f 100644
--- a/worldeditadditions/lib/layers.lua
+++ b/worldeditadditions/lib/layers.lua
@@ -35,11 +35,11 @@ function worldeditadditions.layers(pos1, pos2, node_weights, min_slope, max_slop
 	
 	local node_ids, node_ids_count = wea.unwind_node_list(node_weights)
 	
-	local heightmap, heightmap_size = wea.make_heightmap(
+	local heightmap, heightmap_size = wea.terrain.make_heightmap(
 		pos1, pos2,
 		manip, area, data
 	)
-	local slopemap = wea.calculate_slopes(heightmap, heightmap_size)
+	local slopemap = wea.terrain.calculate_slopes(heightmap, heightmap_size)
 	-- worldeditadditions.format.array_2d(heightmap, heightmap_size.x)
 	-- print_slopes(slopemap, heightmap_size.x)
 	--luacheck:ignore 311
diff --git a/worldeditadditions/lib/noise/apply_2d.lua b/worldeditadditions/lib/noise/apply_2d.lua
index d0cd7eb..1885871 100644
--- a/worldeditadditions/lib/noise/apply_2d.lua
+++ b/worldeditadditions/lib/noise/apply_2d.lua
@@ -3,13 +3,13 @@ local wea = worldeditadditions
 
 --- Applies the given noise field to the given heightmap.
 -- Mutates the given heightmap.
--- @param	heightmap		number[]			A table of ZERO indexed numbers representing the heghtmap - see worldeditadditions.make_heightmap().
+-- @param	heightmap		number[]			A table of ZERO indexed numbers representing the heghtmap - see worldeditadditions.terrain.make_heightmap().
 -- @param	noise			number[]			An table identical in structure to the heightmap containing the noise values to apply.
 -- @param	heightmap_size	{x:number,z:number}	A 2d vector representing the size of the heightmap.
 -- @param	region_height	number				The height of the defined region.
 -- @param	apply_mode		string				The apply mode to use to apply the noise to the heightmap.
 -- @returns	bool[,string]	A boolean value representing whether the application was successful or not. If false, then an error message as a string is also returned describing the error that occurred.
-function worldeditadditions.noise.apply_2d(heightmap, noise, heightmap_size, pos1, pos2, apply_mode)
+function wea.noise.apply_2d(heightmap, noise, heightmap_size, pos1, pos2, apply_mode)
 	if type(apply_mode) ~= "string" and type(apply_mode) ~= "number" then
 		return false, "Error: Expected value of type string or number for apply_mode, but received value of type "..type(apply_mode)
 	end
@@ -17,12 +17,12 @@ function worldeditadditions.noise.apply_2d(heightmap, noise, heightmap_size, pos
 	local region_height = pos2.y - pos1.y
 	
 	
-	print("NOISE APPLY_2D\n")
-	worldeditadditions.format.array_2d(noise, heightmap_size.x)
+	-- print("NOISE APPLY_2D\n")
+	wea.format.array_2d(noise, heightmap_size.x)
 	
 	
 	local height = tonumber(apply_mode)
-	print("DEBUG apply_mode", apply_mode, "as height", height)
+	-- print("DEBUG apply_mode", apply_mode, "as height", height)
 	
 	for z = heightmap_size.z - 1, 0, -1 do
 		for x = heightmap_size.x - 1, 0, -1 do
@@ -53,8 +53,8 @@ function worldeditadditions.noise.apply_2d(heightmap, noise, heightmap_size, pos
 	-- 	heightmap[(z * heightmap_size.x) + x] = z
 	-- end
 	
-	print("HEIGHTMAP\n")
-	worldeditadditions.format.array_2d(heightmap, heightmap_size.x)
+	-- print("HEIGHTMAP\n")
+	-- worldeditadditions.format.array_2d(heightmap, heightmap_size.x)
 	
 	
 	return true
diff --git a/worldeditadditions/lib/noise/run2d.lua b/worldeditadditions/lib/noise/run2d.lua
index 001f8ce..c08657a 100644
--- a/worldeditadditions/lib/noise/run2d.lua
+++ b/worldeditadditions/lib/noise/run2d.lua
@@ -12,25 +12,25 @@ local wea = worldeditadditions
 -- @param	pos1			Vector	pos1 of the defined region
 -- @param	pos2			Vector	pos2 of the defined region
 -- @param	noise_params	table	A noise parameters table.
-function worldeditadditions.noise.run2d(pos1, pos2, noise_params)
+function wea.noise.run2d(pos1, pos2, noise_params)
 	pos1, pos2 = worldedit.sort_pos(pos1, pos2)
 	-- pos2 will always have the highest co-ordinates now
 	
 	-- Fill in the default params
 	-- print("DEBUG noise_params_custom ", wea.format.map(noise_params))
-	noise_params = worldeditadditions.noise.params_apply_default(noise_params)
+	noise_params = wea.noise.params_apply_default(noise_params)
 	-- print("DEBUG noise_params[1] ", wea.format.map(noise_params[1]))
 	
 	-- Fetch the nodes in the specified area
 	local manip, area = worldedit.manip_helpers.init(pos1, pos2)
 	local data = manip:get_data()
 	
-	local heightmap_old, heightmap_size = worldeditadditions.make_heightmap(
+	local heightmap_old, heightmap_size = wea.terrain.make_heightmap(
 		pos1, pos2,
 		manip, area,
 		data
 	)
-	local heightmap_new = worldeditadditions.table.shallowcopy(heightmap_old)
+	local heightmap_new = wea.table.shallowcopy(heightmap_old)
 	
 	local success, noisemap = wea.noise.make_2d(
 		heightmap_size,
@@ -49,7 +49,7 @@ function worldeditadditions.noise.run2d(pos1, pos2, noise_params)
 	if not success then return success, message end
 	
 	local stats
-	success, stats = wea.apply_heightmap_changes(
+	success, stats = wea.terrain.apply_heightmap_changes(
 		pos1, pos2,
 		area, data,
 		heightmap_old, heightmap_new,
diff --git a/worldeditadditions/lib/sculpt/apply.lua b/worldeditadditions/lib/sculpt/apply.lua
index 6887c0f..c2da299 100644
--- a/worldeditadditions/lib/sculpt/apply.lua
+++ b/worldeditadditions/lib/sculpt/apply.lua
@@ -36,7 +36,7 @@ local function apply(pos1, brush_name, height, brush_size)
 	local manip, area = worldedit.manip_helpers.init(pos1_compute, pos2_compute)
 	local data = manip:get_data()
 	
-	local heightmap, heightmap_size = wea.make_heightmap(
+	local heightmap, heightmap_size = wea.terrain.make_heightmap(
 		pos1_compute, pos2_compute,
 		manip, area,
 		data
@@ -52,7 +52,7 @@ local function apply(pos1, brush_name, height, brush_size)
 	if not success2 then return success2, added end
 	
 	-- 3: Save back to disk & return
-	local success3, stats = wea.apply_heightmap_changes(
+	local success3, stats = wea.terrain.apply_heightmap_changes(
 		pos1_compute, pos2_compute,
 		area, data,
 		heightmap_orig, heightmap,
diff --git a/worldeditadditions/utils/terrain.lua b/worldeditadditions/utils/terrain.lua
deleted file mode 100644
index 8649871..0000000
--- a/worldeditadditions/utils/terrain.lua
+++ /dev/null
@@ -1,180 +0,0 @@
-local wea = worldeditadditions
-local Vector3 = wea.Vector3
-
---- Given a manip object and associates, generates a 2D x/z heightmap.
--- Note that pos1 and pos2 should have already been pushed through
--- worldedit.sort_pos(pos1, pos2) before passing them to this function.
--- @param	pos1	Vector		Position 1 of the region to operate on
--- @param	pos2	Vector		Position 2 of the region to operate on
--- @param	manip	VoxelManip	The VoxelManip object.
--- @param	area	area		The associated area object.
--- @param	data	table		The associated data object.
--- @return	table,table			The ZERO-indexed heightmap data (as 1 single flat array), followed by the size of the heightmap in the form { z = size_z, x = size_x }.
-function worldeditadditions.make_heightmap(pos1, pos2, manip, area, data)
-	-- z y x (in reverse for little-endian machines) is the preferred loop order, but that isn't really possible here
-	
-	local heightmap = {}
-	local hi = 0
-	local changes = { updated = 0, skipped_columns = 0 }
-	for z = pos1.z, pos2.z, 1 do
-		for x = pos1.x, pos2.x, 1 do
-			local found_node = false
-			-- Scan each column top to bottom
-			for y = pos2.y+1, pos1.y, -1 do
-				local i = area:index(x, y, z)
-				if not (wea.is_airlike(data[i]) or wea.is_liquidlike(data[i])) then
-					-- It's the first non-airlike node in this column
-					-- Start heightmap values from 1 (i.e. there's at least 1 node in the column)
-					heightmap[hi] = (y - pos1.y) + 1 
-					found_node = true
-					break
-				end
-			end
-			
-			if not found_node then heightmap[hi] = -1 end
-			hi = hi + 1
-		end
-	end
-	
-	local heightmap_size = Vector3.new(
-		(pos2.x - pos1.x) + 1,	-- x
-		0,						-- y
-		(pos2.z - pos1.z) + 1	-- z
-	)
-	
-	return heightmap, heightmap_size
-end
-
---- Calculates a normal map for the given heightmap.
--- Caution: This method (like worldeditadditions.make_heightmap) works on
--- X AND Z, and NOT x and y. This means that the resulting 3d normal vectors
--- will have the z and y values swapped.
--- @param	heightmap		table	A ZERO indexed flat heightmap. See worldeditadditions.make_heightmap().
--- @param	heightmap_size	int[]	The size of the heightmap in the form [ z, x ]
--- @return	Vector[]		The calculated normal map, in the same form as the input heightmap. Each element of the array is a Vector3 instance representing a normal.
-function worldeditadditions.calculate_normals(heightmap, heightmap_size)
-	-- print("heightmap_size: "..heightmap_size.x.."x"..heightmap_size.z)
-	local result = {}
-	for z = heightmap_size.z-1, 0, -1 do
-		for x = heightmap_size.x-1, 0, -1 do
-			-- Algorithm ref https://stackoverflow.com/a/13983431/1460422
-			-- Also ref Vector.mjs, which I implemented myself (available upon request)
-			local hi = z*heightmap_size.x + x
-			-- Default to this pixel's height
-			local up = heightmap[hi]
-			local down = heightmap[hi]
-			local left = heightmap[hi]
-			local right = heightmap[hi]
-			if z - 1 > 0 then up = heightmap[(z-1)*heightmap_size.x + x] end
-			if z + 1 < heightmap_size.z-1 then down = heightmap[(z+1)*heightmap_size.x + x] end
-			if x - 1 > 0 then left = heightmap[z*heightmap_size.x + (x-1)] end
-			if x + 1 < heightmap_size.x-1 then right = heightmap[z*heightmap_size.x + (x+1)] end
-			
-			-- print("[normals] UP	| index", (z-1)*heightmap_size.x + x, "z", z, "z-1", z - 1, "up", up, "limit", 0)
-			-- print("[normals] DOWN	| index", (z+1)*heightmap_size.x + x, "z", z, "z+1", z + 1, "down", down, "limit", heightmap_size.x-1)
-			-- print("[normals] LEFT	| index", z*heightmap_size.x + (x-1), "x", x, "x-1", x - 1, "left", left, "limit", 0)
-			-- print("[normals] RIGHT	| index", z*heightmap_size.x + (x+1), "x", x, "x+1", x + 1, "right", right, "limit", heightmap_size.x-1)
-			
-			result[hi] = wea.Vector3.new(
-				left - right,	-- x
-				2,				-- y - Z & Y are flipped
-				down - up		-- z
-			):normalise()
-			
-			-- print("[normals] at "..hi.." ("..x..", "..z..") normal "..result[hi])
-		end
-	end
-	return result
-end
-
---- Converts a 2d heightmap into slope values in radians.
--- Convert a radians to degrees by doing (radians*math.pi) / 180 for display,
--- but it is STRONGLY recommended to keep all internal calculations in radians.
--- @param	heightmap		table	A ZERO indexed flat heightmap. See worldeditadditions.make_heightmap().
--- @param	heightmap_size	int[]	The size of the heightmap in the form [ z, x ]
--- @return	Vector[]		The calculated slope map, in the same form as the input heightmap. Each element of the array is a (floating-point) number representing the slope in that cell in radians.
-function worldeditadditions.calculate_slopes(heightmap, heightmap_size)
-	local normals = worldeditadditions.calculate_normals(heightmap, heightmap_size)
-	local slopes = {  }
-	
-	local up = wea.Vector3.new(0, 1, 0) -- Z & Y are flipped
-	
-	for z = heightmap_size.z-1, 0, -1 do
-		for x = heightmap_size.x-1, 0, -1 do
-			local hi = z*heightmap_size.x + x
-			
-			-- Ref https://stackoverflow.com/a/16669463/1460422
-			-- slopes[hi] = wea.Vector3.dot_product(normals[hi], up)
-			slopes[hi] = math.acos(normals[hi].y)
-		end
-	end
-	
-	return slopes
-end
-
---- Applies changes to a heightmap to a Voxel Manipulator data block.
--- @param	pos1	vector		Position 1 of the defined region
--- @param	pos2	vector		Position 2 of the defined region
--- @param	area	VoxelArea	The VoxelArea object (see worldedit.manip_helpers.init)
--- @param	data	number[]	The node ids data array containing the slice of the Minetest world extracted using the Voxel Manipulator.
--- @param	heightmap_old	number[]	The original heightmap from worldeditadditions.make_heightmap.
--- @param	heightmap_new	number[]	The new heightmap containing the altered updated values. It is expected that worldeditadditions.table.shallowcopy be used to make a COPY of the data worldeditadditions.make_heightmap for this purpose. Both heightmap_old AND heightmap_new are REQUIRED in order for this function to work.
--- @param	heightmap_size	vector	The x / z size of the heightmap. Any y value set in the vector is ignored.
--- @returns	bool, string|{ added: number, removed: number }	A bool indicating whether the operation was successful or not, followed by either an error message as a string (if it was not successful) or a table of statistics (if it was successful).
-function worldeditadditions.apply_heightmap_changes(pos1, pos2, area, data, heightmap_old, heightmap_new, heightmap_size)
-	local stats = { added = 0, removed = 0 }
-	local node_id_air = minetest.get_content_id("air")
-	local node_id_ignore = minetest.get_content_id("ignore")
-	
-	for z = heightmap_size.z - 1, 0, -1 do
-		for x = heightmap_size.x - 1, 0, -1 do
-			local hi = z*heightmap_size.x + x
-			
-			local height_old = heightmap_old[hi]
-			local height_new = heightmap_new[hi]
-			-- print("[conv/save] hi", hi, "height_old", heightmap_old[hi], "height_new", heightmap_new[hi], "z", z, "x", x, "pos1.y", pos1.y)
-			
-			-- Lua doesn't have a continue statement :-/
-			if height_old == height_new then
-				-- noop
-			elseif height_new < height_old then
-				local node_id_replace = data[area:index(
-					pos1.x + x,
-					pos1.y + height_old + 1,
-					pos1.z + z
-				)]
-				-- Unlikely, but if it can happen, it *will* happen.....
-				if node_id_replace == node_id_ignore then
-					node_id_replace = node_id_air
-				end
-				stats.removed = stats.removed + (height_old - height_new)
-				local y = height_new
-				while y < height_old do
-					local ci = area:index(pos1.x + x, pos1.y + y, pos1.z + z)
-					-- print("[conv/save] remove at y", y, "→", pos1.y + y, "current:", minetest.get_name_from_content_id(data[ci]))
-					if data[ci] ~= node_id_ignore then
-						data[ci] = node_id_replace
-					end
-					y = y + 1
-				end
-			else -- height_new > height_old
-				-- We subtract one because the heightmap starts at 1 (i.e. 1 = 1 node in the column), but the selected region is inclusive
-				local node_id = data[area:index(pos1.x + x, pos1.y + (height_old - 1), pos1.z + z)]
-				-- print("[conv/save] filling with ", node_id, "→", minetest.get_name_from_content_id(node_id))
-				
-				stats.added = stats.added + (height_new - height_old)
-				local y = height_old
-				while y < height_new do
-					local ci = area:index(pos1.x + x, pos1.y + y, pos1.z + z)
-					-- print("[conv/save] add at y", y, "→", pos1.y + y, "current:", minetest.get_name_from_content_id(data[ci]))
-					if data[ci] ~= node_id_ignore then
-						data[ci] = node_id
-					end
-					y = y + 1
-				end
-			end
-		end
-	end
-	
-	return true, stats
-end
diff --git a/worldeditadditions/utils/terrain/apply_heightmap_changes.lua b/worldeditadditions/utils/terrain/apply_heightmap_changes.lua
new file mode 100644
index 0000000..fc1e1d1
--- /dev/null
+++ b/worldeditadditions/utils/terrain/apply_heightmap_changes.lua
@@ -0,0 +1,71 @@
+local wea = worldeditadditions
+local Vector3 = wea.Vector3
+
+--- Applies changes to a heightmap to a Voxel Manipulator data block.
+-- @param	pos1	vector		Position 1 of the defined region
+-- @param	pos2	vector		Position 2 of the defined region
+-- @param	area	VoxelArea	The VoxelArea object (see worldedit.manip_helpers.init)
+-- @param	data	number[]	The node ids data array containing the slice of the Minetest world extracted using the Voxel Manipulator.
+-- @param	heightmap_old	number[]	The original heightmap from worldeditadditions.make_heightmap.
+-- @param	heightmap_new	number[]	The new heightmap containing the altered updated values. It is expected that worldeditadditions.table.shallowcopy be used to make a COPY of the data worldeditadditions.make_heightmap for this purpose. Both heightmap_old AND heightmap_new are REQUIRED in order for this function to work.
+-- @param	heightmap_size	vector	The x / z size of the heightmap. Any y value set in the vector is ignored.
+-- @returns	bool, string|{ added: number, removed: number }	A bool indicating whether the operation was successful or not, followed by either an error message as a string (if it was not successful) or a table of statistics (if it was successful).
+local function apply_heightmap_changes(pos1, pos2, area, data, heightmap_old, heightmap_new, heightmap_size)
+	local stats = { added = 0, removed = 0 }
+	local node_id_air = minetest.get_content_id("air")
+	local node_id_ignore = minetest.get_content_id("ignore")
+	
+	for z = heightmap_size.z - 1, 0, -1 do
+		for x = heightmap_size.x - 1, 0, -1 do
+			local hi = z*heightmap_size.x + x
+			
+			local height_old = heightmap_old[hi]
+			local height_new = heightmap_new[hi]
+			-- print("[conv/save] hi", hi, "height_old", heightmap_old[hi], "height_new", heightmap_new[hi], "z", z, "x", x, "pos1.y", pos1.y)
+			
+			-- Lua doesn't have a continue statement :-/
+			if height_old == height_new then
+				-- noop
+			elseif height_new < height_old then
+				local node_id_replace = data[area:index(
+					pos1.x + x,
+					pos1.y + height_old + 1,
+					pos1.z + z
+				)]
+				-- Unlikely, but if it can happen, it *will* happen.....
+				if node_id_replace == node_id_ignore then
+					node_id_replace = node_id_air
+				end
+				stats.removed = stats.removed + (height_old - height_new)
+				local y = height_new
+				while y < height_old do
+					local ci = area:index(pos1.x + x, pos1.y + y, pos1.z + z)
+					-- print("[conv/save] remove at y", y, "→", pos1.y + y, "current:", minetest.get_name_from_content_id(data[ci]))
+					if data[ci] ~= node_id_ignore then
+						data[ci] = node_id_replace
+					end
+					y = y + 1
+				end
+			else -- height_new > height_old
+				-- We subtract one because the heightmap starts at 1 (i.e. 1 = 1 node in the column), but the selected region is inclusive
+				local node_id = data[area:index(pos1.x + x, pos1.y + (height_old - 1), pos1.z + z)]
+				-- print("[conv/save] filling with ", node_id, "→", minetest.get_name_from_content_id(node_id))
+				
+				stats.added = stats.added + (height_new - height_old)
+				local y = height_old
+				while y < height_new do
+					local ci = area:index(pos1.x + x, pos1.y + y, pos1.z + z)
+					-- print("[conv/save] add at y", y, "→", pos1.y + y, "current:", minetest.get_name_from_content_id(data[ci]))
+					if data[ci] ~= node_id_ignore then
+						data[ci] = node_id
+					end
+					y = y + 1
+				end
+			end
+		end
+	end
+	
+	return true, stats
+end
+
+return apply_heightmap_changes
diff --git a/worldeditadditions/utils/terrain/calculate_normals.lua b/worldeditadditions/utils/terrain/calculate_normals.lua
new file mode 100644
index 0000000..fbadf11
--- /dev/null
+++ b/worldeditadditions/utils/terrain/calculate_normals.lua
@@ -0,0 +1,47 @@
+local wea = worldeditadditions
+local Vector3 = wea.Vector3
+
+
+--- Calculates a normal map for the given heightmap.
+-- Caution: This method (like worldeditadditions.make_heightmap) works on
+-- X AND Z, and NOT x and y. This means that the resulting 3d normal vectors
+-- will have the z and y values swapped.
+-- @param	heightmap		table	A ZERO indexed flat heightmap. See worldeditadditions.terrain.make_heightmap().
+-- @param	heightmap_size	int[]	The size of the heightmap in the form [ z, x ]
+-- @return	Vector[]		The calculated normal map, in the same form as the input heightmap. Each element of the array is a Vector3 instance representing a normal.
+local function calculate_normals(heightmap, heightmap_size)
+	-- print("heightmap_size: "..heightmap_size.x.."x"..heightmap_size.z)
+	local result = {}
+	for z = heightmap_size.z-1, 0, -1 do
+		for x = heightmap_size.x-1, 0, -1 do
+			-- Algorithm ref https://stackoverflow.com/a/13983431/1460422
+			-- Also ref Vector.mjs, which I implemented myself (available upon request)
+			local hi = z*heightmap_size.x + x
+			-- Default to this pixel's height
+			local up = heightmap[hi]
+			local down = heightmap[hi]
+			local left = heightmap[hi]
+			local right = heightmap[hi]
+			if z - 1 > 0 then up = heightmap[(z-1)*heightmap_size.x + x] end
+			if z + 1 < heightmap_size.z-1 then down = heightmap[(z+1)*heightmap_size.x + x] end
+			if x - 1 > 0 then left = heightmap[z*heightmap_size.x + (x-1)] end
+			if x + 1 < heightmap_size.x-1 then right = heightmap[z*heightmap_size.x + (x+1)] end
+			
+			-- print("[normals] UP	| index", (z-1)*heightmap_size.x + x, "z", z, "z-1", z - 1, "up", up, "limit", 0)
+			-- print("[normals] DOWN	| index", (z+1)*heightmap_size.x + x, "z", z, "z+1", z + 1, "down", down, "limit", heightmap_size.x-1)
+			-- print("[normals] LEFT	| index", z*heightmap_size.x + (x-1), "x", x, "x-1", x - 1, "left", left, "limit", 0)
+			-- print("[normals] RIGHT	| index", z*heightmap_size.x + (x+1), "x", x, "x+1", x + 1, "right", right, "limit", heightmap_size.x-1)
+			
+			result[hi] = Vector3.new(
+				left - right,	-- x
+				2,				-- y - Z & Y are flipped
+				down - up		-- z
+			):normalise()
+			
+			-- print("[normals] at "..hi.." ("..x..", "..z..") normal "..result[hi])
+		end
+	end
+	return result
+end
+
+return calculate_normals
diff --git a/worldeditadditions/utils/terrain/calculate_slopes.lua b/worldeditadditions/utils/terrain/calculate_slopes.lua
new file mode 100644
index 0000000..8da99eb
--- /dev/null
+++ b/worldeditadditions/utils/terrain/calculate_slopes.lua
@@ -0,0 +1,30 @@
+local wea = worldeditadditions
+local Vector3 = wea.Vector3
+
+
+--- Converts a 2d heightmap into slope values in radians.
+-- Convert a radians to degrees by doing (radians*math.pi) / 180 for display,
+-- but it is STRONGLY recommended to keep all internal calculations in radians.
+-- @param	heightmap		table	A ZERO indexed flat heightmap. See worldeditadditions.terrain.make_heightmap().
+-- @param	heightmap_size	int[]	The size of the heightmap in the form [ z, x ]
+-- @return	Vector[]		The calculated slope map, in the same form as the input heightmap. Each element of the array is a (floating-point) number representing the slope in that cell in radians.
+local function calculate_slopes(heightmap, heightmap_size)
+	local normals = wea.terrain.calculate_normals(heightmap, heightmap_size)
+	local slopes = {  }
+	
+	local up = wea.Vector3.new(0, 1, 0) -- Z & Y are flipped
+	
+	for z = heightmap_size.z-1, 0, -1 do
+		for x = heightmap_size.x-1, 0, -1 do
+			local hi = z*heightmap_size.x + x
+			
+			-- Ref https://stackoverflow.com/a/16669463/1460422
+			-- slopes[hi] = wea.Vector3.dot_product(normals[hi], up)
+			slopes[hi] = math.acos(normals[hi].y)
+		end
+	end
+	
+	return slopes
+end
+
+return calculate_slopes
diff --git a/worldeditadditions/utils/terrain/init.lua b/worldeditadditions/utils/terrain/init.lua
new file mode 100644
index 0000000..c83577e
--- /dev/null
+++ b/worldeditadditions/utils/terrain/init.lua
@@ -0,0 +1,12 @@
+local wea = worldeditadditions
+local Vector3 = wea.Vector3
+
+
+local terrain = {
+	make_heightmap = dofile(wea.modpath.."/utils/terrain/make_heightmap.lua"),
+	calculate_normals = dofile(wea.modpath.."/utils/terrain/calculate_normals.lua"),
+	calculate_slopes = dofile(wea.modpath.."/utils/terrain/calculate_slopes.lua"),
+	apply_heightmap_changes = dofile(wea.modpath.."/utils/terrain/apply_heightmap_changes.lua")
+}
+
+return terrain
diff --git a/worldeditadditions/utils/terrain/make_heightmap.lua b/worldeditadditions/utils/terrain/make_heightmap.lua
new file mode 100644
index 0000000..1b3c67b
--- /dev/null
+++ b/worldeditadditions/utils/terrain/make_heightmap.lua
@@ -0,0 +1,50 @@
+local wea = worldeditadditions
+local Vector3 = wea.Vector3
+
+
+--- Given a manip object and associates, generates a 2D x/z heightmap.
+-- Note that pos1 and pos2 should have already been pushed through
+-- worldedit.sort_pos(pos1, pos2) before passing them to this function.
+-- @param	pos1	Vector		Position 1 of the region to operate on
+-- @param	pos2	Vector		Position 2 of the region to operate on
+-- @param	manip	VoxelManip	The VoxelManip object.
+-- @param	area	area		The associated area object.
+-- @param	data	table		The associated data object.
+-- @return	table,table			The ZERO-indexed heightmap data (as 1 single flat array), followed by the size of the heightmap in the form { z = size_z, x = size_x }.
+local function make_heightmap(pos1, pos2, manip, area, data)
+	-- z y x (in reverse for little-endian machines) is the preferred loop order, but that isn't really possible here
+	
+	local heightmap = {}
+	local hi = 0
+	local changes = { updated = 0, skipped_columns = 0 }
+	for z = pos1.z, pos2.z, 1 do
+		for x = pos1.x, pos2.x, 1 do
+			local found_node = false
+			-- Scan each column top to bottom
+			for y = pos2.y+1, pos1.y, -1 do
+				local i = area:index(x, y, z)
+				if not (wea.is_airlike(data[i]) or wea.is_liquidlike(data[i])) then
+					-- It's the first non-airlike node in this column
+					-- Start heightmap values from 1 (i.e. there's at least 1 node in the column)
+					heightmap[hi] = (y - pos1.y) + 1 
+					found_node = true
+					break
+				end
+			end
+			
+			if not found_node then heightmap[hi] = -1 end
+			hi = hi + 1
+		end
+	end
+	
+	local heightmap_size = Vector3.new(
+		(pos2.x - pos1.x) + 1,	-- x
+		0,						-- y
+		(pos2.z - pos1.z) + 1	-- z
+	)
+	
+	return heightmap, heightmap_size
+end
+
+
+return make_heightmap