Module 4, Basic analyses with vector data

learning objective: learning how to perform spatial selections, build geometric relations, measure geometry properties, construct and manipulate geometries, and convert coordinate systems

introduction

This module covers the core spatial operations you will use in almost every GeoDMS project that involves vector data. Before diving in, it helps to understand the three conceptually distinct things you can do with vector data in GeoDMS:

Selections filter a domain to a subset based on an attribute condition. The geometry is not involved; only attribute values determine which rows are kept.

Spatial relations use geometry to link two domain units. The result is a relation attribute, just like the code-based relations from Module 1c, but now derived from the spatial position of features rather than from a shared key.

Geometry operations create or transform geometries: buffers, unions, intersections, differences, simplifications. These produce new geometries as their output.

All examples in this module use the same datasets introduced in Module 2A: NS_Stations_2019_RD (train stations, points), OSM_Motorways_NL (motorways, arcs), and CBS_COROP_2012 (COROP regions, polygons). Data is assumed to be in the Dutch national coordinate system RD New (EPSG:28992), with the rdc coordinate unit.

---

selections

A selection creates a new, smaller domain unit from an existing one, keeping only the rows that satisfy a condition.

select_with_org_rel

select_with_org_rel(condition) is the standard selection function. It takes a boolean condition and returns a new domain unit containing only the rows where the condition is true. A derived attribute org_rel is automatically generated, providing a relation back to the original domain.

// Select only COROP regions in the province of Noord-Holland (prov_code 27)
unit<uint32> corop_NH := select_with_org_rel(CBS_COROP/prov_code == 27[uint32])
{
   attribute<rdc>    geometry  (poly) := org_rel -> geometry;
   attribute<string> statnaam         := org_rel -> statnaam;
}

The -> operator is shorthand for a lookup via org_rel. The expression org_rel -> statnaam is equivalent to CBS_COROP/statnaam[org_rel]. Both are correct; the arrow notation is more readable.

Note that select_with_org_rel does not automatically copy any attributes to the new domain. You declare each attribute you need explicitly, using org_rel to look up the values from the original domain.

select_with_attr_by_cond

select_with_attr_by_cond(source_container, condition) is a more convenient variant: it automatically copies all attributes from the source container to the new domain, without you having to declare each one:

// Select stations with more than 50000 passengers per day, copy all attributes
unit<uint32> busy_stations := select_with_attr_by_cond(NS_Stations, NS_Stations/passengers > 50000[uint32]);

All subitems of NS_Stations that share its domain become subitems of busy_stations automatically. The tradeoff is that this function does not produce an org_rel, so you cannot easily look up additional attributes from the original domain afterwards.

When to use which:

	select_with_org_rel	select_with_attr_by_cond
Produces org_rel	Yes	No
Auto-copies attributes	No	Yes
Best for	Small selections, need org_rel for further joins	Large selections, want all attributes immediately

Reading tip: for a full comparison of all selection variants, see selection operator comparison

---

spatial relations

Spatial relations are relations derived from the geometric position of features rather than from shared attribute keys. They work identically to the code-based relations from Module 1c: the result is a relation attribute that can be used in lookups and aggregations.

point_in_polygon

point_in_polygon(point_attribute, polygon_attribute) returns, for each point, the index of the polygon it falls in. If a point falls in no polygon, the result is null. If a point falls in multiple overlapping polygons, the first match is returned (see point_in_ranked_polygon if you need to control this).

// For each train station, find the COROP region it is located in
attribute<CBS_COROP> corop_rel (NS_Stations) :=
    point_in_polygon(NS_Stations/geometry, CBS_COROP/geometry);

// Use the relation to bring COROP attributes to the station domain
attribute<string> corop_name (NS_Stations) := CBS_COROP/statnaam[corop_rel];

This pattern is the spatial equivalent of rlookup from Module 1c: where rlookup matches on equal attribute values, point_in_polygon matches on spatial containment.

point_in_all_polygons

