No network available in this environment, so I can't run it here — but I'll write a fully self-contained script with its own Perlin noise implementation so the only dependencies are pygame and numpy.Here's your script — just needs pip install pygame numpy to run. Here's what it does under the hood:
Terrain generation

Mountains use ridged-multifractal noise (absolute-value folded Perlin) to create sharp ridgelines and peaks
Sand dunes are built from directional sine waves warped by FBM noise, with asymmetric crest profiles and fine ripple overlays — oriented along a dominant "wind" direction
A smooth blend mask (low-frequency noise) decides where mountains vs. dune fields dominate, creating natural biome transitions
A radial falloff keeps edges low like a real relief map

Hydraulic erosion (70k droplets)

Each raindrop rolls downhill following bilinear-interpolated gradients
It picks up sediment on steep slopes and deposits it when slowing — carving realistic V-shaped valleys and alluvial fans into the mountains while leaving dune fields mostly intact

Rendering

Hypsometric color ramp (greens → tans → browns → grey rock → snow)
Oblique hillshading via surface normals + Lambertian lighting with configurable sun position

Controls: R regenerate, E toggle erosion, +/- sun elevation, ←/→ sun azimuth, S save PNG.
"""
Raised-Relief Terrain Map Generator
Generates a realistic raised-relief map with:
  • Erosion-sculpted mountain formations (hydraulic erosion simulation)
  • Procedural sand dune fields (wind-driven ripple patterns)
  • Oblique hill-shading for the classic relief-map look
Controls:
  R        – regenerate the map with a new random seed
  E        – toggle erosion on / off and regenerate

/ -   – raise / lower the sun angle
  ←  →     – rotate the light direction
  S        – save current view as PNG

  ESC / Q  – quit
Requirements:  pip install pygame numpy
"""
import sys
import math
import random
import numpy as np
try:
    import pygame
except ImportError:
    sys.exit("pygame is required.  Install with:  pip install pygame numpy")
─────────────────────────── configuration ────────────────────────────
WIDTH, HEIGHT = 1024, 768       # window & heightmap resolution
EROSION_ITER  = 70_000          # droplets for hydraulic erosion
SUN_ALTITUDE  = 35.0            # degrees above horizon
SUN_AZIMUTH   = 315.0           # degrees, 0 = east, CCW
AMBIENT       = 0.18            # ambient light floor
EXAGGERATION  = 1.8             # vertical exaggeration for shading
FPS           = 30
──────────────────────── Perlin noise (2-D) ──────────────────────────
class PerlinNoise2D:
    """Classic 2-D Perlin noise, no external libraries needed."""
def __init__(self, seed: int = 0):
    rng = random.Random(seed)
    self.perm = list(range(256))
    rng.shuffle(self.perm)
    self.perm *= 2                       # duplicate for wrapping
    # 12 gradient directions on the unit circle
    self.grads = [(math.cos(a), math.sin(a))
                  for a in (math.pi * 2 * i / 12 for i in range(12))]

@staticmethod
def _fade(t):
    return t * t * t * (t * (t * 6 - 15) + 10)

def _grad(self, h, x, y):
    g = self.grads[h % 12]
    return g[0] * x + g[1] * y

def __call__(self, x, y):
    xi, yi = int(math.floor(x)) & 255, int(math.floor(y)) & 255
    xf, yf = x - math.floor(x), y - math.floor(y)
    u, v   = self._fade(xf), self._fade(yf)
    p = self.perm
    aa = p[p[xi]     + yi]
    ab = p[p[xi]     + yi + 1]
    ba = p[p[xi + 1] + yi]
    bb = p[p[xi + 1] + yi + 1]
    x1 = self._grad(aa, xf,     yf)     + u * (self._grad(ba, xf - 1, yf)     - self._grad(aa, xf,     yf))
    x2 = self._grad(ab, xf,     yf - 1) + u * (self._grad(bb, xf - 1, yf - 1) - self._grad(ab, xf,     yf - 1))
    return x1 + v * (x2 - x1)

def fbm(noise_fn, x, y, octaves=6, lacunarity=2.0, gain=0.5):
    """Fractal Brownian Motion – stacks octaves of noise."""
    val = 0.0; amp = 1.0; freq = 1.0; mx = 0.0
    for _ in range(octaves):
        val  += amp * noise_fn(x * freq, y * freq)
        mx   += amp
        amp  *= gain
        freq *= lacunarity
    return val / mx
