[prev in list] [next in list] [prev in thread] [next in thread]
List: postgis-users
Subject: Re: [postgis-users] Issue with ST_Intersection seems related to lines crossing themself
From: magog002 () web ! de
Date: 2022-10-31 15:08:04
Message-ID: trinity-90d1c190-45b4-4ac5-a8b2-c5e58bc379cf-1667228884494 () 3c-app-webde-bs22
[Download RAW message or body]
[Attachment #2 (text/html)]
<html><head></head><body><div style="font-family: Verdana;font-size: 12.0px;"><div>Hi \
Rory,</div>
<div> </div>
<div>thank you for that blog post. That was helpful.</div>
<div> </div>
<div>What I implemented also on friday was the query down below for the most recent \
data catched directly from the stored ais position reports. There are some of the \
checks missing you added into your example.</div>
<div>This is basically the part to get track from today (for older days this can be \
combined using UNION to get the details directly from the stored table with the \
stored MultiLinestrings + MultiPoints created before).</div>
<div> </div>
<div>In general the missing part is still the option to get back to the PointM (the \
stored timestamp) when using ST_Intersection because I need to look into the data of \
a complete month and looking for 1.8 Mio single points per day * 31 days of \
points is too slow. Using the linestring with ST_Intersecs from the daily stored \
data is much faster and I can get the first and last Point on that day from \
it...the the annoying issue of loosing the timestamp.</div>
<div>So basically when I have the mmsi, the points and I know at least the day \
and I need to do another lookup with only this to find that Point in my stored \
MultiPointM column. </div>
<div> </div>
<div> </div>
<div>--<br/>
-- Vessel tracks with separated LineStrings combined into Multilinestring</div>
<div>-- A track is separated if a point is > 7min away from the previous \
entry.</div>
<div>-- What's missing is the check that at least 2 points are required to build \
a linestring.</div>
<div>-- Also additional checks could be added.<br/>
--<br/>
WITH cte_point_m_pos_today AS (<br/>
--<br/>
-- Get data from today (since 00:00)<br/>
--<br/>
SELECT ais_positions_age_tagged.mmsi<br/>
,ais_positions_age_tagged.geom<br/>
,ais_positions_age_tagged.created_on<br/>
,ais_positions_age_tagged.ais_position_to_old<br/>
,SUM(ais_positions_age_tagged.ais_position_to_old) \
OVER (PARTITION BY ais_positions_age_tagged.mmsi, \
CAST(ais_positions_age_tagged.created_on AS DATE) ORDER BY \
ais_positions_age_tagged.created_on) AS ais_position_cluster_group<br/> \
FROM (<br/> SELECT mmsi<br/>
,created_on<br/>
,geo_lat_long_point AS geom<br/>
,COALESCE( CAST( created_on - \
(LAG(created_on) OVER (PARTITION BY mmsi, CAST(created_on AS DATE) ORDER BY \
created_on)) > CAST('00:07:00' AS INTERVAL) AS INTEGER), 0) AS \
ais_position_to_old<br/> FROM \
ais_position_report<br/> WHERE created_on >= CAST( \
CURRENT_DATE AS TIMESTAMPTZ)<br/> AND created_on \
< CAST( CURRENT_DATE AS TIMESTAMPTZ) + INTERVAL '1 day'<br/> \
AND mmsi = :the_mmsi_to_lookup<br/> \
) ais_positions_age_tagged<br/> ) -- CTE-1<br/>
SELECT mmsi<br/>
,MIN(min_created_on) AS min_created_on<br/>
,MAX(max_created_on) AS max_created_on<br/>
,ST_Collect(geom ORDER BY min_created_on) AS \
geo_lat_long_line<br/> FROM (<br/>
SELECT mmsi<br/>
,MIN(created_on) AS min_created_on<br/>
,MAX(created_on) AS max_created_on<br/>
,CAST(created_on AS DATE) AS created_on \
<br/> ,ST_MakeLine(geom ORDER BY mmsi, \
created_on) AS geom<br/> FROM cte_point_m_pos_today<br/>
GROUP BY mmsi, CAST(created_on AS DATE), ais_position_cluster_group<br/>
) AS my_daily_multilinestring_geom<br/>
GROUP BY mmsi, created_on<br/>
;</div>
<div>
<div>
<div name="quote" style="margin:10px 5px 5px 10px; padding: 10px 0 10px 10px; \
border-left:2px solid #C3D9E5; word-wrap: break-word; -webkit-nbsp-mode: space; \
-webkit-line-break: after-white-space;"> <div style="margin:0 0 10px \
0;"><b>Gesendet:</b> Freitag, 28. Oktober 2022 um 10:22 Uhr<br/> \
<b>Von:</b> "Rory Meyer" <rory.meyer@VLIZ.be><br/> \
<b>An:</b> "PostGIS Users Discussion" \
<postgis-users@lists.osgeo.org><br/> <b>Betreff:</b> Re: [postgis-users] \
Issue with ST_Intersection seems related to lines crossing themself</div>
<div name="quoted-content"><!--P {
margin-top: 0;
margin-bottom: 0;
}
-->
<div>
<div class="elementToProof ContentPasted0" style="font-family: Calibri , Arial , \
Helvetica , sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);">Hi Juergen,</div>
<div class="elementToProof ContentPasted0" style="font-family: Calibri , Arial , \
Helvetica , sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);"> </div>
<div class="elementToProof ContentPasted0" style="font-family: Calibri , Arial , \
Helvetica , sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);">I wrote a quick blog post about how I grouped up AIS points into \
lines, but broken up on steps in time, distance, and calculated speed. </div>
<div class="elementToProof ContentPasted0" style="font-family: Calibri , Arial , \
Helvetica , sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);"> </div>
<div class="elementToProof ContentPasted0" style="font-family: Calibri , Arial , \
Helvetica , sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);">It's very short but contains some thinking and some \
code.</div>
<div class="elementToProof ContentPasted0" style="font-family: Calibri , Arial , \
Helvetica , sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);"> </div>
<div class="elementToProof ContentPasted0" style="font-family: Calibri , Arial , \
Helvetica , sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);"><a href="https://openais.xyz/post/2022-10-28-ais-traj/" \
target="_blank">https://openais.xyz/post/2022-10-28-ais-traj/</a></div>
<div class="_Entity _EType_OWALinkPreview _EId_OWALinkPreview \
_EReadonly_1"> </div>
<div class="elementToProof" style="font-family: Calibri , Arial , Helvetica , \
sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);"> </div>
<div class="elementToProof" style="font-family: Calibri , Arial , Helvetica , \
sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);">You could add an additional spatial filter in the first part of \
the blog query to only fetch AIS points that are within the geometry that defines \
your port area.</div>
<div class="elementToProof" style="font-family: Calibri , Arial , Helvetica , \
sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);"> </div>
<div class="elementToProof" style="font-family: Calibri , Arial , Helvetica , \
sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);">Rory</div>
<div id="appendonsend"> </div>
<hr style="display: inline-block;width: 98.0%;"/>
<div id="divRplyFwdMsg"><font color="#000000" face="Calibri, sans-serif" \
style="font-size: 11.0pt;"><b>From:</b> postgis-users \
<postgis-users-bounces@lists.osgeo.org> on behalf of magog002@web.de \
<magog002@web.de><br/> <b>Sent:</b> 24 October 2022 19:40<br/>
<b>To:</b> postgis-users@lists.osgeo.org <postgis-users@lists.osgeo.org><br/>
<b>Subject:</b> Re: [postgis-users] Issue with ST_Intersection seems related to lines \
crossing themself</font>
<div> </div>
</div>
<div>
<div style="font-family: Verdana;font-size: 12.0px;">
<div>
<div><br/>
Hello Rory,</div>
<div>Hello Simon,</div>
<div> </div>
<div>thank you both for the feedback.</div>
<div> </div>
<div>Sorry for the delay. I'm now back on this topic after my crazy last \
week.</div>
<div> </div>
<div><br/>
First some additional explanation:<br/>
==================================<br/>
I use the LineString also to draw the track on a map (for each vessel / MMSI).<br/>
The track LineString from 00:00 of the current date to NOW() is generated \
"on-the-fly" from the stored GPS points (basically the RAW data \
you you get from the "AIS position report").<br/> For today's \
data + 1-2 days older, the "on-the-fly" generated information is combined \
with a UNION from the information retrieved from \
table 'ais_trajectory_mmsi_daily'.<br/> Currently I've got ~2 \
million new GPS position entries each day.</div>
<div> </div>
<div>In general my aggregation of points per day produces lines that go \
"over land" if there is no data for some time until the next Antenna picks \
up something again (or jumping points).<br/> As you Rory pointed out I might \
have to go the Multi-LineString way so if I don't get new GPS points for \
about ~10+ Min (vessel out of antenna range; by AIS definition the max. interval \
is 180 seconds) I start a new LineString (plus some other checks like \
distance).<br/> More on this further down below as this is an extra task.</div>
<div> </div>
<div><br/>
To answer both of your questions:<br/>
=================================<br/>
Yes I can also use MultiPoint data.<br/>
I have now added a new column to store the MultiPoint (type: PointM) to the table \
'ais_trajectory_mmsi_daily'.<br/> To test it I updated the missing data from \
my stored individual point data "AIS position report" (which is only a \
short time storage due to the large amount of data coming in).<br/> I also fixed the \
old LineString column definition which was missing the explicit SRID before.</div>
<div> </div>
<div>--<br/>
-- Add extra column to also store the PointM (Point + Timestamp) data as MultiPointM \
geometry (WGS84/GPS).<br/>
--<br/>
ALTER TABLE ais_trajectory_mmsi_daily ADD COLUMN geom_multipointm \
GEOMETRY(MULTIPOINTM, 4326);</div>
<div> </div>
<div>--<br/>
-- Change column definition to geom type of LineStringM with WGS84 (SRID: 4326).<br/>
--<br/>
ALTER TABLE ais_trajectory_mmsi_daily ALTER COLUMN geo_lat_long_line TYPE \
GEOMETRY(LINESTRINGM, 4326) USING ST_SetSRID(geo_lat_long_line, 4326);</div>
<div> </div>
<div> </div>
<div>The complete table looks like this now:<br/>
=======================================<br/>
CREATE TABLE IF NOT EXISTS aisstream.ais_trajectory_mmsi_daily<br/>
(<br/>
mmsi \
integer,<br/> min_created_on timestamp \
with time zone,<br/> \
max_created_on timestamp with time \
zone,<br/> geom_linestring geometry(LineStringM,4326), \
-- Fixed: Added exact type + SRID 4326!<br/> \
geom_multipointm geometry(MultiPointM,4326)<br/> \
)</div>
<div> </div>
<div> </div>
<div>Index (so far no index on the geometry columns):</div>
<div>================================================<br/>
CREATE INDEX IF NOT EXISTS ais_trajectory_mmsi_daily_mmsi ON \
aisstream.ais_trajectory_mmsi_daily USING btree (mmsi ASC NULLS LAST);</div>
<div> </div>
<div> </div>
<div>Additional indexes didn't really speed things up so far.<br/>
I've checked an GIST index on the 'geom_linestring' used in the \
ST_Intersects query as well as the 'min_created_on' column. I did an ANALYZE \
on the table afterwards to update the statistics. ;)</div>
<div><br/>
So far I still use the "... WHERE ST_Intersects(geom_linestring, geom_polygon) \
AND ..." with the LineString on my Polygon which is so far fast enough (0.5 \
seconds). Once I run "ST_Dump(ST_Intersection(geom_linestring, geom_polygon)) in \
the SELECT part the query times completely explode (30+ minutes). With \
MultiPoints the ST_DUMP is much faster.<br/> With the suggested MultiPoint approach \
(the newly added column) the query times are fine if I run \
"ST_Intersection(geom_multipointm, geom_polygon) on the already \
'pre-filtered' data.</div>
<div>So, I now have my points inside my polygon still in ~0.5 seconds (if I check for \
a range of ~8 days).</div>
<div> </div>
<div><br/>
Here is my current query:<br/>
=====================<br/>
--<br/>
-- Find lines intersecting the given polygon (port/berth area).<br/>
-- From those get the actual points (at least 2 point) that are inside the given \
polygon (port/berth area).<br/>
-- What we want for now is the timestamp of the first and last point inside the \
polygon.<br/>
-- ==> This information is lost by calling ST_Intersection(...)<br/>
-- in the subquery (see alias: geom_intersection)!</div>
<div>--</div>
<div>-- The query also dumps some additional details (now many points do we \
have,..).<br/>
--<br/>
--EXPLAIN ANALYZE<br/>
WITH my_port AS (<br/>
SELECT ST_GeometryFromText('POLYGON((6.972692012786865 \
50.929290943993465,6.971780061721802 50.930535164787784,6.971426010131836 \
50.93042697299645,6.972337961196899 50.92918951148323,6.972692012786865 \
50.929290943993465))', 4326) AS geom<br/> )<br/>
,my_traj AS (<br/>
SELECT ais.mmsi<br/>
,ais.min_created_on AS \
mintime<br/> ,ais.max_created_on \
AS maxtime<br/> \
,ais.geom_linestring<br/> \
,ais.geom_multipointm<br/> FROM ais_trajectory_mmsi_daily AS ais<br/>
)<br/>
SELECT geom_intersection.mmsi<br/>
,geom_intersection.created_on_date<br/>
,ST_NPoints(geom_intersection.port_mpoint_intersect) AS \
point_count<br/> \
,ST_GeometryN(geom_intersection.port_mpoint_intersect, 1) AS first_pointm<br/> \
--<br/> -- The following does not exist any \
more<br/> --<br/>
,ST_M(ST_GeometryN(geom_intersection.port_mpoint_intersect, 1)) \
AS first_pointm_m_value<br/> --<br/>
\
,ST_GeometryType(ST_GeometryN(geom_intersection.port_mpoint_intersect, 1)) AS \
first_pointm_geom_type<br/> \
,geom_intersection.port_mpoint_intersect<br/> \
,ST_GeometryType(geom_intersection.port_mpoint_intersect) AS geom_type<br/> FROM \
(<br/> SELECT my_traj.mmsi<br/>
,CAST(my_traj.mintime AS date) AS \
created_on_date<br/> \
,ST_Intersection(geom_multipointm, my_port.geom) AS \
port_mpoint_intersect<br/> \
--,ST_3DIntersection(geom_multipointm, my_port.geom) AS port_mpoint_3dintersect \
-- function does not exist error<br/> \
,my_traj.geom_linestring<br/> \
--,(ST_Dump(ST_Intersection(my_traj.geom_linestring, my_port.geom) )).geom \
AS geom_dump -- extremely slow (linestring)<br/> \
--,(ST_Dump(ST_Intersection(my_traj.geom_multipointm, \
my_port.geom) )).geom AS geom_dump -- fast \
(multipoint)<br/> FROM my_port, my_traj<br/>
WHERE ST_Intersects(my_traj.geom_linestring, my_port.geom)<br/>
AND my_traj.mintime >= '2022-10-16'<br/>
AND my_traj.maxtime < '2022-10-24'<br/>
GROUP BY mmsi, created_on_date , my_port.geom, my_traj.geom_multipointm \
,my_traj.geom_linestring<br/> ) AS geom_intersection<br/>
WHERE ST_NPoints(geom_intersection.port_mpoint_intersect) > 1<br/>
GROUP BY mmsi, created_on_date, port_mpoint_intersect<br/>
;</div>
<div> </div>
<div><br/>
1 - Issue - How do I get my timestamp for the points inside my polygon?:<br/>
========================================================================<br/>
I now have the issue that I lost the timestamp stored within my Point ("M" \
value) because ST_Intersection(..) removes this (<a \
href="https://postgis.net/docs/ST_Intersection.html" \
target="_blank">https://postgis.net/docs/ST_Intersection.html</a>).</div>
<div><br/>
The documentation kind of suggests 'ST_3DIntersection' (at least the Z is not \
removed...so I assume "M" would be available too after calling this \
function). This function should be available with PostGIS 2.1.0+.<br/> I'm \
running PostgreSQL 13.7 with PostGIS 3.2.1 here (Ubuntu Linux).<br/> The \
problem is that I get the error "function st_3dintersection(geometry, geometry) \
does not exist". I don't get why this is the case.</div>
<div> </div>
<div>So how can I get back the timestamp for my identified points?<br/>
Some kind of LATERAL JOIN on my original "geom_multipointm" column \
data?<br/> Sorry, I'm lost here.</div>
<div> </div>
<div><br/>
2 - How to creating a Multi-LineString instead of a single LineString?<br/>
I assume some kind of window function LAG/LEAD between the points.<br/>
=================================================================<br/>
The other (new) task is break down my single LineString to \
Multi-LineString while aggregating the geom point data.<br/> Currently the data \
is aggreated per day and MMSI. The 'Point' is converted to a 'PointM' \
(with the created_on timestamp) to be added \
to ST_MakeLine(geom_pointm ORDER BY created_on).</div>
<div>It seems I need some help on how to create this \
"multi-linestring".<br/> I assume I need to use a window function to \
compare the current with the next and if >= 10 Min this following point is in the \
next linestring (minimum of two points required in each group) + check distance (to \
remove total nonsense).<br/> For starters the >= 10 Min gap to generate \
Multi-LineStrings would be ok. The rest is for additional fine tuning later on. I \
started experimenting with the duration between two points using a window \
function.</div>
<div> </div>
<div><br/>
Many thanks in advance!<br/>
</div>
</div>
<div>
<div> </div>
<div>Juergen</div>
<div>
<div style="margin: 10.0px 5.0px 5.0px 10.0px;padding: 10.0px 0 10.0px \
10.0px;border-left: 2.0px solid rgb(195,217,229);"> <div style="margin: 0 0 10.0px \
0;"><b>Gesendet:</b> Montag, 17. Oktober 2022 um 15:09 Uhr<br/> \
<b>Von:</b> "Rory Meyer" <rory.meyer@VLIZ.be><br/> \
<b>An:</b> "PostGIS Users Discussion" \
<postgis-users@lists.osgeo.org><br/> <b>Betreff:</b> Re: [postgis-users] \
Issue with ST_Intersection seems related to lines crossing themself</div>
<div>
<div>
<div class="x_elementToProof" style="font-family: Calibri , Arial , Helvetica , \
sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);">Hi Juergen,</div>
<div class="x_elementToProof" style="font-family: Calibri , Arial , Helvetica , \
sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);"> </div>
<div class="x_elementToProof x_ContentPasted0" style="font-family: Calibri , Arial , \
Helvetica , sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);">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: <a href="http://openais.xyz/" id="LPNoLPOWALinkPreview" \
target="_blank">http://openais.xyz/</a> It's still a little immature and the \
database portion of the project is still closed while I rewrite some of the \
containers.</div>
<div class="x_elementToProof x_ContentPasted0" style="font-family: Calibri , Arial , \
Helvetica , sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);"> </div>
<div class="x_elementToProof x_ContentPasted0" style="font-family: Calibri , Arial , \
Helvetica , sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);">Somethings that have helped me: </div>
<div class="x_elementToProof x_ContentPasted0" style="font-family: Calibri , Arial , \
Helvetica , sans-serif;font-size: 12.0pt;color: rgb(0,0,0);background-color: \
rgb(255,255,255);"> <ul>
<li>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. </li> <li>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. </li> \
<li>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. </li> <li>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: <ul style="list-style-type: circle;">
<li>Some kind hysteresis function by having a combination of the port boundaries \
and a buffer of the port boundaries. </li> <li>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 </li> <li \
class="x_elementToProof">Ignoring entrance/exit signals where the reverse signal \
occurred within X seconds.</li> <li class="x_elementToProof">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) </li> <li class="x_elementToProof">It might be worth dumping the Z-value \
of the trajectory to get the first/last time in the intersection segment. </li> \
</ul> </li>
<li class="x_elementToProof">To make your SQL a little easier to read you can use \
Common Table Expressions:</li> </ul>
<div>```</div>
<div class="x_ContentPasted1">WITH my_port AS
<div class="x_ContentPasted1">(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 ),</div>
<div> </div>
<div class="x_ContentPasted1">my_traj AS</div>
<div class="x_ContentPasted1">(SELECT</div>
<div class="x_ContentPasted1">      ais.mmsi,</div>
<div class="x_ContentPasted1">      ais.created_on_date \
as date,</div>
<div class="x_ContentPasted1">      ais.min_created_on \
as mintime,</div>
<div class="x_ContentPasted1">      ais.max_created_on \
as maxtime,</div>
<div class="x_ContentPasted1">      ais.geo_lat_long_line \
as geom --renaming these to just help me out here</div>
<div class="x_ContentPasted1"> FROM ais_trajectory_mmsi_daily AS ais</div>
<div class="x_ContentPasted1">)</div>
<div class="x_ContentPasted1">SELECT</div>
<div class="x_ContentPasted1">      my_traj.mmsi,</div>
<div class="x_ContentPasted1">      my_traj.date,</div>
<div class="x_ContentPasted1">      ST_Intersection(my_traj.geom, \
my_port.geom) as port_traj_intersect,</div>
<div class="x_ContentPasted1">      ST_ZMin(ST_Intersection(my_traj.geom, \
my_port.geom)) as segment_start_time,</div>
<div class="x_ContentPasted1">      ST_ZMax(ST_Intersection(my_traj.geom, \
my_port.geom)) as segment_end_time</div>
<div class="x_ContentPasted1">FROM my_geom, my_traj</div>
<div class="x_ContentPasted1">WHERE ST_Intersects(my_traj.geom, my_port.geom)</div>
<div class="x_ContentPasted1">AND my_traj.mintime > '2022-10-01'</div>
AND my_traj.maxtime < '2022-10-02'</div>
<div> </div>
<div>```</div>
<div>Hope that helps some,</div>
<div>Rory</div>
</div>
<div class="x__Entity x__EType_OWALinkPreview x__EId_OWALinkPreview \
x__EReadonly_1"> </div>
<div id="x_appendonsend"> </div>
<hr style="display: inline-block;width: 98.0%;"/>
<div id="x_divRplyFwdMsg"><font color="#000000" face="Calibri, sans-serif" \
style="font-size: 11.0pt;"><b>From:</b> postgis-users \
<postgis-users-bounces@lists.osgeo.org> on behalf of Simon SPDBA Greener \
<simon@spdba.com.au><br/> <b>Sent:</b> 17 October 2022 01:56<br/>
<b>To:</b> postgis-users@lists.osgeo.org <postgis-users@lists.osgeo.org><br/>
<b>Subject:</b> Re: [postgis-users] Issue with ST_Intersection seems related to lines \
crossing themself</font>
<div> </div>
</div>
<div>
<p>Would using a MultiPoint (or GeometryCollection with Points) for the GPS locations \
rather than a Linestring be possible?</p>
<p>Simon</p>
<div class="x_x_moz-cite-prefix">On 17/10/2022 10:10 am, <a \
class="x_x_moz-txt-link-abbreviated" href="mailto:magog002@web.de" \
onclick="parent.window.location.href='mailto:magog002@web.de'; return false;" \
target="_blank"> magog002@web.de</a> wrote:</div>
<blockquote>
<div style="font-family: Verdana;font-size: 12.0px;">
<div>Hi,</div>
<div> </div>
<div>I've got an issue with an SQL PostGIS query and hope to get a solution \
here.</div>
<div> </div>
<div>First some details.</div>
<div>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.</div>
<div>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.</div>
<div> </div>
<div>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.</div>
<div>So far so good.</div>
<div>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.</div>
<div> </div>
<div> </div>
<div>Here the basic search query that basically works:</div>
<div>================================</div>
<div>--<br/>
-- Return all vessels (MMSI) found in the selected timeframe where the<br/>
-- LineString Intersects the Polygon given as parameter to search.<br/>
--<br/>
SELECT ais_trajectory_mmsi_daily.mmsi<br/>
,ais_trajectory_mmsi_daily.min_created_on<br/>
,ais_trajectory_mmsi_daily.max_created_on<br/>
,ST_SetSRID(ais_trajectory_mmsi_daily.geo_lat_long_line, 4326) \
AS geo_lat_long_line<br/> FROM ais_trajectory_mmsi_daily<br/>
WHERE ST_Intersects(ST_SetSRID(ais_trajectory_mmsi_daily.geo_lat_long_line, \
4326),<br/> \
ST_GeometryFromText('POLYGON((7.1084729427821 50.7352828020046,7.1089771980769 \
50.7353778664669,7.10848259899649 50.7365542732199,7.10801482200623 \
50.7364673592064,7.1084729427821 50.7352828020046))', 4326) )<br/> AND \
min_created_on >= CAST('2022-10-01' AS TIMESTAMPTZ)<br/> AND \
max_created_on < CAST('2022-10-15' AS TIMESTAMPTZ)<br/> ;</div>
<div> </div>
<div> </div>
<div>Here is now the query I use with ST_Intersection and where it starts to go \
wrong.</div>
<div>When a vessel stays a the same place the GPS coordinates are the same + \
drift.</div>
<div> </div>
<div>As those points are all added to the LineString the line sections are running \
over themself and it looks like a "woolen noodle".</div>
<div>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).</div>
<div> </div>
<div>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.</div>
<div> </div>
<div>--<br/>
-- Outer query with union to visualize the Polygon and the <br/>
-- data received from the LineString that Intersects<br/>
-- with the polygon via e.g. pgAdmin.<br/>
--<br/>
SELECT NULL AS mmsi<br/>
,NULL AS created_on_date<br/>
,NULL AS min_created_on<br/>
,NULL AS max_created_on<br/>
,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<br/> ,NULL AS \
geom_start<br/> ,NULL AS geom_end<br/>
UNION<br/>
-- The actual data of interest<br/>
SELECT intersection.mmsi<br/>
,intersection.created_on_date<br/>
,intersection.min_created_on<br/>
,intersection.max_created_on<br/>
,intersection.geom<br/>
,MIN(DATE_TRUNC('second', \
TO_TIMESTAMP(ST_InterpolatePoint(intersection.geo_lat_long_line, \
ST_StartPoint(intersection.geom)) )::timestamptz )) AS geom_start<br/> \
,MAX(DATE_TRUNC('second', \
TO_TIMESTAMP(ST_InterpolatePoint(intersection.geo_lat_long_line, \
ST_EndPoint(intersection.geom)) )::timestamptz )) AS geom_end<br/> FROM (<br/>
SELECT ais_trajectory_mmsi_daily.mmsi<br/>
,CAST(ais_trajectory_mmsi_daily.min_created_on AS \
date) AS created_on_date <br/> \
,ais_trajectory_mmsi_daily.min_created_on<br/> \
,ais_trajectory_mmsi_daily.max_created_on<br/> \
,ST_SetSRID(ais_trajectory_mmsi_daily.geo_lat_long_line, 4326) AS \
geo_lat_long_line</div>
<div> --</div>
<div> -- Here seems to be the problem...</div>
<div> --<br/>
\
,(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)<br/> \
\
,ST_SetSRID(ais_trajectory_mmsi_daily.geo_lat_long_line, 4326)<br/> \
)<br/> \
)).geom AS geom<br/> \
FROM ais_trajectory_mmsi_daily<br/> WHERE \
ST_Intersects(ST_SetSRID(ais_trajectory_mmsi_daily.geo_lat_long_line, 4326),<br/> \
\
ST_GeometryFromText('POLYGON((7.1084729427821 \
50.7352828020046,7.1089771980769 50.7353778664669,7.10848259899649 \
50.7365542732199,7.10801482200623 50.7364673592064,7.1084729427821 \
50.7352828020046))', 4326) )<br/> AND min_created_on \
>= CAST('2022-10-01' AS TIMESTAMPTZ)<br/> AND \
max_created_on < CAST('2022-10-02' AS TIMESTAMPTZ)<br/> ) AS \
intersection<br/> GROUP BY intersection.mmsi<br/>
,intersection.created_on_date<br/>
,intersection.min_created_on<br/>
,intersection.max_created_on<br/>
,intersection.geom<br/>
ORDER BY created_on_date, mmsi, geom_start<br/>
;</div>
<div> </div>
<div>If needed I can try to create some sample with data to put \
on sqlfiddle.</div>
<div> </div>
<div>On Stackoverflow I found this which seems to be related to the same issue:</div>
<div><a class="x_x_moz-txt-link-freetext" \
href="https://stackoverflow.com/questions/62090829/why-st-intersection-from-postgis-returns-multilinestring-for-self-intersecting-l" \
target="_blank">https://stackoverflow.com/questions/62090829/why-st-intersection-from-postgis-returns-multilinestring-for-self-intersecting-l</a></div>
<div> </div>
<div>Many thanks in advance for hints to solve the problem.</div>
<div> </div>
<div>Juergen</div>
<div> </div>
</div>
<fieldset class="x_x_moz-mime-attachment-header"> </fieldset>
<pre class="x_x_moz-quote-pre">_______________________________________________
postgis-users mailing list
<a class="x_x_moz-txt-link-abbreviated" href="mailto:postgis-users@lists.osgeo.org" \
onclick="parent.window.location.href='mailto:postgis-users@lists.osgeo.org'; \
return false;" target="_blank">postgis-users@lists.osgeo.org</a> <a \
class="x_x_moz-txt-link-freetext" \
href="https://lists.osgeo.org/mailman/listinfo/postgis-users" \
target="_blank">https://lists.osgeo.org/mailman/listinfo/postgis-users</a> </pre>
</blockquote>
<pre class="x_x_moz-signature">--
Simon Greener
39 Cliff View Drive, Allens Rivulet, 7150, Tasmania, Australia
(m) +61 418 396 391
(w) <a class="x_x_moz-txt-link-abbreviated" href="http://www.spdba.com.au" \
target="_blank">www.spdba.com.au</a> (m) <a class="x_x_moz-txt-link-abbreviated" \
href="mailto:simon@spdba.com.au" \
onclick="parent.window.location.href='mailto:simon@spdba.com.au'; return \
false;" target="_blank">simon@spdba.com.au</a></pre> </div>
_______________________________________________ postgis-users mailing list \
postgis-users@lists.osgeo.org <a \
href="https://lists.osgeo.org/mailman/listinfo/postgis-users" target="_blank"> \
https://lists.osgeo.org/mailman/listinfo/postgis-users</a></div> </div>
</div>
</div>
</div>
</div>
</div>
_______________________________________________ postgis-users mailing list \
postgis-users@lists.osgeo.org <a \
href="https://lists.osgeo.org/mailman/listinfo/postgis-users" \
target="_blank">https://lists.osgeo.org/mailman/listinfo/postgis-users</a></div> \
</div> </div>
</div>
</div></div></body></html>
_______________________________________________
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic