Mont Blanc — 3D terrain mesh

This example produces two terrain-following meshes for the Mont Blanc massif using a user-defined bounding box and a Copernicus GLO-30 Digital Elevation Model:

  • Surface mesh (geoms_to_msh_3d): a terrain-following surface with z-coordinates sampled from the DEM.
  • Volume mesh (geoms_to_msh_3d_volume): the same surface extruded downward by depth metres to form a solid tetrahedral block.

Features highlighted:

  • prepare_dem: one-call DEM download, mosaic, and reproject
  • geoms_to_msh_3d and geoms_to_msh_3d_volume
  • UTM CRS so that mesh_size and depth are in metres

Aspect ratio

The Mont Blanc massif spans ~37 km × ~27 km horizontally but only ~4,800 m vertically. At true scale the mesh looks flat — apply vertical exaggeration (e.g. 5×) in your visualiser (ParaView, QGIS 3D, Visit).

Mont Blanc (3D terrain)
Mont Blanc mesh
using GeoGmsh

data_dir = joinpath(@__DIR__, "..", "data")
mkpath(data_dir)

Bounding box

We define the domain of interest as a simple rectangular polygon in EPSG:4326. A small padding around the target area avoids edge effects in the mesh.

Original bounding box: 6.785431, 45.785243, 7.044296, 45.956878

lon_min, lat_min = 6.73, 45.75
lon_max, lat_max = 7.10, 45.99
bbox = (lon_min, lat_min, lon_max, lat_max)

bbox_path = joinpath(data_dir, "montblanc_bbox.geojson")
open(bbox_path, "w") do io
  write(io, """
{
  "type": "FeatureCollection",
  "features": [{
    "type": "Feature",
    "geometry": {
      "type": "Polygon",
      "coordinates": [[[$(lon_min),$(lat_min)],[$(lon_max),$(lat_min)],
                        [$(lon_max),$(lat_max)],[$(lon_min),$(lat_max)],
                        [$(lon_min),$(lat_min)]]]
    },
    "properties": {}
  }]
}
""")
end

Prepare DEM

prepare_dem downloads the required Copernicus GLO-30 tiles (skipping any already on disk), stitches them into a seamless mosaic with gdalbuildvrt, and reprojects to UTM zone 32N (EPSG:32632) at 30 m resolution with gdalwarp — all in one call. The domain sits in the N45 latitude band and spans longitude columns E006 and E007, so two tiles are fetched.

dem_tif_utm = prepare_dem(
  bbox, "EPSG:32632",
  joinpath(data_dir, "montblanc_dem_32632.tif");
  source   = :cop30,
  dest_dir = data_dir,
)

3D terrain meshes

geoms_to_msh_3d runs the standard 2D pipeline (reproject → simplify → ingest) and then lifts every mesh node's z-coordinate by bilinearly interpolating the DEM at its (x, y) position. mesh_size = 500.0 gives ~500 m characteristic element length.

geoms_to_msh_3d_volume does the same but also extrudes the surface downward by depth metres to produce a solid tetrahedral mesh. A depth of 1,000 m gives a thin but clearly visible pedestal relative to the ~4,800 m of vertical relief in the massif.

output = joinpath(data_dir, "montblanc")

geoms_to_msh_3d(
  bbox_path, dem_tif_utm, output;
  target_crs   = "EPSG:32632",
  simplify_alg = MinEdgeLength(tol = 500.0),
  mesh_size    = 500.0,
  nodata_fill  = 0.0,
  verbose      = true,
)

geoms_to_msh_3d_volume(
  bbox_path, dem_tif_utm, output * "_volume";
  target_crs   = "EPSG:32632",
  simplify_alg = MinEdgeLength(tol = 500.0),
  mesh_size    = 500.0,
  depth        = 1_000.0,
  nodata_fill  = 0.0,
  verbose      = true,
)

This page was generated using Literate.jl.