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 <[email protected]> on behalf of Simon
SPDBA Greener <[email protected]>
Sent: 17 October 2022 01:56
To: [email protected] <[email protected]>
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, [email protected]<mailto:[email protected]> 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
[email protected]<mailto:[email protected]>
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) [email protected]<mailto:[email protected]>
_______________________________________________
postgis-users mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/postgis-users
