Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
    
I'm assuming your datasets are not 100% topologically correct with eachother, are they from different sources / authorities? Do they have the same SRID? –  Geoist Jul 13 '12 at 6:15
    
@Geoist the two datasets do not need to be 100% topologically correct that's why I am doing this exercise. I think this make sense isn't it? And they both have the same SRID, because they have the same content in respective .prj file. And the content from .prj file is: GEOGCS["GCS_GDA_1994",DATUM["D_GDA_1994",SPHEROID["GRS_1980",6378137,298.2572221‌​01]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]] –  Jeff Jul 13 '12 at 6:22
    
I am not 100% sure what SRID have to do with calculation since I am not projecting the map (like seeing/showing the map). Excuse me that I am not 100% sure what I am talking about here lol. So I may be wrong. –  Jeff Jul 13 '12 at 6:25
1  
Thanks for confirming that. FYI, different SRIDs cause STOverlaps to fail, among other things. But, aside from that, STOverlap will enable you to find what you're looking for. My only comment is that MakeValid actually changes the actual data by removing redundant points and trimming bowties etc. That is why it's not prefixed with 'ST' which signifies conformance to OGC standards. For a large enough dataset this could have a major influence on your results. –  Geoist Jul 13 '12 at 6:26
1  
SELECT DISTINCT geom.STSrid FROM gov and SELECT DISTINCT geom.STSrid FROM poa. And I'm not saying ignore the where clause, I'm just suggesting that your data has invalid geometries that are messing up the queries. Are they publicly available datasets that I could have a look at by chance? –  Geoist Jul 13 '12 at 6:30

2 Answers 2

up vote 3 down vote accepted
+50

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.

share|improve this answer
    
so I better use ST_Intersects to determine what government boundaries and post boundaries "intersect". –  Jeff Jul 17 '12 at 2:20
    
I must admit, I'm not totally clear what you're trying to do. I gather that you have two polygon shapefiles, gov.shp (large) regions, and post.shp of (smaller) postcodes. What result are you looking for? –  Micha Jul 17 '12 at 7:49
    
I have two boundaries (shape files), one is for local government boundary and one is post code boundary, both could be multiple polygon. And government boundary is bigger than post code boundary and post code may across more than one government boundary. I am trying to find out which post codes are within which government boundaries. I hope it's clear about what I am trying to do. Thanks. –  Jeff Jul 18 '12 at 1:10

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?

share|improve this answer
    
ST_Intersects is exactly what I am doing I've tried with both Intersects and Overlaps both in SQL Server 2008 R2 and PostGIS as I mentioned in my question. –  Jeff Jul 19 '12 at 0:30
    
OK, and are you not getting the results you want? –  Micha Jul 19 '12 at 19:28
    
Only because I am trying in both SQL Server/PostgreSQL PostGIS/Google Fusion Tables and I am getting different results back so. But I think I am getting there. Thanks heaps. –  Jeff Jul 19 '12 at 23:24

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.