PostGIS in Production: Hunting Zone Geometry at Scale

What a Zone Is

In Oma Riista, a zone (alue) is a user-defined hunting area. Hunters and hunting clubs compose zones from three distinct geometry sources: parcels from the Finnish National Land Survey cadastral registry (palstaalue), custom drawn polygons (zone_feature), and pre-defined Metsähallitus state forest hunting areas (mh_hirvi). A single zone can include hundreds of parcels, a handful of custom polygons, and one or more state forest blocks. The combined geometry is what gets stored, displayed on maps, and used for permit calculations.

The platform serves roughly 190,000 hunters. At peak — spring moose permit applications — thousands of zones are being read, recalculated, and rendered simultaneously. The queries in fi.riista.feature.gis.zone.query reflect that: they're heavily optimised for the specific geometry shapes and sizes of Finnish land parcels, with chunk-based processing and multi-level caching baked in from the start.

Building the Combined Geometry

CalculateCombinedGeometryQueries is the entry point for zone geometry computation. Its job is to merge all three geometry sources into a single unified polygon and persist it back to the zone table. The algorithm has four distinct phases.

Phase 1: Fetch and Sub-Divide

The first query pulls raw geometries from all three source tables and immediately passes them through ST_SubDivide:

WITH g AS (
  SELECT (ST_Dump(zp.geom)).geom AS geom
  FROM zone_palsta zp WHERE zp.zone_id = :zoneId
  UNION ALL
  SELECT (ST_Dump(zp.geom)).geom AS geom
  FROM zone_feature zp WHERE zp.zone_id = :zoneId
  UNION ALL
  SELECT (ST_Dump(mh.geom)).geom AS geom
  FROM zone_mh_hirvi zh
  JOIN mh_hirvi mh ON (zh.zone_id = :zoneId AND zh.mh_hirvi_id = mh.gid)
)
SELECT ST_AsBinary(ST_SubDivide(geom, 16384)) AS geom FROM g

ST_SubDivide is the key insight here. Finnish forest parcels are legal cadastral units — they follow property boundaries, not geographic simplicity. A single parcel can be an extremely complex polygon with thousands of vertices. Attempting ST_Union directly on a large set of such geometries is a reliable way to exhaust PostgreSQL's memory and produce a timeout. ST_SubDivide(geom, 16384) breaks each geometry into sub-polygons with at most 16,384 vertices each, which makes the subsequent union tractable.

The chunk size of 16,384 was not arbitrary — it reflects the observed vertex complexity of Metsähallitus state forest blocks, which are the most complex geometries in the dataset.

Phase 2: Union in the Application Layer

The sub-divided geometries are loaded into JTS (Java Topology Suite) geometry objects via WKB reader and merged using GISUtils.computeUnionFaster(). Doing the final union in Java rather than in a single ST_Union aggregate query gives finer control over memory usage and allows progress logging during long operations. The JTS union also applies ST_Buffer(geom, 0) on geometries that fail validity checks — a standard PostGIS idiom for repairing ring orientation and self-intersection artifacts that occasionally appear in cadastral data.

Phase 3: Subtract Excluded Areas

Zones support exclusions — geometry that the user explicitly wants removed from the union (for example, a water body or a privately posted area). After the union is computed, excluded geometry is subtracted with a small buffer:

ST_Difference(union_geom, ST_Buffer(excluded_geom, 0.5))

The 0.5-metre buffer on the excluded geometry ensures that adjacent-but-touching parcels don't produce slivers at the exclusion boundary. Without it, floating-point precision in parcel boundary coordinates reliably produces narrow polygon strips that are geometrically valid but semantically meaningless and visually confusing on the map.

Phase 4: Persist and Index

The final geometry is written back to zone.simple_geom in ETRS-TM35FIN (SRID 3067), Finland's national projected coordinate system. All area calculations in the system work in this coordinate system because ST_Area returns square metres directly when the geometry is in a metric projection — no unit conversion required.

Area Calculations Across Land Classifications

CalculateZoneAreaSizeQueries computes the zone's area broken down by land classification. This isn't a single query — it's a pipeline of five independent calculations that get assembled into a GISZoneSizeDTO.

Total and Water Area

Total area is straightforward: SELECT ST_Area(geom) FROM zone WHERE zone_id = :zoneId. Water area requires an intersection with the vesialue (water bodies) dataset:

WITH z AS (
  SELECT ST_SubDivide((ST_Dump(geom)).geom, 8192) AS geom
  FROM zone WHERE zone_id = :zoneId
)
SELECT SUM(ST_Area(ST_Intersection(va.geom, ST_Buffer(z.geom, 0))))
FROM z JOIN vesialue va ON ST_Intersects(va.geom, z.geom)

