API reference

geoms_to_geo(geom, output_name; kwargs...) -> String

Convert any GeoInterface-compatible geometry (or a DataFrame with a geometry column) to a Gmsh .geo script. output_name should be given without the .geo extension.

The pipeline runs these steps before writing:

1. Reproject (GeometryOps.reproject) — if target_crs is set
2. Simplify (GeometryOps.simplify) — if simplify_alg is set
3. Segmentize (GeometryOps.segmentize) — if max_distance is set
4. Ingest (ingest) — convert to Gmsh-ready internal representation
5. Filter (filter_components) — drop degenerate rings
6. Rescale (rescale) — if bbox_size is set

Keyword arguments

• target_crs — destination CRS (e.g. "EPSG:32632"), :auto_utm to detect the UTM zone automatically from the geometry's centroid (default), or nothing to skip reprojection.
• select — when geom is a file path, a predicate row -> Bool passed to read_geodata. Ignored otherwise.
• simplify_alg — any GO.SimplifyAlg (e.g. MinEdgeLength(tol=5_000.0), AngleFilter(tol=20.0), or a composition via ∘). nothing skips simplification.
• max_distance — maximum edge length for segmentization (same units as the reprojected CRS). nothing skips.
• bbox_size — rescale so the largest bounding-box dimension equals this value, origin at (0, 0). nothing to skip.
• mesh_size — characteristic element length (default: 1.0).
• mesh_algorithm — Gmsh algorithm tag (default: nothing).
• split_components — write one .geo file per component (default: false).
• verbose — print progress (default: true).

Returns output_name.

geoms_to_msh(geom, output_name; kwargs...) -> String

Same as geoms_to_geo but generates a 2-D Gmsh mesh (.msh file).

Additional keyword arguments (beyond those of geoms_to_geo):

• order — element order: 1 = linear (default), 2 = quadratic.
• recombine — recombine triangles into quadrilaterals (default: false).

shapefile_to_geo(path, output_name; kwargs...) -> String

Backward-compatible wrapper: read a geospatial file and produce a .geo file. All keyword arguments are forwarded to geoms_to_geo.

shapefile_to_msh(path, output_name; kwargs...) -> String

Backward-compatible wrapper: read a geospatial file and produce a .msh file. All keyword arguments are forwarded to geoms_to_msh.

read_geodata(path; layer = nothing, select = nothing) -> DataFrame

Read any geospatial file supported by GeoDataFrames (Shapefile, GeoJSON, GeoPackage, GeoParquet, GeoArrow, FlatGeobuf, and any GDAL-supported format).

Returns a standard DataFrame with a :geometry column of GeoInterface-compatible geometries. CRS metadata is accessible via GeoInterface.crs(df).

Keyword arguments

• layer — layer index (Int, 0-based) or name (String) for multi-layer formats such as GeoPackage. Defaults to the first layer.
• select — a predicate row -> Bool to filter rows after reading.

list_components(path; layer = nothing) -> DataFrame

Read the geospatial file at path and print a summary table that includes all attribute columns plus per-geometry statistics (:n_pts, :area, :xmin, :xmax, :ymin, :ymax). Returns the augmented DataFrame.

Useful for exploring a file before deciding which features to select. For multi-layer formats, pass layer to inspect a specific layer.

read_shapefile(path; layer = nothing, select = nothing) -> DataFrame

Thin backward-compatible wrapper around read_geodata.

MinEdgeLength(; tol)

Simplification algorithm that removes any vertex whose distance to the previously kept vertex is less than tol.

Unlike RadialDistance (which measures to the last visited point), MinEdgeLength measures to the last retained point, giving a strict guarantee: no edge in the output is shorter than tol.

This is useful for Gmsh mesh quality — very short boundary edges can cause degenerate mesh elements or meshing failures.

Example

import GeometryOps as GO
simplified = GO.simplify(MinEdgeLength(tol = 5_000.0), polygon)

AngleFilter(; tol)

Simplification algorithm that removes any vertex whose interior angle (∠ a–b–c, in degrees) is below tol.

A zig-zag spike has a very acute interior angle at the tip: the two neighbouring segments nearly reverse direction, so the vectors b→a and b→c point almost the same way and the angle between them is close to 0°. Removing such vertices flattens spikes without touching smooth curves or genuine corners.

A single pass may leave new acute angles at formerly-adjacent points, so the algorithm iterates until no further points are removed.

