DHI.Spatial – Internal Developer Guide

This module provides a lightweight spatial core used across services for geometry, features, meshes, and spatially-aware data sources. It's deliberately minimal: easy to embed, no heavy dependencies, and suitable as a building block for MIKE providers and Web APIs.

What's in the box

• Geometry model: Point, LineString, Polygon, Multi* types, GeometryCollection, WKT parsing/serialization, Position, BoundingBox, and CRS helpers.
• Feature model: Feature, FeatureCollection, Attribute, Association.
• Spatial data model: Node, Element, Mesh, SpatialDefinition* (EqD2 and NonEqD2 implemented), plus contracts for data sources (IDataSource<TItemId>).
• Utilities: EPSG:3857 ("Google / Web Mercator") conversions, tolerant string ops, numeric parsing.

  Many interfaces are designed so you can plug in MIKE/DFS-backed implementations without taking a dependency here.

---

Geometry layer

Core types

• Position – mutable coordinate class with X, Y, and optional Z.
  • ToWKT() / FromWKT(string) (supports xyz)
  • ProjectTo(Position origin, double dx, double dy) -> local, grid-index-like coordinates
  • ToGoogle() (lon/lat -> EPSG:3857 meters)
• BoundingBox – axis-aligned rectangle with helpers:
  • Parse("xmin, ymin, xmax, ymax")
  • Intersects(BoundingBox), Contains(Position)
  • GoogleToLonLat(), LonLatToGoogle()
  • GetPositions() returns LL, LR, UR, UL corners
• CoordinateReferenceSystem – holder for a CRS Type and arbitrary Properties.

Geometry objects

All geometries implement: * IGeometry -> Type (e.g., "Polygon"), Coordinates (dynamic; strongly-typed collections in practice), CRS, ToWKT(). * Concrete types: * Point -> Coordinates: Position * LineString -> List<Position> * Polygon -> List<List<Position>> (outer ring + zero or more inner rings) * MultiPoint, MultiLineString, MultiPolygon -> nested lists of Position * GeometryCollection -> List<IGeometry>

WKT parsing/formatting

• Geometry.FromWKT(string) supports:
  • POINT [Z], LINESTRING [Z], POLYGON [Z], MULTIPOINT, MULTILINESTRING, MULTIPOLYGON
• GeometryCollection has ToWKT() (no FromWKT).
• Throws NotSupportedException for types not listed above; throws with friendly messages on malformed WKT.

Web Mercator helpers (EPSG:3857)

DHI.Spatial.ExtensionMethods: * LonToGoogleX(), LatToGoogleY(), GoogleXToLon(), GoogleYToLat() * IsGoogle("EPSG:3857") (or "EPSG:900913") * SameAs(string, string) (case-insensitive) * Numeric helpers: ToDouble(string), ToFloat(double) Tip: prefer storing geometry in a well-defined CRS (WKT in CRS) and convert on IO boundaries. The helpers are simple and assume WGS84 input for lon/lat.

---

Feature model

• Feature = Geometry + bags for:
  • AttributeValues: Dictionary<string, object>
  • Associations: IList<IAssociation> (e.g., link to a domain entity)
• FeatureCollection = Features + schema-level Attributes (IAttribute: name, display name, type, length, default).
• Association<T> carries the associated CLR type and id. Serialization hints: FeatureCollection and Feature expose ShouldSerialize*() to omit empty arrays when JSON-serializing.

---

Mesh & spatial definitions

The mesh model abstracts discretized spatial domains: * INode<TId> -> Id, Position * IElement<TId, TNodeId> -> Id, IList<INode<TNodeId>> * IMesh<TElementId, TNodeId> -> Elements, Nodes (deduplicated via node equality by Id) Concrete implementations: * Node<TId> and Element<TId,TNodeId> (value semantics by Id). * Mesh<TElementId,TNodeId> (holds elements; Nodes computed from elements).

Spatial definitions (ISpatialDefinition)

