[GIS] Comparing two tables, using PostgreSQL Intersect (on PostGIS geometry)

postgis

I am trying to compare two tables, that are suppose to be equal. One is a reference table (with the "correct" values) and the other is a test output. That is, they have same columns and columntypes. I then want to find the rows that are "identical".
I figured that PostgreSQL's' Intersect (or maybe its opposite, Except) would work here, since it selects rows that appear in both tables.

This works fine if I give it rows that are completely the same:

select 2, ST_Polygon(st_geomfromtext('LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)'), 4326) 
intersect
select 2, ST_Polygon(st_geomfromtext('LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)'), 4326)

returns:

row 1: 2;0103000020E6100000010000000500000000000000000000000000000000000000000000000000F03F0000000000000000000000000000F03F000000000000F03F0000000000000000000000000000F03F00000000000000000000000000000000

However, if the start/end vertex is changed in the geometry column in one of the rows, it returns zero rows:

select 2, ST_Polygon(st_geomfromtext('LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)'), 4326) 
intersect
select 2, ST_Polygon(st_geomfromtext('LINESTRING( 1 0, 1 1, 0 1, 0 0, 1 0)'), 4326)

That is, if i used ST_Equals on the geometries, it would return true. But in Intersect, this is not true.

Is there a smarter way to compare entire tables?

Best Answer

You are mixing non-spatial functionality with spatial objects. PostgreSQL's intersect does not know about spatial. Stay in PostGIS for everything spatial. In other words, use the ST_ functions (Spatial Type) when you work with spatial objects.

For intersections there are ST_Intersects for a boolean condition and ST_Intersection for the geometric intersection.

SELECT ST_Intersects(
  ST_Polygon(ST_GeomFromText('LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)'), 4326),
  ST_Polygon(ST_GeomFromText('LINESTRING(1 0, 1 1, 0 1, 0 0, 1 0)'), 4326)
);

Related Solutions

[GIS] PostGIS: ST_Equals false when ST_Intersection=100% of geometry

My guess is that you coordinate transformations have introduced tiny rounding errors (see an example below). As there is no way to set the tolerance in ST_Equals, this is causing ST_Equals to return false for some geometries that only differ in the nth decimal place, as the geometries have to be identical in every respect -- see the intersection matrix definition in libgeos. You can check this with a really extreme example,

SELECT ST_Equals(
      ST_MakePoint(0,0),
      ST_MakePoint(0,0.000000000000000000000000000000000000000000000000000000000001));

which returns false.

If you use ST_SnapToGrid you can impose a given precision, for example, to ten decimal places,

SELECT ST_Equals(
      ST_MakePoint(0,0),
      ST_SnapToGrid(
             ST_MakePoint(0,0.00000000000000000000000000000000000000000000001),
      10));

now returns true.

If you were to run,

CREATE table matchdata AS
SELECT  a.*
FROM gg2014 a, gg2013 b
WHERE ST_Equals(ST_SnapToGrid(a.g1, 5), ST_SnapToGrid(b.g1, 5));

setting an appropriate tolerance, I suspect your problems would disappear.

Here is a link to a Postgis developer discussion about tolerance which suggests it is less than trivial to implement.

I did a couple of conversions between British National Grid (EPSG:27700) and lat/lon to illustrate the point about rounding precision, taking a point somewhere in London,

SELECT ST_AsText(ST_Transform(ST_SetSrid(ST_MakePoint(525000, 190000),27700),4326));

returns POINT(-0.19680497282746 51.5949871603888)

and reversing this,

SELECT ST_AsText(ST_Transform(ST_SetSrid(ST_MakePoint(-0.19680497282746, 51.5949871603888),4326),27700));

returns POINT(525000.000880007 189999.999516211)

which is off by less than a millimeter, but more than enough to make ST_Equals return false.

[GIS] How to intersect two polygon-layers in PostgreSQL/PostGIS

I expected to get a new table(layer) with a multigeometry of 382 multipolygons, but it did not work

Your expectation is incorrect. PostGIS is only returning ANY combination of layer_1 and layer_2 that intersects. If every polygon in layer_2 straddled two polygons in layer_1 the result of

SELECT COUNT(*)
FROM layer_1, layer_2
WHERE ST_Intersects(layer_1.geom, layer_2.geom)

Would be 897 * 2.

If you want 382 multipolygons you need to use a combination of ST_Collect with GROUP BY like so

SELECT layer_1.id, ST_COLLECT(ST_INTERSECTION(layer_1.geom, layer_2.geom))
FROM layer_1, layer_2
WHERE ST_Intersects(layer_1.geom, layer_2.geom)
GROUP BY layer_1.id

after about one hour this query "crashed" my hp-laptop

Are you using a spatial index on either layer? If not then PostGIS is doing 382 * 897 comparisons to check if each layer_1 polygon intersects with each layer_2 polygon, which may take a while.

Best Answer

Related Solutions

[GIS] PostGIS: ST_Equals false when ST_Intersection=100% of geometry

[GIS] How to intersect two polygon-layers in PostgreSQL/PostGIS

Related Question