Re: Spatial join issues
Paul Ramsey via postgis-users <[email protected]> Wed, 29 Oct 2025 10:25:08 -0700
| Newsgroups | gmane.comp.gis.postgis |
|---|---|
| Message-ID | <CACowWR3MDzQcsMLKpa4CdB6mOONggS6K3mj7xV=Eywwr_9mmbQ@mail.gmail.com> |
--000000000000ae594206424f6acc
Content-Type: text/plain; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable
I would imagine that my little example is a bit error prone, but the
classic solution to "do these differently structured things with slightly
different coordinates cover basically the same space" is to inspect the
ratio "area(intersection(a, b)) / area(union(a, b))". You can see
intuitively how the more the same two polygons are, the closer that ratio
will be to 1.0.
WITH p AS (
SELECT 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'::geometry AS p1,
'POLYGON((0 0, 10 0, 10 5, 10 10, 0 10, 0 0))'::geometry AS p2,
0.001 AS shift
),
shifted AS (
SELECT p1, ST_Translate(p2, shift, shift) AS p2, shift
FROM p
),
areas AS (
SELECT ST_Area(ST_Union(p1, p2, shift/10)) AS area_union,
ST_Area(ST_Intersection(p1, p2, shift/10)) AS area_inter
FROM shifted
)
SELECT ST_AsText(shifted.p1) AS p1_orig,
ST_AsText(shifted.p2) AS p2_orig,
area_inter/area_union AS ratio
FROM areas, shifted;
On Wed, Oct 29, 2025 at 10:17=E2=80=AFAM Paul Ramsey <pramsey@cleverelephan=
t.ca>
wrote:
> Greg isn't passing judgement on how hard your problem is, he's saying you
> haven't explained it particularly well. Pictures help. Taking a guess at
> what you mean, here's some SQL that creates two polygons, with slightly
> different structure, and slightly different coordinates, that describe th=
e
> same general space in the universe, and then massages them until they pas=
s
> an equals test.
>
> WITH p AS (
> SELECT 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'::geometry AS p1,
> 'POLYGON((0 0, 10 0, 10 5, 10 10, 0 10, 0 0))'::geometry AS p2
> ),
> shifted AS (
> SELECT p1, ST_Translate(p2, 0.0001, 0.0001) AS p2
> FROM p
> ),
> rp AS (
> SELECT ST_ReducePrecision(p1,0.1) AS p1,
> ST_ReducePrecision(p2,0.1) AS p2
> FROM shifted
> ),
> snap AS (
> SELECT ST_Snap(p1,p2,0.1) AS p1,
> ST_Snap(p2,p1,0.1) AS p2
> FROM rp
> )
> SELECT ST_AsText(shifted.p1) AS p1_orig,
> ST_AsText(shifted.p2) AS p2_orig,
> ST_AsText(snap.p1) AS p1_snap,
> ST_AsText(snap.p2) AS p2_snap,
> ST_Equals(snap.p1, snap.p2)
> FROM snap, shifted;
>
> -[ RECORD 1
> ]------------------------------------------------------------------------=
--------------------------
> p1_orig | POLYGON((0 0,10 0,10 10,0 10,0 0))
> p2_orig | POLYGON((0.0001 0.0001,10.0001 0.0001,10.0001 5.0001,10.0001
> 10.0001,0.0001 10.0001,0.0001 0.0001))
> p1_snap | POLYGON((0 10,10 10,10 5,10 0,0 0,0 10))
> p2_snap | POLYGON((0 10,10 10,10 5,10 0,0 0,0 10))
> st_equals | t
>
> On Wed, Oct 29, 2025 at 10:00=E2=80=AFAM Shaozhong SHI <shishaozhong@gmai=
l.com>
> wrote:
>
>> This is very challenging. Take my words for it. Try on any polygons yo=
u
>> created and modified.
>>
>> On Wed, 29 Oct 2025 at 14:36, Greg Troxel <[email protected]> wrote:
>>
>>> Shaozhong SHI <[email protected]> writes:
>>>
>>> > Visually, there appears some matching polygons. Even if two
>>> geometries
>>> > represent the same shape visually, they might not be considered equal
>>> due
>>> > to tiny differences in precision or metadata. Have you encountered
>>> > problems of failure of spatial join? How did you overcome the
>>> problems?
>>> > Regards, David
>>>
>>> Could you post your example polygons, and the queries you are using?
>>> Your question is much too open ended. It even sounds like it might be =
a
>>> request for help with GIS homework, but it's hard to tell :-)
>>>
>>>
--000000000000ae594206424f6acc
Content-Type: text/html; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable
<div dir=3D"ltr">I would imagine that my little example is a bit error pron=
e, but the classic solution to "do these differently structured things=
with slightly different coordinates cover basically the same space" i=
s to inspect the ratio =C2=A0"area(intersection(a, b)) / area(union(a,=
b))". You can see intuitively how the more the same two polygons are,=
the closer that ratio will be to 1.0.=C2=A0<div><br></div><div>WITH p AS (=
<br>=C2=A0 SELECT 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'::geometr=
y AS p1,<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0'POLYGON((0 0, 10 0, 10 5=
, 10 10, 0 10, 0 0))'::geometry AS p2,<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 =
=C2=A00.001 AS shift<br>),<br>shifted AS (<br>=C2=A0 SELECT p1, ST_Translat=
e(p2, shift, shift) AS p2, shift<br>=C2=A0 FROM p<br>),<br>areas AS (<br>=
=C2=A0 SELECT ST_Area(ST_Union(p1, p2, shift/10)) AS area_union,<br>=C2=A0 =
=C2=A0 =C2=A0 =C2=A0 =C2=A0ST_Area(ST_Intersection(p1, p2, shift/10)) AS ar=
ea_inter<br>=C2=A0 FROM shifted<br>)<br>SELECT ST_AsText(shifted.p1) AS p1_=
orig, <br>=C2=A0 =C2=A0 =C2=A0 =C2=A0ST_AsText(shifted.p2) AS p2_orig, <br>=
=C2=A0 =C2=A0 =C2=A0 =C2=A0area_inter/area_union AS ratio<br>FROM areas, sh=
ifted;<br></div></div><br><div class=3D"gmail_quote gmail_quote_container">=
<div dir=3D"ltr" class=3D"gmail_attr">On Wed, Oct 29, 2025 at 10:17=E2=80=
=AFAM Paul Ramsey <<a href=3D"mailto:[email protected]">pramsey@=
cleverelephant.ca</a>> wrote:<br></div><blockquote class=3D"gmail_quote"=
style=3D"margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:=
solid;border-left-color:rgb(204,204,204);padding-left:1ex"><div dir=3D"ltr"=
>Greg isn't passing judgement on how hard your problem is, he's say=
ing you haven't explained it particularly well. Pictures help. Taking a=
guess at what you mean, here's some SQL that creates two polygons, wit=
h slightly different structure, and slightly different coordinates, that de=
scribe the same general space in the universe, and then massages them until=
they pass an equals test.<div><br></div><div>WITH p AS (<br>=C2=A0 SELECT =
'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'::geometry AS p1,<br>=C2=A0=
=C2=A0 =C2=A0 =C2=A0 =C2=A0'POLYGON((0 0, 10 0, 10 5, 10 10, 0 10, 0 0=
))'::geometry AS p2<br>),<br>shifted AS (<br>=C2=A0 SELECT p1, ST_Trans=
late(p2, 0.0001, 0.0001) AS p2<br>=C2=A0 FROM p<br>),<br>rp AS (<br>=C2=A0 =
SELECT ST_ReducePrecision(p1,0.1) AS p1,<br>=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=
=A0ST_ReducePrecision(p2,0.1) AS p2<br>=C2=A0 FROM shifted<br>),<br>snap AS=
(<br>=C2=A0 SELECT ST_Snap(p1,p2,0.1) AS p1, <br>=C2=A0 =C2=A0 =C2=A0 =C2=
=A0 =C2=A0ST_Snap(p2,p1,0.1) AS p2<br>=C2=A0 FROM rp<br>)<br>SELECT ST_AsTe=
xt(shifted.p1) AS p1_orig, <br>=C2=A0 =C2=A0 =C2=A0 =C2=A0ST_AsText(shifted=
.p2) AS p2_orig, <br>=C2=A0 =C2=A0 =C2=A0 =C2=A0ST_AsText(snap.p1) AS p1_sn=
ap, <br>=C2=A0 =C2=A0 =C2=A0 =C2=A0ST_AsText(snap.p2) AS p2_snap,<br>=C2=A0=
=C2=A0 =C2=A0 =C2=A0ST_Equals(snap.p1, snap.p2)<br>FROM snap, shifted;<br>=
</div><div><br></div><div>-[ RECORD 1 ]------------------------------------=
--------------------------------------------------------------<br>p1_orig =
=C2=A0 | POLYGON((0 0,10 0,10 10,0 10,0 0))<br>p2_orig =C2=A0 | POLYGON((0.=
0001 0.0001,10.0001 0.0001,10.0001 5.0001,10.0001 10.0001,0.0001 10.0001,0.=
0001 0.0001))<br>p1_snap =C2=A0 | POLYGON((0 10,10 10,10 5,10 0,0 0,0 10))<=
br>p2_snap =C2=A0 | POLYGON((0 10,10 10,10 5,10 0,0 0,0 10))<br>st_equals |=
t<br></div></div><br><div class=3D"gmail_quote"><div dir=3D"ltr" class=3D"=
gmail_attr">On Wed, Oct 29, 2025 at 10:00=E2=80=AFAM Shaozhong SHI <<a h=
ref=3D"mailto:[email protected]" target=3D"_blank">shishaozhong@gmail.=
com</a>> wrote:<br></div><blockquote class=3D"gmail_quote" style=3D"marg=
in:0px 0px 0px 0.8ex;border-left-width:1px;border-left-style:solid;border-l=
eft-color:rgb(204,204,204);padding-left:1ex"><div dir=3D"ltr">This is very =
challenging.=C2=A0 Take my words for it.=C2=A0 Try on any polygons you crea=
ted and modified.</div><br><div class=3D"gmail_quote"><div dir=3D"ltr" clas=
s=3D"gmail_attr">On Wed, 29 Oct 2025 at 14:36, Greg Troxel <<a href=3D"m=
ailto:[email protected]" target=3D"_blank">[email protected]</a>> wrote:<br></=
div><blockquote class=3D"gmail_quote" style=3D"margin:0px 0px 0px 0.8ex;bor=
der-left-width:1px;border-left-style:solid;border-left-color:rgb(204,204,20=
4);padding-left:1ex">Shaozhong SHI <<a href=3D"mailto:shishaozhong@gmail=
.com" target=3D"_blank">[email protected]</a>> writes:<br>
<br>
>=C2=A0 =C2=A0Visually, there appears some matching polygons.=C2=A0 Even=
if two geometries<br>
> represent the same shape visually, they might not be considered equal =
due<br>
> to tiny differences in precision or metadata.=C2=A0 Have you encounter=
ed<br>
> problems of failure of spatial join?=C2=A0 How did you overcome the pr=
oblems?<br>
> Regards, David<br>
<br>
Could you post your example polygons, and the queries you are using?<br>
Your question is much too open ended.=C2=A0 It even sounds like it might be=
a<br>
request for help with GIS homework, but it's hard to tell :-)<br>
<br>
</blockquote></div>
</blockquote></div>
</blockquote></div>
--000000000000ae594206424f6acc--