Ticket #450: 450-1.patch

File 450-1.patch, 3.5 KB (added by pramsey, 2 years ago)

Fix for spheroid area

  • liblwgeom/lwspheroid.c

     
    430430 
    431431                LWDEBUGF(4, "in_south %d", in_south); 
    432432 
     433                LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) ); 
     434 
    433435                if ( crosses_dateline(&a, &b) ) 
    434436                { 
    435437                        double shift; 
     
    439441                        else 
    440442                                shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */ 
    441443 
     444                        LWDEBUGF(4, "shift: %.8g", shift); 
     445                        LWDEBUGF(4, "before shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon); 
    442446                        point_shift(&a1, shift); 
    443447                        point_shift(&b1, shift); 
     448                        LWDEBUGF(4, "after shift a1(%.8g %.8g) b1(%.8g %.8g)", a1.lat, a1.lon, b1.lat, b1.lon); 
     449                         
    444450                } 
    445451 
    446                 LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(&a, &b) ); 
    447452 
    448453                delta_lon = fabs(b1.lon - a1.lon); 
    449454 
    450                 LWDEBUGF(4, "(%.18g %.18g) (%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon); 
     455                LWDEBUGF(4, "a1(%.18g %.18g) b1(%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon); 
    451456                LWDEBUGF(4, "delta_lon %.18g", delta_lon); 
    452457                LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance); 
    453458 
  • liblwgeom/cunit/cu_geodetic.c

     
    920920        CU_ASSERT_DOUBLE_EQUAL(a1, 89.7211470368, 0.0001); /* sphere */ 
    921921        CU_ASSERT_DOUBLE_EQUAL(a2, 89.8684316032, 0.0001); /* spheroid */ 
    922922 
    923  
    924923        /* Big-ass polygon */ 
    925924        lwg = lwgeom_from_ewkt("POLYGON((-2 3, -2 4, -1 4, -1 3, -2 3))", PARSER_CHECK_NONE); 
    926925        lwgeom_calculate_gbox_geodetic(lwg, &gbox); 
     
    930929        CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */ 
    931930        CU_ASSERT_DOUBLE_EQUAL(a2, 12286574431.9, 10.0); /* spheroid */ 
    932931 
     932        /* One-degree square */ 
     933        lwg = lwgeom_from_ewkt("POLYGON((8.5 2,8.5 1,9.5 1,9.5 2,8.5 2))", PARSER_CHECK_NONE); 
     934        lwgeom_calculate_gbox_geodetic(lwg, &gbox); 
     935        //a1 = lwgeom_area_sphere(lwg, &gbox, &s); 
     936        a2 = lwgeom_area_spheroid(lwg, &gbox, &s); 
     937        //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); 
     938        CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */ 
     939        CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */ 
     940 
     941        /* One-degree square *near* dateline */ 
     942        lwg = lwgeom_from_ewkt("POLYGON((178.5 2,178.5 1,179.5 1,179.5 2,178.5 2))", PARSER_CHECK_NONE); 
     943        lwgeom_calculate_gbox_geodetic(lwg, &gbox); 
     944        //a1 = lwgeom_area_sphere(lwg, &gbox, &s); 
     945        a2 = lwgeom_area_spheroid(lwg, &gbox, &s); 
     946        //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); 
     947        CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */ 
     948        CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */ 
     949 
     950        /* One-degree square *across* dateline */ 
     951        lwg = lwgeom_from_ewkt(" POLYGON((179.5 2,179.5 1,-179.5 1,-179.5 2,179.5 2))", PARSER_CHECK_NONE); 
     952        lwgeom_calculate_gbox_geodetic(lwg, &gbox); 
     953        a1 = lwgeom_area_sphere(lwg, &gbox, &s); 
     954        a2 = lwgeom_area_spheroid(lwg, &gbox, &s); 
     955        //printf("\nsphere: %.12g\nspheroid: %.12g\n", a1, a2); 
     956        //CU_ASSERT_DOUBLE_EQUAL(a1, 12341436880.1, 10.0); /* sphere */ 
     957        CU_ASSERT_DOUBLE_EQUAL(a2, 12304814950.073, 100.0); /* spheroid */ 
     958 
    933959} 
    934960 
    935961 
  • liblwgeom/lwgeodetic.c

     
    131131        double lon = p->lon + shift; 
    132132        if ( lon > M_PI ) 
    133133                p->lon = -1.0 * M_PI + (lon - M_PI); 
     134        else 
     135                p->lon = lon; 
    134136        return; 
    135137} 
    136138