[PROJ] Synthetic Coordinate System

[Ben Griffin]

(List: Apologies for the earlier double-post: 40 years of emails and I can still double up To/Cc)
Thomas, Javier, Martin, Even — thank you all for engaging with this so carefully.

Thomas, the UTM zone analogy is exactly right. Just as UTM zones share a common projection method but differ by zone parameter, the Hex9 planar projections share a common octahedral base but differ in how that base is unfolded or arranged in 2D. The analogy holds further: you wouldn't normally express UTM Zone 30N as derived from Zone 29N — and similarly, butterfly and mortar are independent renderings of the same base, not transformations of each other. It was the UTM Zone model that empowered me to venture an octant index.

Javier — you're correct that other polyhedral projections in PROJ return (x, y) without face information. The difference is that those projections define a global continuous mapping from sphere to plane, so the face is genuinely an internal detail.
Hex9 doesn't do that — and deliberately so. Each octant face maps to the same 2D coordinate range. A point in NWA and a point in SEP can have identical (x, y) values; without the octant index the coordinate is ambiguous. The octant is not a hint about which face was used internally — it's part of the coordinate itself, exactly as Thomas's UTM zone number is part of the easting. Stripping the octant from the output would be like returning bare eastings from a UTM projection without the zone — technically a number, but not a locatable coordinate.

The planar projections (butterfly, mortar etc.) resolve this by consuming the octant index to place each face unambiguously on the plane via per-face rigid transforms — that's exactly the "algorithm on top" you're describing. The CRS layer supplies (x, y, octant); the planar layer consumes all three and returns a globally unambiguous (u, v). The separation is intentional.

To make this concrete, the two relevant CRS definitions as registered in my current model of PROJ:

  HEX9:Oct — the intermediate octahedral barycentric CRS:
  PROJCRS["Hex9 Octahedral",
    BASEGEOGCRS["WGS 84", ...EPSG:4326...],
    CONVERSION["Hex9 Octahedral", METHOD["PROJ h9_boct", ID["HEX9","h9_boct"]]],
    CS[Cartesian,3],
    AXIS["hex9 x (X)", east,       LENGTHUNIT["metre",1]],
    AXIS["hex9 y (Y)", north,      LENGTHUNIT["metre",1]],
    AXIS["hex9 octant (O)", unspecified, SCALEUNIT["unity",1]],
    ID["HEX9","Oct"]]

  HEX9:butterfly — one of several planar projection CRS that consumes HEX9:Oct
  PROJCRS["Hex9 Butterfly",
    BASEGEOGCRS["WGS 84", ...EPSG:4326...],
    CONVERSION["Hex9 Butterfly", METHOD["PROJ h9_net_butterfly"]],
    CS[Cartesian,2],
    AXIS["net x", east,  LENGTHUNIT["metre",1]],
    AXIS["net y", north, LENGTHUNIT["metre",1]],
    ID["HEX9","butterfly"]]

The architecture is visible in the contrast: HEX9:Oct has a 3-axis Cartesian CS — x, y, and octant as a unity-scale third axis, exactly as  Thomas's UTM zone analogy suggests. HEX9:butterfly has a conventional 2-axis Cartesian CS — the h9_net_butterfly method has consumed the octant internally and produced a globally unambiguous planar coordinate. Both are fully resolvable by PROJ and usable in cs2cs pipelines.

However the WKT as stated above raises further questions:

Martin's points about ISO 19111 prompted me to take a closer look at the HEX9 WKT definitions, and I think it reveals exactly an architectural question he was describing.

Currently HEX9:butterfly is registered as:

  PROJCRS["Hex9 Butterfly",
    BASEGEOGCRS["WGS 84", ...EPSG:4326...],
    CONVERSION["Hex9 Butterfly", METHOD["PROJ h9_net_butterfly"]],
    CS[Cartesian,2], ...]

This works in PROJ — h9_net_butterfly is a single plugin that chains the octahedral barycentric step and the net unfold step internally.But it misrepresents the architecture: Butterfly is not derived from WGS 84, it is derived from HEX9:Oct.

The conceptually correct expression would be a DERIVEDPROJCRS with BASEPROJCRS:

DERIVEDPROJCRS["Hex9 Butterfly",
  BASEPROJCRS["Hex9 Octahedral", ...HEX9:Oct...],
  DERIVINGCONVERSION["Hex9 Butterfly", METHOD["PROJ h9_planar_butterfly"]],
  CS[Cartesian,2], ...]

This would properly express the two-layer architecture: WGS 84 → HEX9:Oct (ISO 19111 CRS layer) → HEX9:butterfly (planar rendering layer).

Does PROJ support DERIVEDPROJCRS with custom methods, and is this the recommended form? 
The reason it matters (beyond tidiness) is that expressing the planar layouts as DERIVEDPROJCRS from HEX9:Oct formally establishes the octahedral CRS as a stable base from which multiple independent projections can be derived — butterfly, mortar, windmill and others are then alternative renderings of the same CRS rather than separate projections from the ellipsoid. Any future projection of the octahedral coordinate space would derive from the same base, without revisiting the ellipsoid relationship, and facilitate a far more elegant process when eg, projecting from ‘butterfly’ to ‘mortar’.

Alternatively, since the planar forms are per-face rigid transforms, does PROJ accept DERIVEDPROJCRS in text_definition? 
If so, PROJ handles the chaining internally and we may not need h9_net_* combined plugins at all — just h9_planar taking a +layout= parameter.

Does the 3D→2D transition work? 
While HEX9:Oct is 3-axis (x, y, octant); the DERIVINGCONVERSION consumes all three and outputs 2D (x, y). 
h9_planar already does exactly this, but PROJ needs to be happy with the CS dimension change.

Best,
- Ben



> On 24 Mar 2026, at 10:00, Thomas Knudsen <knudsen.thomas at gmail.com> wrote:
> 
>> The octant component is ordinal/nominal, not a third spatial axis — there's no
>> meaningful interpolation between octant 3 and octant 4.
>> I'm encoding it as a 3rd Cartesian axis with orientation=unspecified and
>> SCALEUNIT["unity",1] because OrdinalCRS isn't composable with
>> ProjectedCRS in current PROJ. Is there a better-supported pattern for this?
> 
> Perhaps a variant of the German "UTM zone 31N/32N/33N with zone prefix"
> (CRS EPSG:5651/52/53) using false eastings of 31500000, 32500000, 33500000,
> respectively, to force the first two digits of the easting to signify
> the zone number
> 
> /Thomas


> I think I see your point. I have been recently developing some code related to GBT and DGGS.
> What I do not completely understand is how the face is part of the projection, specially in PROJ. Yes, it can be easier to implement. However all the other polyhedral projections in PROJ are doing the forward and inverse transformation without relying on that information.
> I think the projection in PROJ should be like the other ones, taking lat and lon, and return x and y (and viceversa).
> Deciding on with face is the point and what to do with it is a task for the algorithm on top of it.
> 
> Cheers,
> Javier
>>