When points can overlap with multiple polygons (for example, when polygons overlap, or when you want to find all municipalities within a buffer zone rather than just the one containing the centroid), use point_in_all_polygons:

// Find all (station, COROP region) combinations where the point lies within the polygon
// This produces a new domain with one row per matching (station, corop) pair
unit<uint32> station_corop_pairs := point_in_all_polygons(NS_Stations/geometry, CBS_COROP/geometry)
{
   attribute<NS_Stations> station_rel := first_rel;
   attribute<CBS_COROP>   corop_rel   := second_rel;
}

The result is a cross-domain unit where each row represents a point-polygon pair. This is the GeoDMS equivalent of a spatial many-to-many join.

dist: point-to-point distance

dist(point_set_A, point_set_B) calculates the Euclidean distance between corresponding points in two sets. If both arguments are attributes of the same domain, the result is the distance between each pair of points at the same index position:

// Distance from each station to a fixed reference point (e.g. Amsterdam Centraal)
parameter<rdc> amsterdam_cs := point_xy(121400.0, 487400.0, rdc);

attribute<m> dist_to_acs (NS_Stations) :=
    dist(NS_Stations/geometry, const(amsterdam_cs, NS_Stations));

For distances from every point to every other point (a full distance matrix), the impedance functions (such as impedance_matrix) are more appropriate, but those are covered in the network modules (Modules 6a and 6b).

---

geometry properties

These functions compute scalar or point values from existing geometries. They do not create or modify geometries.

area

area(polygon_attribute) returns the surface area of each polygon. The values unit of the result is derived from the coordinate unit:

attribute<m2> corop_area_m2 (CBS_COROP) := area(CBS_COROP/geometry);

// Convert to km2 for readability
attribute<km2> corop_area_km2 (CBS_COROP) := corop_area_m2[km2];

arc_length

arc_length(arc_attribute) returns the length of each arc (line):

attribute<m> road_length_m (OSM_Motorways) := arc_length(OSM_Motorways/geometry);

centroid and centroid_or_mid

centroid(polygon_attribute) returns the geometric centre of mass of each polygon as a point. For convex polygons this always falls inside the polygon. For concave or U-shaped polygons, the centroid can fall outside the polygon boundary, which causes problems if you subsequently use the centroid in a point_in_polygon call.

centroid_or_mid(polygon_attribute) solves this: it returns the centroid if it lies within the polygon, otherwise it returns a guaranteed interior point:

// Use centroid_or_mid to avoid the "centroid outside polygon" problem
attribute<rdc> corop_label_point (CBS_COROP) := centroid_or_mid(CBS_COROP/geometry);

// Now safe to use in point_in_polygon or for map labels

Use centroid_or_mid as your default. Reserve centroid for cases where you explicitly need the true geometric centre, knowing it may fall outside the shape.

bounding box

lower_bound and upper_bound return the minimum and maximum corner point of the axis-aligned bounding box of each geometry:

attribute<rdc> bbox_min (CBS_COROP) := lower_bound(CBS_COROP/geometry);
attribute<rdc> bbox_max (CBS_COROP) := upper_bound(CBS_COROP/geometry);

These are useful for quick spatial filtering before applying more expensive operations.

---

constructing geometry

Sometimes you need to build geometries from raw coordinate data rather than reading them from a file.

point_xy and point_yx

point_xy(x_attribute, y_attribute, coordinate_unit) constructs a point geometry from separate X and Y columns. This is the standard way to turn a CSV file with coordinate columns into spatial point data:

unit<uint32> addresses
:   StorageName     = "%ProjDir%/Data/addresses.csv"
,   StorageType     = "gdal.vect"
,   StorageReadOnly = "True"
{
   attribute<string>  postcode;
   attribute<float64> x_rd;
   attribute<float64> y_rd;

   // Construct a point geometry from the X and Y columns
   attribute<rdc>     geometry := point_xy(x_rd, y_rd, rdc);
}

