How to reduce the clothesline-effect in ST_Transform?

Question

When we use standard ST_Transform (e.g. PostGIS implementation), from a SRID to another that curves straight lines, the original straight lines represented by two points, that must be transformed into curves, will be transformed into straight lines. This is the central problem.

Illustrating below a real-life example (from an equal-area projection to WGS84), a grid transformation, where the large cell, after transform, loses compatibility with the small cells.

I don't know if the problem already had an appropriate name, but everyone knows the curve made by the line of a clothesline supported between two points, and they know that tightening it, it will be a straight line.

The set of original straight lines are like a tight clothesline: if the ST_Transform is made with more points (e.g. by ST_Segmentize) it results in a correct transformation, analog to a loose clothesline.

[Theorem] When not good, that is, when ST_Transform results are bad:
the central distance d (illustrated) between "correct line" and "bad line" is proportional to the length of the original line.
So, seems that a "good segmentization procedure" must use two reasonable assumptions:

• the optimal split is into the middle (where d is maximum);
• the length of the segments must be a fraction of the "characteristic diameter" (cdiam) of the geometry that will be transformed.

[The question] How to do "good segmentization"? If the hypotheses are correct, we have three problems:

1. more segments is more cost: more CPU time and more "pollution" (more points and geometrical complexity), we need a balance between quality and complexity.

2. segmentize is not necessary in some directions: in general a direction (for example lines parallel to longitude) not need to be segmentized.

3. aleatory segmentization is not good: the "split into segments" process must to use (the optimal) middle point.

... A naive solution is to use a fraction of cdiam as parameter in ST_Segmentize.

---

A clue but not the optimal solution

This example was expressed in PostGIS, but is only an illustration, the ideia is to segmentize by "characteristic diameter".

CREATE FUNCTION ST_Transform_resilient(
  g geometry,
  srid integer,
  density float DEFAULT 0.05
) RETURNS geometry AS $f$
  SELECT CASE
    WHEN isnot_point THEN -- using the fraction of cdiam here:
          ST_Transform(  ST_Segmentize(g,density*(a+p)/2.0)  ,srid)
    ELSE  ST_Transform(g,srid)
    END
  FROM (
    SELECT CASE WHEN is_poly THEN SQRT(ST_Area(g))
                ELSE 0 END,
           CASE WHEN is_poly THEN ST_Perimeter(g)/3.5
                WHEN isnot_point THEN 2.0*ST_Length(g) END,
           isnot_point
    FROM ( SELECT GeometryType(g) ) t1(geomtype), LATERAL (
      SELECT
         CASE WHEN geomtype='POLYGON' OR geomtype='MULTIPOLYGON' THEN true ELSE false END AS is_poly,
         CASE WHEN geomtype='POINT'   OR geomtype='MULTIPOINT'   THEN false ELSE true END AS isnot_point
    ) t2
  ) t3(a,p,isnot_point)
$f$ LANGUAGE SQL IMMUTABLE;

PS: a final correction is possible by minimal ST_Simplify, to remove collinear points. I need an optimal solution, one that generates good geometries with less complexity, and uses less CPU time to do so... And, ideally, understands the directions that really need ST_Segmentize.

Answer 1

After 2 years, no answer. Let's use the question's clues.

The answer to "How to do good segmentization?" is the function bellow. It is good for save disk space or reduce bandwidth in APIs... But but it is not the fastest solution for APIs. To make it faster, you need to adapt the function to your SRIDs (depends on source and destination), simplifying it, or developing a custom solution.

Solution for save disk space and bandwidth

The first good clue was to use the "characteristic diameter" of the geometry as reference in the ST_Segmentize() function.

CREATE FUNCTION ST_CharactDiam(g geometry)
 RETURNS float
AS $f$
   SELECT CASE
       WHEN tp IS NULL OR tp IN ('POINT','MULTIPOINT') THEN 0.0
       WHEN is_poly AND poly_p<2*poly_a THEN (poly_a+poly_p)/2.0 -- perimeter
       WHEN is_poly THEN (2*poly_a+SQRT(poly_p))/3.0  -- fractal perimeter
       ELSE ST_Length(g)/2.0  -- or use buffer
     END
   FROM (
     SELECT tp, is_poly,
            CASE WHEN is_poly THEN SQRT(ST_Area(g)) ELSE 0 END AS poly_a,
            CASE WHEN is_poly THEN ST_Perimeter(g)/3.5 ELSE 0 END AS poly_p
     FROM (
       SELECT tp, CASE
                WHEN tp IN ('POLYGON','MULTIPOLYGON') THEN true
                ELSE false
              END is_poly
       FROM (SELECT GeometryType(g)) t1(tp)
   ) t2
) t3;
$f$ language SQL IMMUTABLE;

The second good clue was to generalize the solution by a ST_Transform() extension, that here was adapted to apply the ST_SimplifyPreserveTopology() function after transformation.

CREATE FUNCTION ST_Transform_Resilient(
  g    geometry,     -- the input geometry
  srid integer,      -- target SRID, to transform g
  size_fraction float DEFAULT 0.05,  -- 1/density. 
   -- Density of points per charactDiam (or negative for absolute fraction).
  tolerance     float DEFAULT 0      -- E.g. on srid=4326 use 0.00000005. 
   -- ZERO=0.000000001. Good from 0.000000005 to 0.00000002.
) RETURNS geometry
AS $f$
 SELECT CASE
      WHEN COALESCE(size_fraction,0.0)>0.0 AND COALESCE(tolerance,0)>0 THEN
           ST_SimplifyPreserveTopology(geom,tolerance) -- ST_Simplify enough for grid 
      ELSE geom
      END
 FROM (
  SELECT CASE
    WHEN size>0.0 THEN  ST_Transform(  ST_Segmentize(g,size)  , srid  )
    ELSE  ST_Transform(g,srid)
    END geom, size
  FROM (
    SELECT CASE
         WHEN size_fraction IS NULL THEN 0.0
         WHEN size_fraction<0       THEN -size_fraction
         ELSE         ST_CharactDiam(g) * size_fraction
         END
  ) t1(size)
 ) t2
$f$ language SQL IMMUTABLE
;
COMMENT ON FUNCTION ST_Transform_Resilient IS
  'See problem/solution discussed at https://gis.stackexchange.com/q/444441'
;

---

Example appling in a equal-area-Albers grid, using ST_Transform_Resilient(geom,4326,0.05,0.00000005)... We can obtain the number of points of a geometry, before and after transform, by ST_NPoints(geom).

Here all geometries has 4 points (cells of square grid) and densified to 749 points. Starting from a 14-cell grid (hLevel=0), and subdividing each cell by 4 to form the next level, we obtain a table with the folowing columns: hLevel is the hierarchical level of the grid, n_cells the number of grid-cells, side_km the side-size in kilometers and nPts_simpl the number of points after simplification.

hLevel	n_cells	side_km	nPts_simpl
0	14	262	293
1	56	131	152
2	224	66	77
3	896	33	39
4	3584	16	20
5	14336	8	11