[postgis-devel] [PostGIS] #503: ST_Area on small polygons in EPSG:31467 gives incorrect results depending on 32-/64-system-architecture

[PostGIS]

#503: ST_Area on small polygons in EPSG:31467 gives incorrect results depending
on 32-/64-system-architecture
----------------------+-----------------------------------------------------
  Reporter:  mjansen  |       Owner:  pramsey      
      Type:  defect   |      Status:  new          
  Priority:  medium   |   Milestone:  PostGIS 1.5.2
 Component:  postgis  |     Version:  1.5.X        
Resolution:           |    Keywords:               
----------------------+-----------------------------------------------------
Comment (by nicklas):

 I think I have found something.
 I can see the same calculation gets different results from inside
 postgresql/postgis than if I run it in a standalone c program in the same
 environment.

 I tried with the two first coordinates in Paul's example above:
 In the area function the calculation is:
 ringarea += ( p1.x * p2.y ) - ( p1.y * p2.x );

 I put a lwnotice just before this calculation that looks like:
 lwnotice("ringarea before:%f, first part:%f, second part:%f, all:%f\n",
 ringarea, ( p1.x * p2.y ), ( p1.y * p2.x ), (( p1.x * p2.y ) - ( p1.y *
 p2.x )));

 in Ubuntu (9.10 32 bit) this gives:
 NOTICE:  ringarea before:0.000000, first part:20344027056218.062500,
 second part:20344027398598.898438, all:-342380.838655

 in windows (XP 32 bit on the same machine) this gives:
 NOTICE:  ringarea before:0.000000, first part:20344027056218.062000,
 second part:20344027398598.898000, all:-342380.835938


 Then I tried to run it in windows from a program like this:

 {{{
 #include <stdio.h>

 int
 main(){
 double x1 = 3498246.2563;
 double x2 = 5815493.2389;
 double y1 = 3498246.2472;
 double y2 = 5815493.1259;
   printf("first part=%f, second part = %f, all=%f\n", ( x1 * y2 ),( y1 *
 x2 ),(( x1 * y2 ) - ( y1 * x2 ))  );
 }
 }}}

 first part=20344027056218.062000, second part = 20344027398598.898000,
 all=-342380.838655

 So here I get exactly the same result in the end like ubuntu.

 Then I did put the same code in the end of the area function inside
 postgis and out putted through lwnotice.

 {{{
 .........
                 poly_area += ringarea;
         }
 /*My extra code*/
 double x1 = 3498246.2563;
 double x2 = 5815493.2389;
 double y1 = 3498246.2472;
 double y2 = 5815493.1259;
   lwnotice("first part=%f, second part = %f, all=%f\n", ( x1 * y2 ),( y1 *
 x2 ),(( x1 * y2 ) - ( y1 * x2 ))  );

 /*Back to original*/
         return poly_area;
 }
 }}}

 Then I had the wrong result again :
 NOTICE:  first part=20344027056218.062000, second part =
 20344027398598.898000, all=-342380.835938


 So something is different inside postgresql/postgis than outside.
 Can the truncation in the first and second part affect or is that only a
 limitation in the output.


 /Nicklas

-- 
Ticket URL: <http://trac.osgeo.org/postgis/ticket/503#comment:12>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.