Example

import GeometryOps as GO
simplified = GO.simplify(AngleFilter(tol = 20.0), polygon)

ComposedAlg(first, second)

Applies two simplification algorithms in sequence: first is run first, its output is passed to second.

Use the ∘ operator as syntactic sugar (note: a ∘ b applies b first):

alg = MinEdgeLength(tol = 5_000.0) ∘ AngleFilter(tol = 20.0)
# equivalent to: ComposedAlg(AngleFilter(tol = 20.0), MinEdgeLength(tol = 5_000.0))
simplified = GO.simplify(alg, polygon)

AdaptivityAlgorithm

Abstract base type for spatially-varying mesh-size strategies.

Pass an instance as the mesh_size argument of generate_mesh, geoms_to_msh, geoms_to_msh_3d, etc. Passing a plain Real uses a uniform characteristic length (the original behaviour).

Implementing a custom strategy

struct MyAdaptivity <: AdaptivityAlgorithm
    coarse :: Float64
    fine   :: Float64
end

GeoGmsh.nominal_size(alg::MyAdaptivity) = alg.coarse

function GeoGmsh.apply_adaptivity!(alg::MyAdaptivity, dem, curve_tags)
    # build Gmsh fields; return the background field tag (or nothing)
    f = gmsh.model.mesh.field.add("MathEval")
    gmsh.model.mesh.field.setString(f, "F", "$(alg.fine)")
    return f
end

nominal_size(alg::AdaptivityAlgorithm) -> Float64

Return the largest (coarsest) mesh size that alg may produce.

This value is used as:

• the characteristic length (lc) passed to Gmsh point entities, and
• the Mesh.CharacteristicLengthMax global bound.

Must be implemented for every concrete AdaptivityAlgorithm subtype.

apply_adaptivity!(alg::AdaptivityAlgorithm, dem, curve_tags) -> Union{Int,Nothing}

Build Gmsh background mesh-size fields for alg inside an active Gmsh session (i.e. after gmsh.initialize() and geometry synchronisation).

Returns the Gmsh field tag to register as the background mesh, or nothing if the strategy needs no field (e.g. because it degenerates to a uniform size).

Arguments

• dem — DEMRaster or nothing if no DEM is available (2D pipeline). Strategies that require a DEM should call error(...) when this is nothing.
• curve_tags — Vector{Int} of all curve entity tags in the current model. Strategies that need boundary distance should use these.

Must be implemented for every concrete AdaptivityAlgorithm subtype.

SlopeAdaptivity(; mesh_size_min, mesh_size_max, percentile = 0.95)

Slope-based spatially-varying mesh size: steep terrain gets small elements, flat terrain gets large elements.

The DEM gradient magnitude |∇h| is computed via finite differences on the raster grid. Values are normalised by the percentile-th percentile of the gradient distribution (default: 95th), so outlier cliffs do not compress the whole scale. The resulting mesh size at each DEM pixel is:

t    = clamp(|∇h| / slope_at_percentile, 0, 1)
size = mesh_size_min + (mesh_size_max - mesh_size_min) * (1 − t)

This produces mesh_size_max on perfectly flat terrain and mesh_size_min at slopes at or above the percentile threshold.

Requires a DEM — only usable with geoms_to_msh_3d, geoms_to_msh_3d_volume, or by passing a DEMRaster directly to generate_mesh.

Example

geoms_to_msh_3d("region.geojson", "dem.tif", "output";
    mesh_size  = SlopeAdaptivity(mesh_size_min = 100.0, mesh_size_max = 500.0),
    target_crs = "EPSG:32632",
)

BoundaryLayerAdaptivity(; mesh_size, mesh_size_min, width)

Boundary-layer mesh refinement: small elements right at polygon edges, smoothly coarsening to mesh_size over a distance of width (in the mesh CRS units, typically metres for a UTM projection).

Uses Gmsh's Distance + Threshold fields:

size(d) = mesh_size_min   for d ≤ 0
         ↗ (linear ramp)
        mesh_size          for d ≥ width

Works in both 2D and 3D pipelines (curve tags are available in both).

Example

geoms_to_msh("region.geojson", "output";
    mesh_size = BoundaryLayerAdaptivity(
        mesh_size     = 500.0,
        mesh_size_min =  50.0,
        width         = 2_000.0,
    ),
)