A spatial definition describes a domain compactly and can create a mesh: * Type, DimensionCount * ElementCount() / ElementCount(int dimension) * ToMesh() -> IMesh<int,int> Implemented: * SpatialDefinitionEqD2 (2D equidistant grid): * ctor: (x0, y0, dx, dy, nx, ny) * ToMesh() creates a structured quadrilateral mesh with: * Element IDs: 1..(nx*ny) left->right, bottom->top * Node IDs: 1..((nx+1)*(ny+1)) left->right, bottom->top Upper nodes are offset by nx+1 from the lower nodes. * SpatialDefinitionNonEqD2 (2D non-equidistant grid): * ctor: (xCoordinates, yCoordinates) — node coordinate arrays, each strictly increasing, length nx+1 and ny+1 respectively. * ToMesh() creates the same quadrilateral structure as EqD2 but derives node coordinates directly from the provided arrays instead of a uniform step. * ElementCount() returns nx * ny; ElementCount(int dimension) returns nx (1) or ny (2).

Not yet implemented (throw NotImplementedException): SpatialDefinitionCurveLinearD2.

Element & node numbering (EqD2)

For cell (ix, iy) (0-based), lower-left node id = iy*(nx+1) + ix + 1. The four corners are:

LL = base
LR = base + 1
UR = base + (nx + 1) + 1
UL = base + (nx + 1)

ToMesh() follows exactly this scheme. Example (nx=2, ny=2): Nodes (ids): Element 1 nodes: 1–2–5–4 Element 4 nodes: 5–6–9–8

---

Data sources

IDataSource<TItemId> defines a contract for time-varying spatial values:

public interface IDataSource<TItemId>
{
    IList<DateTime> DateTimes { get; }
    DateTime? FirstDateTime { get; }
    DateTime? LastDateTime { get; }
    IEnumerable<Item<TItemId>> Items { get; }
    IEnumerable<string> ItemIds { get; }
    bool ContainsDateTime(DateTime dateTime);
    bool ContainsItem(string itemId);
    IMesh<int,int> Mesh { get; }              // typically derived from SpatialDefinition
    ISpatialDefinition SpatialDefinition { get; }
    double[] Get(TItemId itemId, DateTime dateTime, IGeometry geometry);
    IDictionary<TItemId, IDictionary<DateTime, double>> Get(
        IDictionary<TItemId, IEnumerable<DateTime>> ids,
        IGeometry geometry);
    double[] GetFirstAfter(TItemId itemId, DateTime dateTime, IGeometry geometry);
    double[] GetLastBefore(TItemId itemId, DateTime dateTime, IGeometry geometry);
}

BaseDataSource<TItemId> provides a concrete base with several members already implemented: * Mesh — delegates to SpatialDefinition.ToMesh(). * FirstDateTime / LastDateTime — derived from DateTimes. * ItemIds — derived from Items. * ContainsDateTime / ContainsItem — simple lookups. * GetFirstAfter / GetLastBefore — scan DateTimes and delegate to Get.

You must implement DateTimes, Items, SpatialDefinition, and both Get overloads. Everything else is ready to use.

Recommended semantics (for implementers)

Because the interface leaves some freedom, keep these conventions consistent within our codebase: * Items: scalar variables (e.g., "WaterLevel", "U-velocity"). Use Item<TId> to carry Quantity, Unit, PhysicalDimension. * Get(itemId, dateTime, geometry): * Point -> return a single-element array: sampled value at the point (interpolated if relevant). * Polygon / LineString / Multi* -> return an array of one value per intersecting element, ordered by element id. Aggregation (average/min/max) should be part of the calling service, unless your data source is explicitly "aggregated". * GetFirstAfter / GetLastBefore: * Same shape as Get, but time snapped to first > or last < the given dateTime. * Bulk Get(ids, geometry): * Return a dictionary ItemId -> (DateTime -> AggregateValue). For non-point geometries, use average over intersecting elements unless a different rule is documented in the source.

If your data product has a well-defined aggregation (e.g., area-weighted average per polygon), document it on the repository/service that wraps the data source.

Minimal concrete example (EqD2 grid)

