Search code examples
geospatialpostgis

most precise way to calculate areas (of buildings) in postgis

i have a postgis table buildings containing house perimeters located in Germany, with srid 3857, and i would like to calculate the area of each house.

my assumption is that

select st_area(st_transform(geom, utmzone(st_centroid(geom)))) as area from buildings

returns the most precise results, but i would like to have that confirmed by you.

running

select utmzone(st_centroid(geom)), count(*) from buildings group by utmzone(st_centroid(geom));

gives me:

enter image description here

in the following query, i employ a few more ways to calculate the areas:

select
    rel_err(st_area(geom), area) a,
    rel_err(st_area(st_transform(geom, 4326)), area) b,
    rel_err(st_area(st_transform(geom, 4326)::geography), area) c,
    rel_err(st_area(st_transform(geom, 25832)), area) d,
    rel_err(st_area(st_transform(geom, 25833)), area) e,
    rel_err(st_area(st_transform(geom, 31467)), area) f,
    rel_err(st_area(st_transform(geom, 32631)), area) g,
    rel_err(st_area(st_transform(geom, 32632)), area) h
from
(select st_area(st_transform(geom, utmzone(st_centroid(geom)))) as area, geom from buildings) q;

The resulting output looks like this:

enter image description here

(where rel_err(x, a) := (x-a)/a)

so here is my question: is the calculated value area really the best approximation to the area of the house perimeters, or is the value underlying one of the columns c to h even more precise? or is there another even more precise calculation method?

Solution

  • The answer is c, geography-based distance over a spheroid.

    SRID 3857 is pretty bad at distances and areas the further north you go. It thinks the poles are lines as long as the equator. See what happens when you ask about a location 200 meters east from the North Pole:

    with the_north_pole as (select st_point(0,90,4326) geom)
select ST_Distance( st_transform(geom,3857)
                   ,ST_Translate(st_transform(geom,3857),200,0))--moving 200m East
from the_north_pole ;
    st_distance
    200

    In reality, 200 meters east from the North Pole is still the North Pole - you're just spinning in place, around the point. If you use proper geography, it knows that if you go 180 degrees east, you haven't moved:

    select ST_Distance( st_point(0  ,90,4326)::geography
                   ,st_point(180,90,4326)::geography
                   ,use_spheroid=>true );
    st_distance
    0

    Your area of interest isn't as extreme, but this effect still applies: all your flat projections are stretching things further north, which is why you see their surface areas shrink when you deproject them to geography(polygon,4326). The larger zone you use, the more distorted everything typically gets, so

    • 3857 is your most distorted result,
    • that's followed by the "wide" UTM zones 2583x,
    • narrow strips of 3263x
    • 31467 strikes a good balance as a local projection

    Additionally, in reality all large "flat" surfaces are actually curved to match the curvature of the globe and geography-based calculations account for this, while projected geometries think it's all perfectly flat.

    The price is of course performance: geography calculations over a sphere are expensive, over a spheroid even more so. This expense might not really be worth it if you consider that the locations of vertices of your buildings are likely further from the ground truth than the differences between the ST_Area() results in these cases are from each other.