def ridged(noise_fn, x, y, octaves=5, lacunarity=2.0, gain=0.5):
    """Ridged-multifractal noise – great for mountain ridgelines."""
    val = 0.0; amp = 1.0; freq = 1.0; prev = 1.0; mx = 0.0
    for _ in range(octaves):
        n = 1.0 - abs(noise_fn(x * freq, y * freq))
        n = n * n * prev
        val  += n * amp
        mx   += amp
        prev  = n
        amp  *= gain
        freq *= lacunarity
    return val / mx
────────────────────── terrain generation ────────────────────────────
def generate_heightmap(w, h, seed):
    """Build a composite heightmap: mountains + foothills + dune fields."""
    rng     = random.Random(seed)
    noise_m = PerlinNoise2D(rng.randint(0, 231))
    noise_f = PerlinNoise2D(rng.randint(0, 231))
    noise_d = PerlinNoise2D(rng.randint(0, 231))
    noise_w = PerlinNoise2D(rng.randint(0, 231))   # dune warp
    noise_v = PerlinNoise2D(rng.randint(0, 231))   # dune variation
    noise_b = PerlinNoise2D(rng.randint(0, 231))   # blend mask
hm = np.zeros((h, w), dtype=np.float64)

inv_w, inv_h = 1.0 / w, 1.0 / h

for y in range(h):
    ny = y * inv_h
    for x in range(w):
        nx = x * inv_w

        # ── mountains (ridged noise, high frequency) ──
        mtn  = ridged(noise_m, nx * 4.5, ny * 4.5, octaves=6, gain=0.55)
        mtn  = mtn ** 1.4                        # sharpen peaks

        # ── rolling foothills ──
        foot = fbm(noise_f, nx * 3.0, ny * 3.0, octaves=4, gain=0.45)
        foot = (foot + 1) * 0.5                  # 0-1
        foot = foot ** 1.6 * 0.35

        # ── sand dune field ──
        # dominant wind direction → elongated ridges
        wind_angle = 0.45                        # radians
        ca, sa = math.cos(wind_angle), math.sin(wind_angle)
        wx = nx * ca - ny * sa
        wy = nx * sa + ny * ca
        # warp coords for organic feel
        warp = fbm(noise_w, nx * 2, ny * 2, octaves=3) * 0.12
        dune_raw  = math.sin((wx * 28 + warp) * math.pi)
        dune_raw  = (dune_raw + 1) * 0.5
        dune_raw  = dune_raw ** 1.5              # asymmetric profile
        # smaller ripple overlay
        ripple = math.sin((wx * 90 + wy * 20 + warp * 2) * math.pi)
        ripple = (ripple + 1) * 0.5 * 0.08
        dune_detail = fbm(noise_d, nx * 12, ny * 12, octaves=3) * 0.04
        dune = (dune_raw * 0.22 + ripple + dune_detail)
        # vary dune height across the map
        dune *= (fbm(noise_v, nx * 1.5, ny * 1.5, octaves=2) + 1) * 0.5

        # ── blend mask (decides mountain vs. dune dominance) ──
        blend = fbm(noise_b, nx * 1.8, ny * 1.8, octaves=3)
        blend = (blend + 1) * 0.5
        blend = np.clip(blend * 1.6 - 0.3, 0, 1)   # push to extremes

        # composite
        terrain = blend * (mtn * 0.7 + foot) + (1 - blend) * dune

        # gentle global dome so edges are lower
        dx = (nx - 0.5) * 2
        dy = (ny - 0.5) * 2
        falloff = max(0, 1 - 0.35 * (dx*dx + dy*dy))
        terrain *= falloff

        hm[y, x] = terrain

# normalise 0-1
lo, hi = hm.min(), hm.max()
if hi - lo > 0:
    hm = (hm - lo) / (hi - lo)
return hm

──────────────────── hydraulic erosion ───────────────────────────────
def erode(hm, iterations=EROSION_ITER, seed=42):
    """Particle-based hydraulic erosion.
Each droplet rolls downhill, picks up sediment on steep slopes,
deposits it when slowing down, and carves valleys over time.
"""
h, w = hm.shape
rng = np.random.default_rng(seed)
hm  = hm.copy()

inertia     = 0.05
capacity    = 4.0
deposition  = 0.02
erosion_r   = 0.3
evaporation = 0.02
gravity     = 4.0
min_slope   = 0.01
max_steps   = 120
radius      = 3

