ogr2ogr -segmentize mishandles Z values

[jidanni]

Feature description

One would never know this from the ogr2ogr man page, but ogr2ogr -segmentize mishandles Z values.

[jratike80]

The issue exlained without links is this:

ogr2ogr -segmentize does not interpolate Z values. Probably not M values either, and that behavior is not documented in https://gdal.org/en/stable/programs/ogr2ogr.html#cmdoption-ogr2ogr-segmentize.

LINESTRING Z (0 0 0,8 8 8) with -segmentize 3 becomes LINESTRING Z (0 0 0,2 2 8,4 4 8,6 6 8,8 8 8)

It is hard to say if ogr2ogr mishandles Z because the intented behavior is not documented.

With -segmentize 1 this linestring LINESTRING Z (0 0 0,2 2 2,4 4 4,8 8 8) becomes LINESTRING Z (0 0 0,0.666666666666667 0.666666666666667 2,1.33333333333333 1.33333333333333 2,2 2 2,2.66666666666667 2.66666666666667 4,3.33333333333333 3.33333333333333 4,4 4 4,4.66666666666667 4.66666666666667 8,5.33333333333333 5.33333333333333 8,6 6 8,6.66666666666667 6.66666666666667 8,7.33333333333333 7.33333333333333 8,8 8 8)

So it seems that -segmentize selects a closest Z value by some algorithm. It is not necessarily wrong, interpolating it not always the best thing to do with measurement data that is stored into Z.

However, PostGIS does interpolate Z

select st_astext(st_segmentize(
st_geomfromtext('LINESTRING Z (0 0 0,8 8 8)'),3));
"LINESTRING Z (0 0 0,2 2 2,4 4 4,6 6 6,8 8 8)"

And so does SpatiaLite, with a bit different strategy for creating new points:

LINESTRING Z(0 0 0, 2.12132 2.12132 2.12132, 4.242641 4.242641 4.242641, 6.363961 6.363961 6.363961, 8 8 8)

There is a workaround that is presented in the gis.stackexchange question, and it is to use the -SQLite SQL dialect and ST_Segmentize function. -dialect SQLite -sql "SELECT ST_Segmentize(geometry, 3)..."

[jidanni]

It must be wrong, because it is introducing bends in what should be a straight line.

The bends are not present in the other leading brands' solutions, but are revealed once one looks from any other angle than straight down or up.

How to fix

The program surely uses the same algorithm for both X and Y. It is a good algorithm. Simply also use it for Z. Err, well, something like that.

[jratike80]

I agree that users, including me, would prefer interpolated Z coordinates. I would still not say it "must" be wrong. For example, in the height contour data from our agency the interpolated points have Z value set to -9999, because the DEM specialists say that there is no guarantee that those interpolated heights are accurate, like the ones which have been measured with stereo workstations, and users have right to know it. That naturally makes it hard to use the densified contours with standard 3D viewers but I would still not say that the experts must be wrong. Their way to think does have a point, too.

In case of GDAL, maybe the algorithm is very old. PostGIS did not interpolate Z values until year 2006 but set them into 0.0 https://postgis.net/docs/ST_LineInterpolatePoint.html.

Also I think that the algorithm for interpolating Z in segmetizing/densification is not the same than the one used for X and Y. See the PostGIS documentation https://postgis.net/docs/ST_Segmentize.html.

Returns a modified geometry/geography having no segment longer than max_segment_length. Length is computed in 2D.

It does not really matter what the algorithm is, but it might make the developers feel better if they do not need to read comments like "Simply also use it for Z" if it is something that cannot be done that way.

EDIT:

Adding historical notes:

• GDAL 1.1.1 introduced this change: " * OGRLineString::segmentize() : do not set 0 as z for interpolated points, but the z from the previous point"
• the description in https://github.com/OSGeo/gdal/blob/c8bc3bb594773d03d45ad55971418f569497b743/ogr/ogrgeometry.cpp#L870 gives outdated information "Interpolated points will have Z and M values (if needed) set to 0"

[theroggy]

For information, the st_segmentize spatialite function uses the rttopo library under the hood, which is +- a fork of the relevant code of postgis of some years ago.

Spatialite also has a function GEOSDensify that does the same but uses the GEOS implementation. In that implementation handling of Z coordinates has recently been improved: in earlier versions they were +- ignored.