point_yx(y, x, unit) is the equivalent for data sources where latitude comes before longitude (common in GPS / WGS84 data).

---

geometry operations

Geometry operations transform or combine geometries to produce new ones. Each operation creates a new geometry attribute on the same domain as the input, or in the case of overlay operations, on a new combined domain.

library choice: geos, bg, bp

The GeoDMS offers four geometry libraries: geos_ (GEOS), bg_ (Boost Geometry), cgal_ (CGAL), and bp_ (Boost Polygon). In general:

• Use geos_ functions as your first choice: they are the fastest, most reliable, and handle float64 coordinates natively.
• Fall back to bg_ if a specific geos_ variant is not yet available.
• CGAL is very precise, but much slower.
• Avoid bp_ functions in new configurations unless you specifically need to work with integer coordinates (spoint/ipoint). They are slower and limited to integer precision.

This module uses geos_ throughout; wherever a geos_ function does not yet exist for a specific operation, the bg_ equivalent is noted.

buffer

A buffer creates a new polygon at a specified distance around each input geometry. The buffer size is in the same unit as the coordinate system (metres for RD New):

// 500-metre buffer around each train station
attribute<rdc> station_buffer_500m (poly, NS_Stations) :=
    geos_buffer_multi_point(NS_Stations/geometry, 500.0, 16b);

The third argument (16b) is the number of angles used to approximate the circular arc. More angles give a smoother buffer but larger geometry. 16 is a practical default; use 8 for rough calculations on large datasets, 32 for presentation-quality output.

For polygon buffers:

// 1 km outward buffer around COROP regions
attribute<rdc> corop_buffer_1km (poly, CBS_COROP) :=
    geos_buffer_multi_polygon(CBS_COROP/geometry, 1000.0, 16b);

For arc buffers:

// 100-metre buffer corridor around motorways
attribute<rdc> motorway_buffer_100m (poly, OSM_Motorways) :=
    geos_buffer_linestring(OSM_Motorways/geometry, 100.0, 8b);

union and dissolve

geos_union_polygon merges all polygons in a domain into a single polygon. This is the equivalent of a “dissolve all” operation:

// Merge all COROP regions into one outline of the Netherlands
parameter<rdc> nl_outline (poly) :=
    geos_union_polygon(CBS_COROP/geometry);

To dissolve by group (for example, merge COROP regions to provincial level), use a relation:

// Dissolve COROP regions to province level, grouped by province_rel
attribute<rdc> province_geometry (poly, province) :=
    geos_union_polygon(CBS_COROP/geometry, CBS_COROP/province_rel);

difference

The difference operation subtracts one geometry from another. Using the - operator with float64 coordinates (fpoint/dpoint) automatically selects the geos_difference implementation:

// Cut the station buffer out of the COROP geometry (same-domain operation, void second domain)
attribute<rdc> corop_minus_buffer (poly, CBS_COROP) :=
    CBS_COROP/geometry - corop_buffer_for_corop;

For cross-domain operations, use geos_difference explicitly and ensure the domains are aligned or one is void (a parameter):

// Remove a single exclusion zone (parameter) from each COROP polygon
parameter<rdc> exclusion_zone (poly) := ...;

attribute<rdc> corop_clipped (poly, CBS_COROP) :=
    CBS_COROP/geometry - exclusion_zone;

intersection

The * operator (with float64 coordinates) automatically selects geos_intersect:

// Intersect COROP regions with a study area polygon
parameter<rdc> study_area (poly) := ...;

attribute<rdc> corop_in_study_area (poly, CBS_COROP) :=
    CBS_COROP/geometry * study_area;

For a full polygon overlay that creates a new domain with all unique polygon fragments from two layers, use geos_overlay_polygon:

// Overlay COROP with province boundaries to get all unique fragments
unit<uint32> corop_x_province := geos_overlay_polygon(CBS_COROP/geometry, province/geometry)
{
   attribute<CBS_COROP> corop_rel    := first_rel;
   attribute<province>  province_rel := second_rel;
}