public sealed class GridDataSource : BaseDataSource<string>
{
    private readonly SpatialDefinitionEqD2 _grid;
    private readonly List<DateTime> _times;
    private readonly Dictionary<(string item,int t,int elem), double> _data;
    public GridDataSource(SpatialDefinitionEqD2 grid,
                          IEnumerable<DateTime> dateTimes,
                          IEnumerable<Item<string>> items)
    {
        _grid = grid;
        _times = dateTimes.OrderBy(t => t).ToList();
        Items = items.ToArray();
        ItemIds = Items.Select(i => i.Id);
        _data = new();
    }
    public override ISpatialDefinition SpatialDefinition => _grid;
    public override IEnumerable<Item<string>> Items { get; }
    public override IEnumerable<string> ItemIds { get; }
    public override IList<DateTime> DateTimes => _times;
    public override DateTime? FirstDateTime => _times.FirstOrDefault();
    public override DateTime? LastDateTime  => _times.LastOrDefault();
    public override bool ContainsDateTime(DateTime dt) => _times.BinarySearch(dt) >= 0;
    public override bool ContainsItem(string id) => ItemIds.Contains(id);
    public void Set(string item, DateTime dt, int elementId, double value)
    {
        var ti = _times.IndexOf(dt);
        _data[(item, ti, elementId)] = value;
    }
    public override double[] Get(string item, DateTime dateTime, IGeometry geometry)
    {
        var ti = _times.IndexOf(dateTime);
        if (ti < 0) throw new ArgumentOutOfRangeException(nameof(dateTime));
        var elementIds = IntersectingElementIds(geometry);
        return elementIds.Select(eid => _data.TryGetValue((item, ti, eid), out var v) ? v : double.NaN)
                         .ToArray();
    }
    public override IDictionary<string, IDictionary<DateTime, double>> Get(
        IDictionary<string, IEnumerable<DateTime>> ids, IGeometry geometry)
    {
        var elem = IntersectingElementIds(geometry).ToArray();
        var result = new Dictionary<string, IDictionary<DateTime, double>>();
        foreach (var (item, dts) in ids)
        {
            result[item] = dts.ToDictionary(
                dt => dt,
                dt => elem.Average(eid =>
                {
                    var ti = _times.IndexOf(dt);
                    return _data.TryGetValue((item, ti, eid), out var v) ? v : double.NaN;
                }));
        }
        return result;
    }
    public override double[] GetFirstAfter(string item, DateTime dt, IGeometry g)
        => Get(item, _times.First(t => t > dt), g);
    public override double[] GetLastBefore(string item, DateTime dt, IGeometry g)
        => Get(item, _times.Last(t => t < dt), g);
    private IEnumerable<int> IntersectingElementIds(IGeometry g)
    {
        // Simplest: map point to nearest cell by (x0,dx,y0,dy).
        // For polygons, compute bbox and include all elements whose cell-centroid falls inside.
        // Implementations may use robust polygon/mesh intersection.
        var mesh = Mesh;
        if (g is Point pt) return new[] { FindElementIdForPoint(pt.Coordinates) };
        if (g is Polygon poly) return FindElementIdsForPolygon(poly.Coordinates);
        // Add other cases as needed
        throw new NotSupportedException($"Geometry {g.Type} not supported.");
    }
    private int FindElementIdForPoint(Position p) { /* … */ throw new NotImplementedException(); }
    private IEnumerable<int> FindElementIdsForPolygon(IList<List<Position>> rings) { /* … */ throw new NotImplementedException(); }
}

This illustrates how the mesh from SpatialDefinitionEqD2 anchors all spatial queries.

---

Working with features

Build a FeatureCollection with attributes

var fc = new FeatureCollection();
fc.Attributes.Add(new Attribute("Value", typeof(double), 1));
fc.Attributes.Add(new Attribute("Quality", typeof(string), 16));
var poly = new Polygon();
poly.Coordinates.Add(new List<Position> {
  new(0,0), new(10,0), new(10,10), new(0,10), new(0,0)
});
var f = new Feature(poly);
f.AttributeValues["Value"] = 12.34;
f.AttributeValues["Quality"] = "Good";
f.Associations.Add(new Association<MyDomain>("ABC-123"));
fc.Features.Add(f);

Parse WKT -> Geometry -> Feature