Two things to note. First, ST_Intersects is used as a pre-filter before ST_Intersection: ST_Intersects uses the spatial index and is fast; ST_Intersection computes the actual overlap geometry and is expensive. Always filter before you compute. Second, ST_SubDivide appears again with a chunk size of 8,192 — smaller than the 16,384 used in geometry combination because intersection queries are more memory-intensive than union queries on the same geometry.

State Land: A Six-Step Query

State land area (valtionmaa) is more involved because the requirement is to calculate state land area excluding water — i.e., the huntable state land. This requires a three-way intersection:

1. Sub-divide zone geometry into chunks
2. Intersect chunks with valtionmaa (state land dataset) — produces state land portions of the zone
3. Intersect those state land portions with vesialue (water bodies) — produces water within state land
4. Sum total state land area
5. Sum water within state land
6. Subtract: state_land_area - state_land_water_area = huntable_state_land

Private land area is then derived: total_area - water_area - state_land_area = private_land_area. This matters for hunting permit calculations because moose permit quotas in Finland are computed against the total huntable land area, with state and private land tracked separately for administrative reporting.

Area by Regional Hunting District

The most complex calculation is area broken down by RHY (riistanhoitoyhdistys, regional hunting district). A zone can span multiple RHY boundaries — a large club area frequently crosses district lines. The query uses a CTE chain to intersect the zone with each RHY polygon and then further intersect each piece with water and state land datasets:

WITH z AS (...subdivided zone...),
z_block AS (
  SELECT rhy.official_code,
    ST_MakeValid((ST_Dump(
      ST_Intersection(ST_SetSRID(rhy.geom, 3067), z.geom)
    )).geom) AS geom
  FROM z JOIN rhy ON ST_Intersects(rhy.geom, ST_SetSRID(z.geom, 3047))
)
-- further intersections with vesialue and valtionmaa, grouped by official_code

There's a subtle coordinate system issue in this query: the rhy.geom join condition uses SRID 3047, while the intersection itself uses 3067. This reflects a data provenance quirk — the RHY boundaries were originally ingested with a slightly different ETRS-based projection and the join predicate was written against that, while the computation uses the canonical national projection. The ST_SetSRID calls are cast-only (no re-projection), which works because the numeric difference between 3047 and 3067 coordinates in Finland is small enough that the spatial index lookup remains correct. It's the kind of thing that accumulates over years of real-world data management.

Inverted Geometry: Everywhere Except Here

GetInvertedGeometryQuery computes the complement of a zone — the geometry of Finland minus the zone polygon. This is used in the map UI to render a shaded overlay covering everything outside the hunting area, a common cartographic technique for focusing attention on a bounded region.

WITH finland AS (
  SELECT ST_Buffer(ST_Envelope(ST_GeomFromText(:worldBounds, 3067)), 0) AS geom
),
d AS (
  SELECT ST_Difference(
    ST_Transform(finland.geom, 4326),
    ST_Buffer(simple_geom, 0)
  ) AS geom
  FROM finland CROSS JOIN zone
  WHERE zone_id = :zoneId AND simple_geom IS NOT NULL
)
SELECT ST_AsBinary(ST_Transform(d.geom, :crs)) AS geom FROM d

The Finland bounding box is a hardcoded ETRS-TM35FIN envelope: LINESTRING(50199.4814 6582464.0358, 761274.6247 7799839.8902). Using a hardcoded national boundary rather than the actual Finnish coastline polygon is a deliberate trade-off — computing ST_Difference against the full 1:10,000 coastline geometry would be prohibitively slow. The rectangular approximation works fine for the cartographic purpose: the shaded area just needs to cover everything outside the zone. What it looks like in the Gulf of Bothnia is not a user concern.

State Land Geometry Extraction

GetStateGeometryQuery extracts the state-owned portions of a zone as a geometry object. It exists both as a database query (loading zone geometry from the zone table) and as a direct geometry operation (taking a JTS geometry parameter for real-time preview without a committed zone). The core CTE is reused across both variants:

z_state AS (
  SELECT (ST_Dump(
    ST_Intersection(vm.geom, ST_Buffer(z.geom, 0))
  )).geom AS geom
  FROM z JOIN valtionmaa vm ON ST_Intersects(vm.geom, z.geom)
  WHERE GeometryType(z.geom) IN ('POLYGON', 'MULTIPOLYGON')
)
SELECT ST_AsBinary(ST_Transform(
  ST_CollectionHomogenize(ST_CollectionExtract(ST_Collect(z_state.geom), 3)),
  :crs)) AS geom
FROM z_state

ST_CollectionExtract(geom, 3) filters the intersection result to polygon types only (type code 3). Geometry intersection can produce lower-dimensional results — lines and points where polygon edges or corners touch exactly — and those need to be discarded before the result is returned to the application layer. ST_CollectionHomogenize then normalises the collection to the most specific type that fits: if all pieces are polygons, you get a MultiPolygon; if there's only one, you get a Polygon. This prevents the application layer from having to handle GeometryCollection types, which JTS handles less cleanly.

Vector Tile Generation

