3

I have two tables in Postgres, both of type geometry(MultiPolygon,27700) and I would like to find any that intersect with an area of more than 100 square metres.

I know how to get any that intersect above an arbitrary threshold:

SELECT * 
FROM aTable a 
JOIN bTable b
WHERE ST_Area(ST_Intersection(a.geom, b.geom)) > 1;

However, the area returned here is in (I think) decimal degrees. How can I ensure it is always more than 100 square metres?

Perhaps if I do something like:

WHERE ST_Area(ST_Transform(ST_Intersection(i.wkb_geometry, w.wkb_geometry), 3857)) > 100

that would work, since the unit of EPSG:3857 is metres?

I've read this and this but I don't know what projection to pick.

Richard
  • 3,289
  • 6
  • 35
  • 66

1 Answers1

4

St_area returns cartesian area for geometries and geodesic areas for geographies.

If your coordinate system is in lat/lan, simply type cast your geometry to geography and your result will be in sq meters

ST_Area(ST_Intersection(a.geom, b.geom)::geography)

EDIT: For projected data, st_area() will always return areas in the linear unit of the projection itself. EPSG 27700 is the British National Grid projection. It is a Transverse Mercator projection based on the OSGB36 datum. Its linear unit is meter. More info here.

The function

ST_Area(ST_Intersection(a.geom, b.geom))

should give you directly an area in sq meters. If you want to double check, test this

ST_Area(ST_Transform(ST_Intersection(a.geom, b.geom),4326)::geography)

This transforms your coordinate from BNG to WGS84 and calculates the geodesic area. Results should be somewhat similar to the simple st_area() call above.

Thomas
  • 1,197
  • 5
  • 11
  • Thanks! Unfortunately I get ERROR 1: ERROR: Only lon/lat coordinate systems are supported in geography.. Table a has geometry of type geometry(Polygon,27700) and table b has geometry of type geometry(MultiPolygon,27700). – Richard Jan 14 '18 at 03:06
  • Maybe I need to transform the intersection before casting it to geography? ST_Area(ST_Transform(ST_Intersection(a.geom, b.geom), 3857)::geography) -- nope, that gives the same error. – Richard Jan 14 '18 at 03:10
  • Sorry, I wasn't familiar with the British National Grid and assumed from your question that it was in lat/lon. It is a Transverse Mercator projection and its linear unit is in meters. a simple st_area() call should do. I edited my answer as well. – Thomas Jan 14 '18 at 04:31