var g = Geometry.FromWKT("LINESTRING (0 0, 10 0, 10 5)");
var f = new Feature(g);

Design notes

• Equality semantics:
  • Node<TId> and Element<TId,TNodeId> compare equal by Id only. If you create duplicate ids, they'll dedupe in Mesh.Nodes.
• SpatialDefinition coverage:
  • EqD2 and NonEqD2 are both fully implemented. CurveLinearD2 is a stub; calling its methods throws.
• BaseDataSource:
  • DateTimes, Items, SpatialDefinition, and both Get overloads are abstract and must be implemented. All other members have default implementations in the base class.
• Geometry coercion:
  • BaseGeometry.Coordinates is dynamic for flexibility; in your code, cast to the concrete type's expected collection (e.g., List<Position>).
• CRS handling:
  • CRS object is not enforced or interpreted by the core. It's a carrier; projection math is up to you (or your service).
• Web Mercator conversions:
  • Helpers assume WGS84 lon/lat inputs. Guard against invalid latitudes |lat| > 85 in production code to avoid infinities.

---

Patterns to reuse

• EqD2 centroids For a cell (ix, iy): centroid = (x0 + (ix + 0.5)*dx, y0 + (iy + 0.5)*dy).
• Polygon containment (simple) Quick filtering: use bounding boxes to preselect elements, then a point-in-polygon test on centroids before more expensive clipping.
• Area-weighted aggregation If you need polygon averages on unstructured meshes, precompute (elementIndex, weight) pairs once per polygon (weight = intersected area / polygon area), then apply across timesteps (see DFSU codebase for inspiration).

---

Quick recipes

Create a 100×100 m grid (20×10) and produce its mesh

var def = new SpatialDefinitionEqD2(x0: 500000, y0: 6200000, dx: 100, dy: 100, nx: 20, ny: 10);
var mesh = def.ToMesh();                // IMesh<int,int>
var elem = mesh.Elements[0];            // Element 1 (lower-left)
var nodeIds = elem.Nodes.Select(n => n.Id);   // e.g., 1,2,23,22

Convert Google tile bbox to lon/lat

var web = new BoundingBox(1.23e6, 6.54e6, 1.24e6, 6.55e6);
var geo = web.GoogleToLonLat();

WKT roundtrip

var poly = Geometry.FromWKT("POLYGON ((0 0, 5 0, 5 2, 0 2, 0 0))");
var wkt  = poly.ToWKT();

Extending the module

1. Implement more spatial definitions Add SpatialDefinitionCurveLinearD2:
   • Provide ElementCount, ElementCount(int), and ToMesh() with sensible node/element id schemes.
   • Keep node id monotonic and stable across ToMesh() calls to aid caching.
2. Richer geometry parsing Support GEOMETRYCOLLECTION in Geometry.FromWKT (delegate to existing parsers inside).
3. Robust spatial operations If you need precise intersections/containment, wrap bindings to NetTopologySuite in a higher layer; keep this core dependency-light.
4. DataSource helpers Provide common utilities (centroid lookup, bbox index, polygon rasterization) in a DHI.Spatial.Data.Helpers namespace to standardize how Get(...) is implemented across sources.

---

Reference: API surface (by area)

Geometry & helpers

• Position, BoundingBox, CoordinateReferenceSystem
• Point, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon, GeometryCollection
• Geometry.FromWKT, Geometry.TrimForOuterParenthesis, Geometry.SplitArray
• ExtensionMethods (Google/EPSG:3857, parsing)

Feature model

• IAttribute, Attribute
• IFeature, Feature
• IFeatureCollection, FeatureCollection
• IAssociation, Association<T>

Mesh & spatial definitions

• INode<TId>, Node<TId>
• IElement<TId,TNodeId>, Element<TId,TNodeId>
• IMesh<TElementId,TNodeId>, Mesh<TElementId,TNodeId>
• ISpatialDefinition, SpatialDefinitionType
• SpatialDefinitionEqD2, SpatialDefinitionNonEqD2, SpatialDefinitionCurveLinearD2 (stub)

Data sources

• IItem<TId>, Item<TId>, Item
• IDataSource<TItemId>
• BaseDataSource<TItemId>