The fi.riista.feature.gis.vector package handles real-time Mapbox Vector Tile (MVT) generation for zone rendering. When the map client requests tile z/x/y, VectorTileService clips the zone geometry to the tile's bounding box and encodes it as a protobuf binary.

Tile Mathematics

VectorTileUtil converts tile coordinates to a WGS84 bounding envelope using the standard Web Mercator inverse projection:

// Tile to longitude: straightforward linear interpolation
lon = x / (1 << z) * 360 - 180

// Tile to latitude: inverse Mercator
lat = Math.toDegrees(Math.atan(Math.sinh(Math.PI - 2 * Math.PI * y / (1 << z))))

The tile resolution in metres per pixel is derived from Earth's circumference: 40075017.0 / 256 / (1 << zoom). This resolution value drives the buffer applied around the tile envelope before clipping — the query clips to a slightly larger box than the tile itself to avoid rendering artifacts at tile edges where geometry crosses the boundary.

Geometry Clipping and Encoding

The database query clips zone geometry to the tile envelope:

SELECT ST_AsBinary(ST_ClipByBox2D(simple_geom,
  ST_MakeEnvelope(:xmin, :ymin, :xmax, :ymax, 4326))) AS geom
FROM zone WHERE zone_id = :zoneId

ST_ClipByBox2D is preferred over ST_Intersection with a box geometry for tile clipping specifically because it's faster: it uses a simplified clipping algorithm tuned for rectangular clip regions rather than the general polygon intersection algorithm. The trade-off is that it can produce invalid output geometries at the clip boundary, but for tile rendering this is acceptable — the client renderer handles edge cases gracefully.

After clipping, the Java code applies an affine transformation to convert geographic coordinates into tile pixel space (0–4096 extent, the MVT standard):

// Translate to origin, scale geographic to tile pixels, flip Y axis
AffineTransformation affine = new AffineTransformation()
    .translate(-xmin, -ymin)
    .scale(4096.0 / (xmax - xmin), -4096.0 / (ymax - ymin))
    .translate(0, 4096);

The Y-axis flip (-4096.0 / ... scale followed by a +4096 translate) is necessary because geographic coordinates increase northward while tile pixel coordinates increase downward. This transformation is applied in Java rather than SQL because it's cheaper to do linear math on already-loaded geometry objects than to round-trip through the database.

The transformed geometry is then serialised to Mapbox protobuf format via the no.ecc.vectortile encoder, placed in a single layer named "all", and returned as a byte array. The service logs a warning if the full pipeline exceeds 100ms or if the JDBC query alone exceeds 50ms — these thresholds were set based on observed p95 latencies under normal load and serve as an early warning for query plan degradation.

Caching Architecture

Zone ID lookups and authorization checks are cached in Guava LoadingCache instances with two distinct TTLs. The zone ID cache (mapping area ID to zone ID) expires after one minute — zone assignments change infrequently but aren't static. The authorization cache (mapping session ID + area ID to a read permission boolean) expires after five minutes. Authorization is checked against EntityPermission.READ per the standard AbstractEntityAuthorization pattern used throughout the system.

There are five distinct cache classes for different area types: personal areas, personal area unions (from harvest permit applications), hunting club areas, moderator-managed areas, and permit application areas. Each has its own cache instance rather than sharing a single cache keyed by area type, because their TTL and population patterns differ — permit application zones change rarely and can tolerate a longer TTL (five minutes vs. one minute for personal areas).

The ST_SubDivide + application-layer union pattern is the single most impactful architectural decision in the zone pipeline. Without it, combining a large Metsähallitus block with 300 cadastral parcels would reliably timeout. With it, the same operation completes in under two seconds.

What Makes This Stack Unusual

Most GIS applications at this scale either push all spatial operations into PostGIS or do all geometry processing in the application layer with a library like GDAL. The Oma Riista approach is deliberately hybrid: PostGIS handles set operations and intersections (where its R-tree spatial index makes a decisive performance difference), while JTS handles union operations (where iterative processing with progress control is more important than raw throughput).

The coordinate system discipline is strict throughout. Everything stored and computed uses ETRS-TM35FIN (SRID 3067). Transformations happen only at ingestion (external data coming in as WGS84) and at serving (vector tiles going out as WGS84). The one exception is the legacy SRID 3047 join condition in the RHY area query — a detail that will eventually be cleaned up, but which works correctly today because the spatial index tolerance covers the numerical difference.

The geometry validity repair pattern (ST_Buffer(geom, 0) in SQL, geometry.buffer(0) in JTS) appears throughout the codebase. Finnish cadastral data, like most real-world surveyed polygon datasets, contains geometries that are technically invalid under the OGC simple features specification — self-touching rings, duplicate points, incorrect winding order. Rather than rejecting these at ingestion, the system repairs them silently. This is pragmatic: the alternative is user-visible errors on valid land parcels, which is worse than a quiet zero-distance buffer.