8

I have, in PostGIS, a bunch of LINESTRING Z objects (which represent poles), and I want to find which ones are inside an area (represented by a POLYGON). For the purpose of this exercise, we can safely assume that a pole is pretty much vertical, so it won't intersect the boundary of the area.

The problem is that sometimes the pole is exactly vertical.

This query, the one I would like to do, does not succeed:

SELECT ST_Intersects(ST_GeomFromText('LINESTRING Z (544483.525 6849134.28 104.1098,544483.525 6849134.28 114.6)',28356),
ST_GeomFromText('POLYGON((543907.636214323 6848710.84802846,543909.787417164 6849286.92923919,544869.040437688 6849283.30837091,544866.842236582 6848707.22673193,543907.636214323 6848710.84802846))',28356));

Exploring a few variants, this succeeds:

SELECT ST_Intersects(ST_GeomFromText('POINT (544483.525 6849134.28)',28356),
ST_GeomFromText('POLYGON((543907.636214323 6848710.84802846,543909.787417164 6849286.92923919,544869.040437688 6849283.30837091,544866.842236582 6848707.22673193,543907.636214323 6848710.84802846))',28356));

This does not succeed:

SELECT ST_Intersects(ST_GeomFromText('LINESTRING (544483.525 6849134.28,544483.525 6849134.28)',28356),
ST_GeomFromText('POLYGON((543907.636214323 6848710.84802846,543909.787417164 6849286.92923919,544869.040437688 6849283.30837091,544866.842236582 6848707.22673193,543907.636214323 6848710.84802846))',28356));

If you do a 3DIntersects query, it succeeds:

SELECT ST_3DIntersects(ST_GeomFromText('LINESTRING Z (544483.525 6849134.28 104.1098,544483.525 6849134.28 114.6)',28356),
ST_GeomFromText('POLYGON((543907.636214323 6848710.84802846,543909.787417164 6849286.92923919,544869.040437688 6849283.30837091,544866.842236582 6848707.22673193,543907.636214323 6848710.84802846))',28356));

however it raises a NOTICE that:

One or both of the geometries is missing z-value. The unknown z-value will be regarded as "any value"

So that works, sort of, but it fills my logs with noise and I'd prefer not to turn warnings off.

My reading of the OpenGIS Simple feature access standard (see part 1, section 4.14) is that even though the LINESTRING is a degenerate one, it should still be interpreted as topologically closed, so they should intersect.

Is this correct?

EDIT

After some experimentation, I can answer part of my own question. The reason why the above queries return false is this because the LINESTRING Z is believed not to be valid:

SELECT ST_IsValidReason(ST_GeomFromText('LINESTRING Z (544483.525 6849134.28 104.1098,544483.525 6849134.28 114.6)',28356));
                           st_isvalidreason                           
----------------------------------------------------------------------
 Too few points in geometry component[544483.525 6849134.28 104.1098]
(1 row)

This is a polite way of saying that OGC doesn't support 3D, and as such, GDAL/PostGIS only kind-of supports it.

I can accept that, although I still haven't found anything in the OGC specification which indicates that the corresponding LINESTRING is not valid.

So I guess my question really is: Is there an officially blessed way to find the footprint of a piece of 3D geometry which always returns valid 2D geometry?

5
+50

A really simple fix on the first query is to use ST_MakeValid to create a "valid" geometry to return true:

SELECT ST_Intersects(ST_MakeValid(A), B), ST_AsText(ST_MakeValid(A)) AS ST_MakeValid
FROM (
  SELECT ST_GeomFromText('LINESTRING Z (544483.525 6849134.28 104.1098,544483.525 6849134.28 114.6)',28356) AS A,
         ST_GeomFromText('POLYGON((543907.636214323 6848710.84802846,543909.787417164 6849286.92923919,544869.040437688 6849283.30837091,544866.842236582 6848707.22673193,543907.636214323 6848710.84802846))',28356) AS B
) AS d;

-[ RECORD 1 ]-+-----------------------------------------
st_intersects | t
st_makevalid  | POINT Z (544483.525 6849134.28 104.1098)

The original query failed because the LineString object was not valid, as discovered by either ST_IsValid or ST_IsValidReason. This is because ST_Intersects is a 2D operator, so the geometry was forced to two dimensions, LINESTRING(544483.525 6849134.28,544483.525 6849134.28), which is invalid.

ST_MakeValid removes the repeated coordinate (in 2D space), and changes it to a point type since the resulting geometry has only one coordinate.

You may see that the Z-coordinate is partially still there, but it's ignored. If this bothers you, you can also use ST_Force2D to just use 2D geometries for things like ST_Intersects:

SELECT ST_Intersects(ST_MakeValid(ST_Force2D(A)), B),
  ST_AsText(ST_Force2D(A)) AS ST_Force2D,
  ST_AsText(ST_MakeValid(ST_Force2D(A))) AS ST_MakeValid
FROM (
  SELECT ST_GeomFromText('LINESTRING Z (544483.525 6849134.28 104.1098,544483.525 6849134.28 114.6)',28356) AS A,
         ST_GeomFromText('POLYGON((543907.636214323 6848710.84802846,543909.787417164 6849286.92923919,544869.040437688 6849283.30837091,544866.842236582 6848707.22673193,543907.636214323 6848710.84802846))',28356) AS B
) AS d;

-[ RECORD 1 ]-+--------------------------------------------------------
st_intersects | t
st_force2d    | LINESTRING(544483.525 6849134.28,544483.525 6849134.28)
st_makevalid  | POINT(544483.525 6849134.28)
  • Thanks! This was a good explanation of what's going on, and the workaround is also good. – Pseudonym May 8 at 2:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.