ComposedAdaptivity(alg1, alg2, ...)
ComposedAdaptivity(algorithms::Vector{AdaptivityAlgorithm})

Combine multiple AdaptivityAlgorithms using element-wise minimum (the finest size wins at every point in the domain). Implemented with Gmsh's Min field.

Example — slope refinement + boundary layer together

geoms_to_msh_3d("region.geojson", "dem.tif", "output";
    mesh_size = ComposedAdaptivity(
        SlopeAdaptivity(mesh_size_min = 100.0, mesh_size_max = 500.0),
        BoundaryLayerAdaptivity(mesh_size = 500.0, mesh_size_min = 50.0, width = 1_000.0),
    ),
    target_crs = "EPSG:32632",
)

rescale(geoms, L) -> Vector{Geometry2D}

Uniformly scale and translate geoms so that the largest dimension of the global bounding box equals L and the minimum corner is at the origin.

filter_components(geoms; min_points = 4) -> Vector{Geometry2D}

Remove geometrically degenerate components:

• Drop any Geometry2D whose exterior ring has fewer than min_points vertices.
• Strip any hole with fewer than min_points vertices (the exterior is kept).

The default min_points = 4 removes 3-point (triangular) rings that survive simplification and can cause Gmsh meshing failures.

ingest(geom) -> Vector{Geometry2D}

Convert any GeoInterface-compatible geometry to our internal representation.

Accepted inputs:

• Any GI.PolygonTrait geometry (single polygon)
• Any GI.MultiPolygonTrait geometry (split into one Geometry2D each)
• Any GI.FeatureTrait (geometry is extracted; properties are dropped)
• Any GI.FeatureCollectionTrait (all features ingested and concatenated)
• A DataFrame with a geometry column (as returned by read_geodata)

Ring orientation is normalised: exterior rings → CCW, hole rings → CW. Coordinates are always converted to Float64. Closing points (repeated first point) are removed automatically.

write_geo(geoms, name; mesh_size=1.0, mesh_algorithm=nothing,
          split_components=false)

Write geoms to a Gmsh .geo script. name should be given without the .geo extension — it is always appended internally.

When split_components = false (default), all geometries are written into a single file name.geo.

When split_components = true, one .geo file per Geometry2D is written into a directory named name/. Files inside are named 1.geo, 2.geo, … .

Keyword arguments

• mesh_size — characteristic element length assigned to every point.
• mesh_algorithm — optional integer passed as Mesh.Algorithm (e.g. 5 = Delaunay, 6 = Frontal-Delaunay, 8 = Frontal-Quad). Omitted from the file when nothing.
• split_components — write one file per geometry component (default false).

write_geo(geoms, name; kwargs...)

Geometry3D overload: writes boundary points with their terrain elevation (non-zero z).

When dem is provided, interior points are sampled on a regular grid at spacing max(mesh_size, DEM pixel size), filtered to lie inside each polygon, and embedded in the surface via Point{...} In Surface{...};. This encodes the full terrain curvature in the .geo file so that gmsh name.geo -2 produces a terrain-following surface mesh without any additional post-processing.

When dem is nothing (default), only boundary elevation is written.

Keyword arguments are identical to the 2D version, plus:

• dem — optional DEMRaster for interior elevation sampling.
• nodata_fill — elevation used for nodata / out-of-bounds cells (default 0.0).

generate_mesh(geoms, name; mesh_size=1.0, mesh_algorithm=nothing,
              order=1, recombine=false, split_components=false)

Build the geometry with the Gmsh API, generate a 2-D mesh, and write a .msh file. name should be given without the .msh extension.

When split_components = true, one .msh file per Geometry2D is written into a directory named name/.

Keyword arguments

• mesh_size — characteristic element length (sets both the min and max bounds passed to Gmsh).
• mesh_algorithm — Gmsh 2-D algorithm tag (e.g. 5 = Delaunay, 6 = Frontal-Delaunay, 8 = Frontal-Quad). Uses Gmsh's default when nothing.
• order — element order: 1 = linear (default), 2 = quadratic.
• recombine — recombine triangles into quadrilaterals (default false).
• split_components — write one file per geometry component (default false).

generate_mesh(geoms, dem, name; kwargs...) -> String

Geometry3D overload: generates a terrain-following surface mesh.

Meshes the 2D domain flat, then lifts every node's z-coordinate by sampling dem at its (x, y) position. The DEM and the 2D geometries must be in the same CRS.

