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.

1 Answer 1

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.

3
  • 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
    Commented Jan 14, 2018 at 3: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
    Commented Jan 14, 2018 at 3: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
    Commented Jan 14, 2018 at 4:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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