Regina,
comments inline below
Obe, Regina wrote:
On a slightly related note, I'm still confused about how differences in
PostGIS work though particularly at the boundary level.
E.g. why when you take a polygon difference out another polygon - that
your new difference intersects with both original polygons.
There's a couple of reasons this will happen. One is roundoff error -
if the difference poly contains new vertices which are the result of
intersections of line segments, then it's very likely that the new
vertex will intersect one or other parent polygon.
The other reason is that the boundaries of the difference and input
polys will definitely intersect, since they both contain some of the
same linework.
The condition which difference actually should theoretically satisfy is
that the Interior of the difference poly will not intersect the Interior
of the poly it was subtracted from. If you use inputs where
intersections can be computed exactly (e.g. rectangles) you should find
this holds.
Well rather frankly I don't see how you would make it so it wouldn't,
but seems to violate point set definition of difference to me (if you
remove an intersecting part how can the new geometry intersect with the
originals?). Seems like it should just touch but not intersect.
Note - touches() implies intersects(). As I mentioned above, you need
to be thinking in terms of testing Interior Intersection. There are no
named predicates which do this - you have to work with the DE-9IM
pattern directly. (I have added an interiorsIntersect predicate to JTS
to do just this, since it's useful in some situations.
More disturbing is why when I cut a polygon with a line I get (shall we
say bizarre results). A slightly larger polygon.
This is a weird one. Can you post the actual geometries? In JTS when I
try differencing a line from a circular point buffer, I get a single
polygon (with extra vertices where the line cut it). This is the
expected behaviour. I don't know why PostGIS/GEOS is doing something
different.
To answer a question you posed on your blog, the reason that when you
subtract a line from a polygon you get basically the same polygon,
rather than say the polygon split into two halves, is that if the latter
was provided as a MultiPolygon it would be invalid, because the halves
would share line segments down the middle. Also, if the line did not
fully overlap the polygon there's no choice - you have to return the
original poly. The behaviour is thus consistent between the two cases.
Unfortunately the SFS (and no other standard I'm aware of) doesn't
define the precise semantics of the overlay operations. So I made 'em
up! Hopefully they are consistent and reasonable. There's no doubt
alternative definitions which might be useful in some cases - but you
have to choose one definition for any given function. (For the
situation above, many people would like a "cut polygon by line"
operation - hopefully that will get provided as a new function sometime
soon).
HTH - Martin
Thanks,
Regina
-----Original Message-----
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of
Martin Davis
Sent: Tuesday, June 03, 2008 12:08 PM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] How to tell if 2 geometries are spatially
equal
Andy,
Can you send me the PDF that you found? That link appears to be dead.
I realize that the SFS doesn't actually define equal. I was the
designer of JTS and GEOS, so the definition is due to me. I followed
what I thought was a logical extension of the other definitions, and
something that was expressible in terms of the DE-9IM. I also provided
"equalsExact", which corresponds to the other proposed definition that
you give.
I will note that in fact according to the SFS, in the geometry
LINESTRING(0 0, 5 5, 10 10)' the point (5,5) is *not* a boundary point.
In fact, even in the geometry MULTILINESTRING((0 0, 5 5), (5 5, 10 10))
the point (5 5) is not on the boundary, due to the (slightly bizarre)
'mod-2' rule used by the SFS.
There's plenty of references that support the definition I chose (in
fact, I suspect I chose it based on other references I scanned at the
time):
* The IBM Spatial Datablade manual says: Using the *ST_Equals()*
function is functionally equivalent to using
*ST_IsEmpty*(*ST_SymDifference*(/a/,/b/))
(http://publib.boulder.ibm.com/infocenter/idshelp/v10/index.jsp?topic=/c
om.ibm.spatial.doc/spat122.htm)
ESRI gives the DE-9IM pattern of T*F**FFF* for ST_Equals, which is what
JTS uses
http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Spatial_re
lationships
Egenhofer has a paper on "Point-Set Topological Spatial Relations"
(http://www.spatial.maine.edu/~max/pointset.pdf)
There's also all the references given in the SFS paper.
In the end it all comes down to naming. There are various kinds of
equality, which are useful for different things. Different systems name
them differently, which is ok as long as the semantics are documented.
Of course, it's nice to have some standard names - and it looks to me
like ST_equals is pretty well defined.
Andy Anderson wrote:
On Jun 2, 2008, at 4:29 PM, Martin Davis wrote:
I use "topologically equal" because the OGC SFS specification uses
the term "topology" extensively in their discussion of the meaning of
the DE-9IM model, on which the semantics of ST_equals is based.
No doubt they do, because the spatial relationships that can be
determined by the DE-9IM are, generally speaking, topological in
nature:
Disjoint, Touches, Crosses, Within and Overlaps
These are the ones defined here:
http://portal.opengeospatial.org/files/?artifact_id=829
Oddly, in this document they also mention "Equals" but they never
explicitly define it (though they claim to later in the document).
However, I did find a very complete description of DE-9IM here:
http://mlblog.osdir.com/gis.postgis/2004-02/pdfHaVE9FZPMj.pdf
and they say that for "equals", the "Geometries must be identical:
- Same dimension
- Same geometry type
- Same number of vertices
- All x,y coordinates must be identical"
So this definition would seem to agree with me that LINESTRING(0 0, 10
10) and LINESTRING(0 0, 5 5, 10 10) are not topologically equal. The
best way to think of this is that the second is a *polyline*, and the
vertex (5, 5) is a boundary between two lines, and that changes the
DE-9IM.
As you point out, however, ST_equals('LINESTRING(0 0, 10
10)','LINESTRING(0 0, 5 5, 10 10)') returns true; perhaps this was by
design to distinguish it from 'LINESTRING(0 0, 10 10)'::geometry ~=
'LINESTRING(0 0, 5 5, 10 10)'::geometry .
If you can point to further discussion on this issue, I would be
interested.
I don't like the term "spatially-equal", because I think "spatially"
is too vague and overloaded. How about "point-set equal"? The idea
is that A = B iff every point of A is in B and every point of B is
in A.
Personally, I don't have a problem with "spatially equal", it says to
me "they fill space in the same way", which in fact they do because a
point has no extent.
-- Andy
Andy Anderson wrote:
I wouldn't call this example "topologically" equal; one has two
vertices and the other has three, and that's the only characteristic
that's relevant in topology (not even their positions :-)
"Coincident" is probably a better term, though "spatially equal" is
probably just as good, and contrasts well with the term
"geometrically equal" that the manual uses to describe the ~=
operator.
-- Andy
On May 30, 2008, at 7:54 PM, Martin Davis wrote:
Will ST_equals do what you want? It reports whether two geometries
are topologically equal.
(So for example, ST_equals('LINESTRING(0 0, 10 10)','LINESTRING(0
0, 5 5, 10 10)') is true)
Obe, Regina wrote:
I recall this having come up before. I always thought that ~=
would
tell me if 2 geometries are spatially equal but it doesn't seem
to.
The only way I can figure to determine spatial equality is if
ST_Within(A,B) And ST_Within(B,A) (or ST_Difference(A,B) AND
ST_Difference(B,A) both return an empty geometry collection)
--So case in point
SELECT geom1 ~= geom2 as what, ST_Within(geom1, geom2) AND
ST_Within(geom2, geom1) As spatial_equal,
ST_AsText(ST_Difference(geom1, geom2)) as diffgeom12,
ST_AsText(ST_Difference(geom2, geom1)) as diffgeom21 FROM (SELECT
'LINESTRING(1 1, 1 2, 1 3)'::geometry As geom1, 'LINESTRING(1
1, 1 3)'::geometry As geom2) As foo
Results:
what | spatial_equal | diffgeom12 |
diffgeom21
------+---------------+--------------------------+----------------------
----
f | t | GEOMETRYCOLLECTION EMPTY |
GEOMETRYCOLLECTION
EMPTY
Is there a function / operator that does that (also what does
geom1 =
geom2 compare - is it just bounding boxes or is that the spatially
equal
operator I am looking for?)
Thanks,
Regina
-----------------------------------------
The substance of this message, including any attachments, may be
confidential, legally privileged and/or exempt from disclosure
pursuant to Massachusetts law. It is intended
solely for the addressee. If you received this in error, please
contact the sender and delete the material from any computer.
_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users
--
Martin Davis
Senior Technical Architect
Refractions Research, Inc.
(250) 383-3022
_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users
_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users
--
Martin Davis
Senior Technical Architect
Refractions Research, Inc.
(250) 383-3022
_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users
_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users
--
Martin Davis
Senior Technical Architect
Refractions Research, Inc.
(250) 383-3022
_______________________________________________
postgis-users mailing list
[email protected]
http://postgis.refractions.net/mailman/listinfo/postgis-users