Keyword arguments are identical to the 2D version.

download_cop30_tiles(bbox, dest_dir; verbose=true) -> Vector{String}

Download Copernicus GLO-30 DEM tiles covering bbox = (lon_min, lat_min, lon_max, lat_max) from the public AWS S3 bucket into dest_dir. Tiles already present on disk are skipped (cached).

Returns the paths of all successfully downloaded (or cached) tiles.

download_opentopo_dem(bbox, output_path;
                      demtype  = :gebco,
                      api_key  = ENV["OPENTOPO_API_KEY"],
                      verbose  = true) -> String

Download a DEM from the OpenTopography global DEM API and save it as a GeoTIFF to output_path. If the file already exists it is returned immediately (cached — safe to call on every run).

Getting an API key

Register for free at https://opentopography.org/developers and generate a key in your account dashboard. Store it in your environment:

export OPENTOPO_API_KEY="your_key_here"   # add to ~/.bashrc or ~/.zshrc

Keyword arguments

• demtype — DEM source (default :gebco). Supported symbols:

  Symbol	Dataset	Resolution
  :gebco	GEBCO topo-bathy (ocean + land)	~500 m
  :cop30	Copernicus Global DEM	30 m
  :cop90	Copernicus Global DEM	90 m
  :srtm	SRTM (land only)	90 m
  :srtm1	SRTM (land only)	30 m
  :alos	ALOS World 3D (land only)	30 m

• api_key — OpenTopography API key. Defaults to ENV["OPENTOPO_API_KEY"].

• verbose — print progress (default true).

Example

bbox    = (4.0, 60.0, 6.0, 62.0)   # Sognefjord, Norway
raw_tif = download_opentopo_dem(bbox, "sognefjord_gebco.tif")
dem_tif = prepare_dem([raw_tif], "EPSG:32632", "sognefjord_32632.tif")

prepare_dem(tile_paths, target_crs, output_path; resolution=30.0, verbose=true)
prepare_dem(bbox, target_crs, output_path; source=:cop30, dest_dir=..., ...)

Mosaic raster tiles and reproject to target_crs, writing a GeoTIFF to output_path. If output_path already exists it is returned immediately (no work done — safe to call on every run).

The first form is generic: tile_paths can be any list of GDAL-readable raster files (GeoTIFF, VRT, HGT, NetCDF, …).

The second form accepts a bbox = (lon_min, lat_min, lon_max, lat_max) and downloads tiles automatically using the given source before mosaicking.

Keyword arguments (both forms)

• resolution — output pixel size in target_crs units (default 30.0; metres for a projected UTM CRS, which matches the GLO-30 native grid).
• verbose — print progress (default true).

Additional keyword arguments (bbox form)

• source — tile provider symbol; currently :cop30 (Copernicus GLO-30, freely available from AWS S3).
• dest_dir — directory for downloaded tiles (default: directory of output_path, or the current working directory if output_path has no directory component).

DEMRaster

An in-memory Digital Elevation Model raster.

Fields:

• data — elevation matrix, shape (ncols, nrows) (ArchGDAL convention). data[col, row] where col=1 is west, row=1 is north.
• transform — GDAL 6-element geotransform [x_origin, x_pixel, x_rot, y_origin, y_rot, y_pixel]. For north-up rasters y_pixel < 0.
• crs_wkt — WKT CRS string, or nothing if unavailable.
• nodata — nodata sentinel value, or nothing if not set.

read_dem(path) -> DEMRaster

Read a Digital Elevation Model raster file (GeoTIFF, SRTM .hgt, NetCDF, …) using ArchGDAL. Only the first band is read.

Returns a DEMRaster with the elevation data, geotransform, CRS, and nodata value.

sample_elevation(pts, dem; nodata_fill = 0.0) -> Vector{Float64}

Bilinearly interpolate elevation values at pts (a Vector{NTuple{2,Float64}}) from dem. Points that fall outside the raster extent or coincide with nodata cells are filled with nodata_fill (default 0.0).

lift_to_3d(geom,  dem; nodata_fill = 0.0) -> Geometry3D
lift_to_3d(geoms, dem; nodata_fill = 0.0) -> Vector{Geometry3D}

Sample elevation from dem at every boundary point of geom (or each element of geoms) and return the corresponding Geometry3D object(s).

The DEM and the 2D geometries must be in the same CRS.

