4

I'm running ST_Area on two different NaturalEarth datasets but am getting different results when the geographies are clipped before performing the area calculation.

First, I've loaded the 10m "admin0" dataset (http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_countries.zip) into a table named countries.

Second, I've loaded the 10m "land" dataset (http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_land.zip) into a table named land.

The area calculation is basically the same for both:

SELECT sum(ST_Area(geom::geography)) AS area FROM land;
     area     
-----------------
 146826440837593


SELECT sum(ST_Area(geom::geography)) AS area FROM countries;
       area       
-----------------
 146826603523460

Examining the data, it appears that both datasets are the same except that the land table has one large multipolygon, whereas countries has many individual multipolygons (ie one for each country in the world).

The issue arises when I clip both to a latitude of +/-52 degrees, in which case the areas I calculate for each diverge quite a bit. This is the query:

SELECT
sum(
  ST_Area(
    ST_Intersection(
      ST_GeomFromEWKT(
        'SRID=4326;POLYGON((-180 -52, -180 52, 180 52, 180 -52, -180 -52))'),
      geom
    )::geography
  )
) AS land_area
FROM
 land;

And it produces these results:

    land_area       
-----------------
 119910798321370

  countries_area       
-----------------
 106619873643318

Any ideas as to why this would be the case? I've also visualised the clipped geometries for both tables and they appear the same to me and I don't see the same area differences when I clip to latitudes like +/-80 degrees.

All the geometries appear to be valid before the intersection as well as after the intersection. For example, the following query produces no results:

SELECT
 *
FROM
 countries
WHERE
NOT ST_IsValid(
      ST_Intersection(
        ST_GeomFromEWKT(
          'SRID=4326;POLYGON((-180 -52, -180 52, 180 52, 180 -52, -180 -52))'),
        geom
      )
    );

UPDATE: If I do the calculation without converting to geography, the numbers do not diverge like they do above:

    land_area     
------------------
 9942.37708696352

   country_area   
------------------
 9942.38832452307
8
  • 1
    Natural Earth data can be quite dirty, have you run isvalid and makevalid on them? Commented Aug 15, 2016 at 13:04
  • @PaulRamsey - I didn't find anything invalid when checking the original geometries or the geometries post-intersection (updated question). Thanks!
    – bosth
    Commented Aug 15, 2016 at 14:32
  • 1
    You might also look at the presence (or not) of vertices, particularly along long, straight lines.
    – mkennedy
    Commented Aug 15, 2016 at 17:28
  • @mkennedy Could you explain a little more? Once the clip is applied along the +/-52 latitude, there will be long straight lines without vertices. But how would this affect the results of the area calculation?
    – bosth
    Commented Aug 15, 2016 at 17:40
  • 1
    Yeah, if there are lines w/ extra vertices, they'll line up differently as great circles than as cartesian lines. You'll looking at the shapes in cartesian and saying "they are the same", but maybe they aren't the same in great circle world... Commented Aug 15, 2016 at 18:34

1 Answer 1

5

The problem was around calculating the area of a large polygon with straight lines that had been converted to a geography type.

Using ST_Segmentize will ensure that there are more intermediary segments in the lines that will provide a better representation of the cartesian lines once you are working with geography:

SELECT
sum(
  ST_Area(
    ST_Segmentize(
      ST_Intersection(
        ST_GeomFromEWKT(
          'SRID=4326;POLYGON((-180 -52, -180 52, 180 52, 180 -52, -180 -52))'),
        geom
      ), 0.1
    )::geography
  )
) AS land_area
FROM
 land;

The reason that the countries table wasn't affected is that it has political boundaries doing some rudimentary "segmentation" of the long lines.

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.