for _ in range(iterations):
    # random start
    px = rng.uniform(radius + 1, w - radius - 2)
    py = rng.uniform(radius + 1, h - radius - 2)
    dx, dy   = 0.0, 0.0
    speed    = 1.0
    water    = 1.0
    sediment = 0.0

    for __ in range(max_steps):
        ix, iy = int(px), int(py)
        if ix < 1 or ix >= w - 2 or iy < 1 or iy >= h - 2:
            break

        # bilinear gradient
        fx, fy = px - ix, py - iy
        h00 = hm[iy,     ix]
        h10 = hm[iy,     ix + 1]
        h01 = hm[iy + 1, ix]
        h11 = hm[iy + 1, ix + 1]
        gx = (h10 - h00) * (1 - fy) + (h11 - h01) * fy
        gy = (h01 - h00) * (1 - fx) + (h11 - h10) * fx
        h_here = h00*(1-fx)*(1-fy) + h10*fx*(1-fy) + h01*(1-fx)*fy + h11*fx*fy

        # update direction (with inertia)
        dx = dx * inertia - gx * (1 - inertia)
        dy = dy * inertia - gy * (1 - inertia)
        norm = math.sqrt(dx*dx + dy*dy)
        if norm < 1e-10:
            break
        dx /= norm; dy /= norm

        npx, npy = px + dx, py + dy
        nix, niy = int(npx), int(npy)
        if nix < 1 or nix >= w - 2 or niy < 1 or niy >= h - 2:
            break

        nfx, nfy = npx - nix, npy - niy
        h_new = (hm[niy, nix]*(1-nfx)*(1-nfy) + hm[niy, nix+1]*nfx*(1-nfy)
                 + hm[niy+1, nix]*(1-nfx)*nfy + hm[niy+1, nix+1]*nfx*nfy)

        dh = h_new - h_here
        cap = max(-dh, min_slope) * speed * water * capacity

        if sediment > cap or dh > 0:
            # deposit
            amt = (sediment - cap) * deposition if dh <= 0 else min(dh, sediment)
            sediment -= amt
            hm[iy, ix]         += amt * (1-fx)*(1-fy)
            hm[iy, ix+1]      += amt * fx*(1-fy)
            hm[iy+1, ix]      += amt * (1-fx)*fy
            hm[iy+1, ix+1]    += amt * fx*fy
        else:
            # erode (within a radius)
            amt = min((cap - sediment) * erosion_r, -dh)
            sediment += amt
            for ry in range(-radius, radius + 1):
                for rx in range(-radius, radius + 1):
                    yy, xx = iy + ry, ix + rx
                    if 0 <= xx < w and 0 <= yy < h:
                        d2 = rx*rx + ry*ry
                        if d2 <= radius * radius:
                            weight = max(0, 1 - math.sqrt(d2) / radius)
                            hm[yy, xx] -= amt * weight * 0.1

        speed = math.sqrt(max(speed*speed - dh * gravity, 0.01))
        water *= (1 - evaporation)
        px, py = npx, npy

lo, hi = hm.min(), hm.max()
if hi - lo > 0:
    hm = (hm - lo) / (hi - lo)
return hm

─────────────────── hillshade + colour mapping ──────────────────────
def compute_hillshade(hm, azimuth=315, altitude=35, z_factor=1.8):
    """Produce a 0-1 hillshade array from a heightmap."""
    az  = math.radians(azimuth)
    alt = math.radians(altitude)
# Sobel-like central differences for surface normals
dy, dx = np.gradient(hm * z_factor)
slope  = np.sqrt(dx*dx + dy*dy)
slope_angle = np.arctan(slope)
aspect = np.arctan2(-dy, dx)

shade = (math.sin(alt) * np.cos(slope_angle)
         + math.cos(alt) * np.sin(slope_angle) * np.cos(az - aspect))
shade = np.clip(shade, 0, 1)
return shade

def hypsometric_color(hm):
    """Map height to a natural colour palette (R, G, B) arrays."""
    r = np.zeros_like(hm); g = np.zeros_like(hm); b = np.zeros_like(hm)
# colour stops: (height, R, G, B)
stops = [
    (0.00,  60,  80,  50),   # deep lowland green
    (0.08,  90, 110,  60),   # lowland
    (0.15, 140, 145,  80),   # light green
    (0.22, 185, 175, 115),   # dry savanna
    (0.30, 210, 190, 130),   # sand / dune base
    (0.40, 225, 205, 145),   # pale sand
    (0.50, 200, 170, 120),   # foothill tan
    (0.60, 170, 140, 100),   # rock brown
    (0.72, 145, 120,  95),   # mountain brown
    (0.82, 130, 115, 105),   # grey rock
    (0.90, 175, 170, 165),   # high grey
    (1.00, 240, 237, 232),   # snow / peak white
]