geoms_to_geo_3d(geom, dem_path, output_name; kwargs...) -> String

Terrain-aware variant of geoms_to_geo: runs the standard 2D pipeline, reads the DEM at dem_path, lifts boundary points to terrain elevation, and writes a .geo file with 3D point coordinates.

The DEM must be in the same CRS as the (reprojected) vector data.

All 2D keyword arguments are supported, except bbox_size which is silently ignored: rescaling the geometry after ingest would make its (x, y) coordinates inconsistent with the DEM's geotransform. Use mesh_size in the native CRS units (metres for a projected CRS) instead.

Additional argument:

• nodata_fill — elevation substituted for out-of-bounds or nodata cells (default 0.0).

geoms_to_msh_3d(geom, dem_path, output_name; kwargs...) -> String

Terrain-aware variant of geoms_to_msh: runs the standard 2D pipeline to generate a flat mesh, then lifts every mesh node's z-coordinate by sampling dem_path at its (x, y) position.

The DEM must be in the same CRS as the (reprojected) vector data.

All 2D keyword arguments are supported, except bbox_size which is silently ignored: rescaling the geometry after ingest would make its (x, y) coordinates inconsistent with the DEM's geotransform. Use mesh_size in the native CRS units (metres for a projected CRS) instead.

Additional argument:

• nodata_fill — elevation substituted for out-of-bounds or nodata cells (default 0.0).

geoms_to_msh_3d_volume(geom, dem_path, output_name; kwargs...) -> String

Generate a volumetric (tetrahedral) terrain mesh.

Runs the same pipeline as geoms_to_msh_3d to produce a terrain surface mesh, then extrudes it to flat bounding planes, splitting each triangular prism into 3 tetrahedra.

The extrusion direction is controlled by depth:

• depth = d — extrude downward by d below min(z_terrain).
• depth = (d_below, d_above) — extrude both directions; the terrain surface becomes an Interface physical group.
• depth = (0.0, h) — extrude upward only.

z_bottom / z_top pin the flat planes to absolute elevations.

The DEM must be in the same CRS as the (reprojected) vector data. bbox_size is silently ignored for the same reason as in geoms_to_msh_3d.

See generate_mesh_volume for the full physical group reference.

generate_mesh_volume(geoms, dem, name; depth=nothing, z_bottom=nothing,
                     z_top=nothing, mesh_size=1.0, mesh_algorithm=nothing,
                     split_components=false, verbose=false)

Generate a volumetric tetrahedral mesh by extruding the terrain surface to flat bounding planes.

The extrusion direction is controlled by depth:

• depth = d (plain Real) — extrude downward by d below min(z_terrain).
• depth = (d_below, d_above) (2-tuple) — extrude both directions simultaneously. The terrain surface becomes an interior Interface physical group shared by Volume_Below and Volume_Above.
• depth = (0.0, h) — extrude upward only.

z_bottom and z_top pin the flat planes to absolute elevations and override the corresponding component of depth.

Physical groups written to the .msh file:

Arguments

• geoms — Vector{Geometry3D} (only the 2D base geometry is used).
• dem — DEMRaster for elevation sampling.
• name — output path without the .msh extension.

Keyword arguments

• depth — see above.
• z_bottom — absolute bottom elevation; overrides the below component.
• z_top — absolute top elevation; overrides the above component.
• mesh_size — characteristic element length (default 1.0).
• mesh_algorithm — Gmsh 2D algorithm tag.
• split_components — write one .msh per geometry component (default false).
• verbose — print progress (default false).

Geometry2D

A single polygon: an exterior ring plus zero or more hole rings. All rings are stored as Contour with closed = true. Exterior rings are oriented CCW (positive signed area); hole rings are oriented CW (negative signed area).

Geometry3D

A terrain-aware 2+1D polygon: a Geometry2D base (2D boundary topology) combined with per-point elevation sampled from a DEM raster.

z_exterior holds one elevation value per point in base.exterior. z_holes[k] holds one elevation value per point in base.holes[k].

Geometry3D objects are always derived from a Geometry2D + a DEM raster via lift_to_3d. They are never constructed from native 3D vector data (which does not exist in standard geospatial formats).

Contour

A closed or open sequence of 2D points. Closed contours represent polygon rings; open contours represent polylines.

Points are stored as NTuple{2,Float64} (x, y).

Number of points in a contour.

Total number of points across exterior + holes.

Number of edges in a contour.