simplify

Simplification reduces the number of vertices in a polygon or arc while preserving its general shape. This is useful for faster calculations or smaller output files when full precision is not needed:

// Simplify COROP boundaries with a 100-metre tolerance
attribute<rdc> corop_simple (poly, CBS_COROP) :=
    geos_simplify_multi_polygon(CBS_COROP/geometry, 100.0);

The tolerance argument is in the coordinate unit (metres for RD New). A larger tolerance means more simplification and fewer vertices.

---

coordinate conversions

GeoDMS models typically work in a single coordinate system throughout. When source data arrives in a different system, convert it at the point of reading.

The most common conversion in Dutch projects is between RD New (EPSG:28992) and WGS84 (latitude/longitude):

// Convert RD New coordinates to WGS84 latitude/longitude
attribute<latlong> station_wgs84 (NS_Stations) :=
    rd2latlongwgs84(NS_Stations/geometry);

// Convert WGS84 back to RD New
attribute<rdc> station_rd (NS_Stations) :=
    latlongwgs842rd(station_wgs84);

The result of rd2latlongwgs84 uses the latlong coordinate unit (a point type with Y = latitude, X = longitude, in degrees). Make sure to define this unit in your configuration:

unit<float64> degrees        := baseunit('deg', float64);
unit<fpoint>  latlong        := fpoint;  // Y = lat, X = lon, in degrees

For other conversions (RD to ED50, RD to Google Maps projection), see rd2latlonged50 and rd2latlongge. Next to using these dedicated coordinate conversion operators, you can convert like you would value units with []. For example:

// Convert WGS84 back to RD New
attribute<rdc> station_rd (NS_Stations) := station_wgs84[rdc];

---

see also

The following functions are related but not covered in detail here:

• point_in_ranked_polygon: like point_in_polygon but lets you specify a ranking attribute to control which polygon is returned when multiple overlap.
• join_near_values: creates a new domain unit pairing points with all nearby points within a given distance. Useful for proximity analyses between two point sets.
• sequence2points: extracts all vertices of each arc or polygon as individual points.
• points2sequence: constructs arc or polygon geometry from a sequence of points.
• bg_simplify_linestring: simplifies arc geometries (the line equivalent of polygon simplification).

---

try it yourself! [TO BE ADDED!]

• Download the project here if you have not done so already, and unzip it to C:/prj/GeoDMSAcademy.
• Open exercise.dms (in the GeoDMS_Academy/basic_analyses_vector/cfg subfolder) in your text editor. The configuration already includes the source data from Module 2A.

Work through the following exercises:

1. Selection: select all COROP regions whose area exceeds 2000 km². Compute the area first using area, then apply select_with_org_rel. How many regions remain? Before you compare against 2000, check your units: open the area attribute in the Detail Pages and verify that its values unit and metric read km², not m². Comparing the m² attribute against a bare 2000 keeps every region, and the metric on the generic tab is the fastest way to spot that.
2. Spatial relation: use point_in_polygon to find the COROP region for each train station. Then aggregate: count the number of stations per COROP region using sum and a constant attribute.
3. Constructing geometry: the CSV file addresses.csv in the data folder has columns x_rd and y_rd. Use point_xy to add a geometry attribute, then map the addresses in the GUI.
4. Buffer and overlay: create a 5 km buffer around each train station. Use point_in_polygon to find which COROP region each buffer centroid falls in (use centroid_or_mid on the buffer geometry). Compare this with the direct point_in_polygon result from exercise 2.
5. Dissolve: merge all COROP regions that belong to the same province into a single province polygon using geos_union_polygon. Visualise the result alongside the original COROP boundaries.

The reference solution is in result.dms in the same cfg subfolder.

---

For other examples, see:

• Spatial Overlap
• Buffer processing example
• Polygon processing example
• convex hull example script

---

Go to previous module: Module 3, Meta scripting ‐ templates and for_each

Go to next module: Module 5, Basic analyses with grid data