[postgis-users] Issue with ST_Intersection seems related to lines crossing themself
Simon SPDBA Greener
simon at spdba.com.au
Sun Oct 16 16:56:36 PDT 2022
Would using a MultiPoint (or GeometryCollection with Points) for the GPS
locations rather than a Linestring be possible?
Simon
On 17/10/2022 10:10 am, magog002 at web.de wrote:
> Hi,
> I've got an issue with an SQL PostGIS query and hope to get a solution
> here.
> First some details.
> I'm storing GPS locations of vessels by their MMSI identification
> number (AIS) that are received in 3 Sek to 3 Minutes ranges in a table.
> As this becomes quite unhandy in a very short time, each day the basic
> details mmsi, received_at (Timestamp), lat, long are converted into a
> Points with the Timestamp integrated and this is then converted into a
> LineString for each MMSI, per day, with the first and last timestamp
> of the day for which a point exists as direct columns.
> Now I want to know when some vessel is/was inside a specific area
> (harbor, berth) so I have a Polygon which the LineString has to intersect.
> So far so good.
> This works but can create false positives (two points outside the
> polygon are connected and this created line goes through the Polygon
> but there is no actual Point of that line in the Polygon)...so I need
> to re-check again by pulling out the points from my LineString.
> Here the basic search query that basically works:
> ================================
> --
> -- Return all vessels (MMSI) found in the selected timeframe where the
> -- LineString Intersects the Polygon given as parameter to search.
> --
> SELECT ais_trajectory_mmsi_daily.mmsi
> ,ais_trajectory_mmsi_daily.min_created_on
> ,ais_trajectory_mmsi_daily.max_created_on
> ,ST_SetSRID(ais_trajectory_mmsi_daily.geo_lat_long_line, 4326)
> AS geo_lat_long_line
> FROM ais_trajectory_mmsi_daily
> WHERE
> ST_Intersects(ST_SetSRID(ais_trajectory_mmsi_daily.geo_lat_long_line,
> 4326),
> ST_GeometryFromText('POLYGON((7.1084729427821
> 50.7352828020046,7.1089771980769 50.7353778664669,7.10848259899649
> 50.7365542732199,7.10801482200623 50.7364673592064,7.1084729427821
> 50.7352828020046))', 4326) )
> AND min_created_on >= CAST('2022-10-01' AS TIMESTAMPTZ)
> AND max_created_on < CAST('2022-10-15' AS TIMESTAMPTZ)
> ;
> Here is now the query I use with ST_Intersection and where it starts
> to go wrong.
> When a vessel stays a the same place the GPS coordinates are the same
> + drift.
> As those points are all added to the LineString the line sections are
> running over themself and it looks like a "woolen noodle".
> This seems to cause issues with ST_Intersection as I use it because
> instead of getting only the entering and exiting details as normally 1
> row where the LineString hits my Polygon with times I now get hundreds
> of rows as a result and the visual result shown on the map is weird
> and sometimes wrong (the line goes out of the polygon in a way the
> original LineString did never exist).
> Here is the query in question (basically the part above) with ST_Dump
> and ST_Intersection in the SELECT part + plus an outer query to get
> more details on what was the first and last time there was an
> intersection.
> --
> -- Outer query with union to visualize the Polygon and the
> -- data received from the LineString that Intersects
> -- with the polygon via e.g. pgAdmin.
> --
> SELECT NULL AS mmsi
> ,NULL AS created_on_date
> ,NULL AS min_created_on
> ,NULL AS max_created_on
> ,ST_GeometryFromText('POLYGON((7.1084729427821
> 50.7352828020046,7.1089771980769 50.7353778664669,7.10848259899649
> 50.7365542732199,7.10801482200623 50.7364673592064,7.1084729427821
> 50.7352828020046))', 4326) AS geom
> ,NULL AS geom_start
> ,NULL AS geom_end
> UNION
> -- The actual data of interest
> SELECT intersection.mmsi
> ,intersection.created_on_date
> ,intersection.min_created_on
> ,intersection.max_created_on
> ,intersection.geom
> ,MIN(DATE_TRUNC('second',
> TO_TIMESTAMP(ST_InterpolatePoint(intersection.geo_lat_long_line,
> ST_StartPoint(intersection.geom)) )::timestamptz )) AS geom_start
> ,MAX(DATE_TRUNC('second',
> TO_TIMESTAMP(ST_InterpolatePoint(intersection.geo_lat_long_line,
> ST_EndPoint(intersection.geom)) )::timestamptz )) AS geom_end
> FROM (
> SELECT ais_trajectory_mmsi_daily.mmsi
> ,CAST(ais_trajectory_mmsi_daily.min_created_on AS date) AS
> created_on_date
> ,ais_trajectory_mmsi_daily.min_created_on
> ,ais_trajectory_mmsi_daily.max_created_on
> ,ST_SetSRID(ais_trajectory_mmsi_daily.geo_lat_long_line, 4326) AS
> geo_lat_long_line
> --
> -- Here seems to be the problem...
> --
> ,(ST_Dump(ST_Intersection(ST_GeometryFromText('POLYGON((7.1084729427821
> 50.7352828020046,7.1089771980769 50.7353778664669,7.10848259899649
> 50.7365542732199,7.10801482200623 50.7364673592064,7.1084729427821
> 50.7352828020046))', 4326)
> ,ST_SetSRID(ais_trajectory_mmsi_daily.geo_lat_long_line, 4326)
> )
> )).geom AS geom
> FROM ais_trajectory_mmsi_daily
> WHERE
> ST_Intersects(ST_SetSRID(ais_trajectory_mmsi_daily.geo_lat_long_line,
> 4326),
> ST_GeometryFromText('POLYGON((7.1084729427821
> 50.7352828020046,7.1089771980769 50.7353778664669,7.10848259899649
> 50.7365542732199,7.10801482200623 50.7364673592064,7.1084729427821
> 50.7352828020046))', 4326) )
> AND min_created_on >= CAST('2022-10-01' AS TIMESTAMPTZ)
> AND max_created_on < CAST('2022-10-02' AS TIMESTAMPTZ)
> ) AS intersection
> GROUP BY intersection.mmsi
> ,intersection.created_on_date
> ,intersection.min_created_on
> ,intersection.max_created_on
> ,intersection.geom
> ORDER BY created_on_date, mmsi, geom_start
> ;
> If needed I can try to create some sample with data to put on sqlfiddle.
> On Stackoverflow I found this which seems to be related to the same issue:
> https://stackoverflow.com/questions/62090829/why-st-intersection-from-postgis-returns-multilinestring-for-self-intersecting-l
> Many thanks in advance for hints to solve the problem.
> Juergen
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
--
Simon Greener
39 Cliff View Drive, Allens Rivulet, 7150, Tasmania, Australia
(m) +61 418 396 391
(w)www.spdba.com.au
(m)simon at spdba.com.au
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20221017/bff48927/attachment.htm>
More information about the postgis-users
mailing list