Why does ST_Area(g.geom) + ST_Area(p.geom) - ST_Area(ST_Union(g.geom,p.geom)) produce negative result when using ST_Overlaps(g.geom, p.geom) = true?

Question

I have two set of shape data. Let's say file gov.shp contains list of post codes and their boundaries. File post.shp contains list of government area name and their boundaries. I want to find out the the union (may not be the correct word here) between two files (i.e. post code 1, 2 and 3 overlaps with government area a; post code 1 also overlaps with government area b). Bear in mind that government area is bigger than post code area. So I have a query like below:

select  g.id, g.name, p.id, p.code
from    gov g
inner join
        poa p
on      ST_Overlaps(g.geom, p.geom) = true

However if I apply this where clause, to my surprise the count(*) actually returns results. I am confused by why total area of g + p - union_area of(g,p) can be negative...

where   ST_Area(g.geom) + ST_Area(p.geom) - ST_Area(ST_Union(g.geom, p.geom)) < 0

Am I using the wrong function (ST_Overlaps)? Should I be using ST_Intersecs/ST_Equals/other functions?

I am new to this area so I am not too sure what the differences are in terms of solving my particular problem.

I imported shape file as geometry and MULTIPOLYGON object and I tried different results in both PostGIS and SQL Server 2008 R2 and both produced slightly different results. The above sql query syntax is in PostGIS though.

--SQL Server 2008 -- STOverlaps -- DIFF < 0 = 27 -- DIFF > 0 = 1646

--SQL Server 2008 -- STEquals() -- 42

--PostGIS -- DIFF < 0 = 730 -- DIFF > 0 = 5312

--PostGIS -- STEquals() -- 0

SQL Server 2008 syntax

select  g.id, g.name, p.id, p.code
from    gov g
inner join
        poa p
on    g.geom.MakeValid().STOverlaps(p.geom.MakeValid()) = 1
--for some reason I had to use *.geom.MakeValid()
where g.geom.MakeValid().STArea() + p.geom.MakeValid().STArea() - g.geom.MakeValid().STUnion(p.geom.MakeValid()).STArea() < 0

I used Shape2SQL to import into SQL Server 2008, Planer Geometry, didn't check SRID, create spatial index.

And I used PostGIS Shapefile Import/Export Manager which came with PostGIS 2.0 to import into PostGIS, and all default settings.

Answer 1

In PostGIS terms, ST_Overlaps(g1, g2) means that the two geometries intersect but neither totally contains the other. If, as you mentioned, many of the poa polygons are totally contained within the gov polygons, they will not be included in the ST_Overlaps() function. ST_Union() on the other hand, just combines all areas, removing intersections. So the area of a union in the total area of both geometries minus intersections. This could certainly be greater than the overlap.

A simple example: if gov contains one large polygon, and poa contains one smaller polygon, totally enclosed in gov, then Overlaps(gov, pos) is false, so it's area would not be included in your query, but ST_Area(ST_Union()) would return the area of gov which is > 0, thus the negative result.

Answer 2

If you choose to work with PostGIS, then you'll need to get the two shapes into a PostGIS database, and create spatial indices for both. The shp2pgsql tool can do that:

shp2pgsql -I -W <encoding> -s <EPSG code> gov.shp | psql -d <your database> 
shp2pgsql -I -W <encoding> -s <EPSG code> post.shp | psql -d <your database>

Now connect to your DB and have a look at the following queries: (I'm assuming the government areas have a 'name' column and the post codes have a 'zip' column)

SELECT g.name AS 'Govt Area', p.zip AS 'Zip code' 
FROM gov AS g, post as p
WHERE ST_Intersects(g.the_geom, p.the_geom)
GROUP BY g.name, p.zip;

THis should return a list of all the government areas, and each zip code that intersects that area. For each gov area there will be multiple rows for each zip that intersects that area. Also, post codes will also appear multiple times when a post code extends across to a second gov area.

You might also consider:

SELECT g.name AS 'Govt Area', p.zip AS 'Zip code' 
FROM gov AS g, post as p
WHERE ST_Within(g.the_geom, p.the_geom)
GROUP BY g.name, p.zip;

This will list only those post codes that are totally within one gov area. Those post codes that overlap more than one gov area will not appear at all. Furthermore, ST_Crosses() will find only those post codes that intersect at least two gov areas.

Are we getting any closer?