ST_Area NaN on some polygons

[rfc2616]

In one of our African datasets, we need to calculate ST_Area on the sphere for all polygons. ST_Area is returning NaN for some, but not all, polygons. The affected polygons appear to be equator crossing but otherwise there is no obvious connection and not all equator crossing polygons are affected.

We're confirming the issue on running PostGIS instances as far back as 2.1.2, but chronologically this only surfaced recently for us -- at some point in the past year or two these calculations used to pass without any NaN.

The failing polygons vary depending on the exact system and PostGIS version and seem somewhat bizarre and random.

The results below come from 2 systems, my Mac laptop and a Linux server and are slightly different but both broken.

Laptop:

PostgreSQL 9.4.5 on x86_64-apple-darwin15.0.0, compiled by Apple LLVM version 7.0.0 (clang-700.1.76), 64-bit, POSTGIS="2.2.0 r14208" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.9.2, 08 September 2015" GDAL="GDAL 1.11.3, released 2015/09/16" LIBXML="2.9.3" LIBJSON="0.12" RASTER

Server:

PostgreSQL 9.3.10 on x86_64-unknown-linux-gnu, compiled by gcc (Ubuntu 4.8.2-19ubuntu1) 4.8.2, 64-bit, POSTGIS="2.1.2 r12389" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.1, released 2013/08/26" LIBXML="2.9.1" LIBJSON="UNKNOWN" RASTER

This test polygon is taken from #2918 and still yields NaN on both systems:

select ST_Area('POLYGON((0 78.703946026663,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026663))'::geography,false);

Here I made 1 digit difference in the starting point (...662 vs ...663)

select ST_Area('POLYGON((0 78.703946026662,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026662))'::geography,false);

This gets a result, 127516467322130, on the server ... but NaN on the laptop.

A different one digit change, ...664:

select ST_Area('POLYGON((0 78.703946026664,0 0,179.999997913235 0,179.999997913235 -33.0888306884702,0 78.703946026664))'::geography,false);

yields 127516467322130 on the *laptop* ... but NaN on the *server*.

I have no idea how to even go about troubleshooting this, but it's causing us a major headache. Happy to step up with whatever support we can offer to help isolate it! Thanks in advance for any help.

[pramsey]

Any reason you're calculating on the sphere in particular? (Not that it matters, the old version would be sphere only anyways, probably.) The results even when you don't get NaN are pretty bad, even back when using postgis 2.0, for example:

select ST_Area('POLYGON((0 78.703946026662,0 0,179.999997913235
 0,179.999997913235 -33.0888306884702,0
 78.703946026662))'::geography,false);
     st_area     
-----------------
 127430010458336

I'd expect changes from 2.1 to 2.2, as we changed to new spherical area code. I'm surprised to see changes from 2.0 to 2.1, but I see that also. Anyways, I'm not sure the answers are right, even when you get non-NaN answers. Will have to confirm that.

[rfc2616]

We were using the sphere for continuity with old published results. I just changed all the queries to use the spheroid and they seem to be running fine on 2.2.0 and within a small tolerance of our old sphere-based results, so that will get us out of short-term trouble. You're also correct that all the sphere based calculations are super crazy right now, and not just on the equator.

select analysis_name, analysis_year, region, range_quality, category, country, ST_Area(st_intersection::geography,true) spheroid, ST_Area(st_intersection::geography,false) sphere from survey_range_intersections where analysis_year=2007 order by region, country, category;
   analysis_name   | analysis_year |     region      | range_quality | category |           country            |      spheroid       |      sphere       
-------------------+---------------+-----------------+---------------+----------+------------------------------+---------------------+-------------------
 2013_africa_final |          2007 | Central Africa  | Known         | C        | Cameroon                     |    44861.9693270706 |  3595587.39688914
 2013_africa_final |          2007 | Central Africa  | Known         | C        | Cameroon                     |    824073902.362827 |               NaN
 2013_africa_final |          2007 | Central Africa  | Known         | C        | Cameroon                     |    678834776.222033 |               NaN
 2013_africa_final |          2007 | Central Africa  | Known         | D        | Cameroon                     |    503424735.671734 |  504085360.780887
 2013_africa_final |          2007 | Central Africa  | Known         | D        | Cameroon                     |    430138925.122056 |               NaN
 2013_africa_final |          2007 | Central Africa  | Known         | D        | Cameroon                     |    653881544.632783 |               NaN
 2013_africa_final |          2007 | Central Africa  | Possible      | D        | Cameroon                     |    734248715.634702 |               NaN
 2013_africa_final |          2007 | Central Africa  | Known         | D        | Cameroon                     |    459.157410163432 |  1815169.69392417
 2013_africa_final |          2007 | Central Africa  | Known         | D        | Cameroon                     |     80418.939783203 |               NaN
 2013_africa_final |          2007 | Central Africa  | Possible      | E        | Cameroon                     |    1977273.09717699 |               NaN
 2013_africa_final |          2007 | Central Africa  | Possible      | E        | Cameroon                     |    1530669683.44256 |               NaN
 2013_africa_final |          2007 | Central Africa  | Known         | E        | Cameroon                     |    3.43081077933311 |  5.26343760021617

[dbaston]

I get 127516467322130 for all three polygons if I make a small tweak in lwgeodetic.c, line 919:

double sphere_distance_cartesian(const POINT3D *s, const POINT3D *e)
{
	return acos(FP_MIN(1.0, dot_product(s, e)));
}

[pramsey]

Make it so.

[dbaston]

I applied the above change back to 2.1 (r14485, r14486, r14487). (If I'm not mistaken, the spheroid area routines changed in 2.2, but the sphere routines did not.)

This resolves the NaN issues, or at least any that I can find. I'm not sure about "all the sphere-based calculations are super-crazy" though. Maybe you can post some test polygons to a separate ticket, rfc2616?

[dbaston]

Also fixed in 2.0 at r14514