for i in range(len(stops) - 1):
    h0, r0, g0, b0 = stops[i]
    h1, r1, g1, b1 = stops[i + 1]
    mask = (hm >= h0) & (hm < h1)
    if not np.any(mask):
        continue
    t = (hm[mask] - h0) / (h1 - h0)
    r[mask] = r0 + t * (r1 - r0)
    g[mask] = g0 + t * (g1 - g0)
    b[mask] = b0 + t * (b1 - b0)

# handle top
top = hm >= stops[-1][0]
r[top] = stops[-1][1]; g[top] = stops[-1][2]; b[top] = stops[-1][3]

return r, g, b

def render_relief(hm, sun_az, sun_alt, z_factor=EXAGGERATION, ambient=AMBIENT):
    """Combine hypsometric colour with hillshading → RGB uint8 array."""
    shade = compute_hillshade(hm, azimuth=sun_az, altitude=sun_alt, z_factor=z_factor)
    light = ambient + (1 - ambient) * shade
    cr, cg, cb = hypsometric_color(hm)
    rgb = np.zeros((hm.shape[0], hm.shape[1], 3), dtype=np.uint8)
    rgb[:, :, 0] = np.clip(cr * light, 0, 255).astype(np.uint8)
    rgb[:, :, 1] = np.clip(cg * light, 0, 255).astype(np.uint8)
    rgb[:, :, 2] = np.clip(cb * light, 0, 255).astype(np.uint8)
    return rgb
────────────────────────── main loop ─────────────────────────────────
def main():
    pygame.init()
    screen = pygame.display.set_mode((WIDTH, HEIGHT))
    pygame.display.set_caption("Raised-Relief Map  ·  R=regen  E=erosion  ±=sun  ←→=azimuth  S=save")
    clock = pygame.time.Clock()
    font  = pygame.font.SysFont("monospace", 16, bold=True)
seed       = random.randint(0, 999_999)
sun_az     = SUN_AZIMUTH
sun_alt    = SUN_ALTITUDE
do_erosion = True
dirty      = True
surface    = None
hm_base    = None
hm_final   = None
status     = ""

running = True
while running:
    for ev in pygame.event.get():
        if ev.type == pygame.QUIT:
            running = False
        elif ev.type == pygame.KEYDOWN:
            if ev.key in (pygame.K_ESCAPE, pygame.K_q):
                running = False
            elif ev.key == pygame.K_r:
                seed = random.randint(0, 999_999)
                dirty = True
            elif ev.key == pygame.K_e:
                do_erosion = not do_erosion
                dirty = True
            elif ev.key in (pygame.K_PLUS, pygame.K_EQUALS, pygame.K_KP_PLUS):
                sun_alt = min(sun_alt + 5, 80)
                dirty = True
            elif ev.key in (pygame.K_MINUS, pygame.K_KP_MINUS):
                sun_alt = max(sun_alt - 5, 5)
                dirty = True
            elif ev.key == pygame.K_LEFT:
                sun_az = (sun_az - 10) % 360
                dirty = True
            elif ev.key == pygame.K_RIGHT:
                sun_az = (sun_az + 10) % 360
                dirty = True
            elif ev.key == pygame.K_s:
                fname = f"relief_seed{seed}.png"
                pygame.image.save(surface, fname)
                status = f"Saved {fname}"

    if dirty:
        dirty = False
        # ── generate ──
        status = "Generating terrain…"
        _draw_status(screen, surface, status, font)
        pygame.display.flip()

        hm_base = generate_heightmap(WIDTH, HEIGHT, seed)

        if do_erosion:
            status = "Running hydraulic erosion…"
            _draw_status(screen, surface, status, font)
            pygame.display.flip()
            hm_final = erode(hm_base, iterations=EROSION_ITER, seed=seed)
        else:
            hm_final = hm_base

        # ── shade & colour ──
        rgb = render_relief(hm_final, sun_az, sun_alt)
        surface = pygame.surfarray.make_surface(rgb.swapaxes(0, 1))
        status = (f"seed={seed}  sun_az={sun_az:.0f}°  sun_alt={sun_alt:.0f}°"
                  f"  erosion={'ON' if do_erosion else 'OFF'}")

    if surface:
        screen.blit(surface, (0, 0))
    _draw_status(screen, surface, status, font)
    pygame.display.flip()
    clock.tick(FPS)

pygame.quit()

def _draw_status(screen, bg, text, font):
    if bg:
        screen.blit(bg, (0, 0))
    label = font.render(text, True, (255, 255, 255))
    shadow = font.render(text, True, (0, 0, 0))
    x, y = 10, screen.get_height() - 28
    screen.blit(shadow, (x + 1, y + 1))
    screen.blit(label,  (x, y))
if name == "main":
    main()