How to Create a Smooth Heightmap

heightmap_smooth.py

# https://www.reddit.com/r/proceduralgeneration/comments/d5gde4/heightmap_with_strict_neighbour_cell_constraints/
import numpy as np
import scipy.ndimage as ndimage

def generate(w, h, d):
    return np.random.randint(0, d, (w, h))

def smooth(arr):
    for z in range(np.min(arr), np.max(arr) + 1):
        # Set of all tiles at height Z.
        ztiles = arr == z
        # Find tiles adjacent to tiles at height Z.
        adjacent = ndimage.binary_dilation(ztiles)
        # Set them to a maximum height of Z + 1.
        arr[adjacent] = np.minimum(arr[adjacent], z + 1)

def main():
    arr = generate(20, 20, 5)
    smooth(arr)
    print(arr)

if __name__ == '__main__':
    main()

out.txt

$ python heightmap_smooth.py
[[1 1 2 2 2 2 1 0 1 2 1 0 1 0 1 2 2 1 1 1]
 [2 2 1 1 2 1 1 1 0 1 2 1 2 1 0 1 2 2 2 1]
 [2 1 0 1 1 0 1 1 1 2 1 2 1 2 1 2 2 2 1 0]
 [1 1 0 1 2 1 2 1 2 1 0 1 0 1 1 1 1 2 2 1]
 [1 2 1 1 2 2 2 1 2 2 1 0 1 1 1 2 2 3 2 2]
 [1 2 1 2 3 2 1 0 1 2 1 1 1 0 1 1 2 2 1 2]
 [0 1 2 1 2 1 2 1 1 1 2 1 0 0 0 1 2 1 0 1]
 [1 1 1 0 1 0 1 2 1 0 1 2 1 1 0 1 1 0 1 2]
 [1 0 1 1 1 0 1 2 1 1 2 1 0 1 1 2 1 0 1 2]
 [2 1 1 0 1 1 2 1 2 2 2 2 1 0 0 1 1 1 1 1]
 [2 1 2 1 1 1 1 2 2 2 1 2 1 1 1 1 0 1 0 1]
 [2 1 1 0 0 1 0 1 1 1 0 1 0 0 0 0 1 0 1 1]
 [1 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 2 1 1 0]
 [1 0 0 1 1 0 1 1 0 1 1 1 1 0 1 2 1 0 1 1]
 [0 1 0 1 2 1 1 0 1 2 1 1 2 1 2 1 2 1 1 2]
 [1 1 1 0 1 1 1 1 1 1 0 0 1 0 1 1 1 0 1 2]
 [1 0 1 1 2 1 1 1 0 1 1 1 1 0 1 0 1 1 1 2]
 [2 1 2 2 1 0 1 2 1 2 2 1 0 1 1 1 0 1 1 1]
 [2 2 2 2 1 1 1 2 1 1 1 1 0 1 0 1 1 0 1 1]
 [1 2 1 1 2 2 2 1 0 1 2 1 1 0 0 1 1 0 1 2]]