[postgis-users] Issue with ST_Intersection seems related to lines crossing themself
Rory Meyer
rory.meyer at VLIZ.be
Mon Oct 17 06:09:29 PDT 2022
Hi Juergen,
I've been helping on an open source project for processing AIS data in postgis with a focus on scientific applications instead of commercial ones: http://openais.xyz/ It's still a little immature and the database portion of the project is still closed while I rewrite some of the containers.
Somethings that have helped me:
* Since AIS is a timeseries, geospatial dataset it often helps to have a timeseries plugin to work with. The TimeScaleDB HA version contains both PostGIS 3+ and TimescaleDB plugins. Timescaledb lets you automatically partition a table based on time and to automatically create binned data from live data. You could have a first/last position per ship MMSI per hour/day/week.
* The GPS on AIS is generally pretty accurate but can suffer from poor GPS lock when in ports that have high buildings, cranes etc near the harbour. This causes the vessels to "jump" around much more than expected. There are also instances of the vessels jumping several hundred km's in either the latitude or longitude direction. This might be due to bit-flip errors in the AIS VHF message.
* I would say that having a line intersect with your port geometry, even if no points are within it, is the kind of behaviour that you want. If this is not the case, then I would just check for points within the port area instead of the trajectory.
* For handling a trajectory that is on the edge of the port area and causing multiple, close in time, entrance and exit signals it might be worth looking into some of the following:
* Some kind hysteresis function by having a combination of the port boundaries and a buffer of the port boundaries.
* Filtering the entrance/exit signals through time. TimeScaleDB gives aggregate functions like "first(<column>, timestamp)" that would select the first item in the column by timestamp. First grouping the data by some sensible time window (daily) and then selecting the first/last entrance/exit signal per vessel
* Ignoring entrance/exit signals where the reverse signal occurred within X seconds.
* You might need more complex Trajectory generation code, something that uses window functions to break trajectories on obviously bad points (gaps between messages longer than 1 hour, sudden jumps of 100km+, points where calculated speed > 2* reported speed)
* It might be worth dumping the Z-value of the trajectory to get the first/last time in the intersection segment.
* To make your SQL a little easier to read you can use Common Table Expressions:
```
WITH my_port AS
(SELECT 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 ),
my_traj AS
(SELECT
ais.mmsi,
ais.created_on_date as date,
ais.min_created_on as mintime,
ais.max_created_on as maxtime,
ais.geo_lat_long_line as geom --renaming these to just help me out here
FROM ais_trajectory_mmsi_daily AS ais
)
SELECT
my_traj.mmsi,
my_traj.date,
ST_Intersection(my_traj.geom, my_port.geom) as port_traj_intersect,
ST_ZMin(ST_Intersection(my_traj.geom, my_port.geom)) as segment_start_time,
ST_ZMax(ST_Intersection(my_traj.geom, my_port.geom)) as segment_end_time
FROM my_geom, my_traj
WHERE ST_Intersects(my_traj.geom, my_port.geom)
AND my_traj.mintime > '2022-10-01'
AND my_traj.maxtime < '2022-10-02'
```
Hope that helps some,
Rory
________________________________
From: postgis-users <postgis-users-bounces at lists.osgeo.org> on behalf of Simon SPDBA Greener <simon at spdba.com.au>
Sent: 17 October 2022 01:56
To: postgis-users at lists.osgeo.org <postgis-users at lists.osgeo.org>
Subject: Re: [postgis-users] Issue with ST_Intersection seems related to lines crossing themself
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<mailto: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<mailto: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<http://www.spdba.com.au>
(m) simon at spdba.com.au<mailto:simon at spdba.com.au>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20221017/6e7e1958/attachment.htm>
More information about the postgis-users
mailing list