[postgis-tickets] [PostGIS] #2918: Use GeographicLib functions for ST_Azimuth, ST_Distance and ST_Project
PostGIS
trac at osgeo.org
Tue Sep 16 04:25:06 PDT 2014
#2918: Use GeographicLib functions for ST_Azimuth, ST_Distance and ST_Project
---------------------+------------------------------------------------------
Reporter: mwtoews | Owner: pramsey
Type: patch | Status: new
Priority: medium | Milestone: PostGIS 2.2.0
Component: postgis | Version: trunk
Keywords: |
---------------------+------------------------------------------------------
Comment(by mwtoews):
More performance and precision results ....
= Test data =
First, the `testgeod` table has been updated from yesterday.
{{{
drop table testgeod;
create unlogged table testgeod (
id serial primary key,
lat1 numeric, -- degrees
lon1 numeric, -- always 0
azi1 numeric, -- azimuth, in degrees
lat2 numeric,
lon2 numeric,
azi2 numeric,
s12 numeric, -- distance, in metres
a12 numeric,
m12 numeric,
"S12" numeric, -- area under the geodesic, in square metres
geog1 geography(Point,4326),
geog2 geography(Point,4326),
polyg geography(Polygon,4326), -- Karney (2013) Fig. 1, quadrilateral
AFHB
area double precision,
azimuth double precision,
distance double precision,
e real, -- first eccentricity: circle is zero, one is thin
antipode_fraction real
);
copy testgeod (lat1, lon1, azi1, lat2, lon2, azi2, s12, a12, m12, "S12")
from '/path/to/GeodTest.dat'
delimiter ' ';
update testgeod set
geog1 = ST_MakePoint(lon1, lat1)::geography,
geog2 = ST_MakePoint(lon2, lat2)::geography,
polyg = ST_MakePolygon(ST_MakeLine(
ARRAY[ST_MakePoint(lon1, lat1), -- A
ST_MakePoint(lon1, 0), -- F
ST_MakePoint(lon2, 0), -- H
ST_MakePoint(lon2, lat2), -- B
ST_MakePoint(lon1, lat1)])), -- close with A
area = abs("S12"),
azimuth = radians(azi1),
distance = s12,
e = case when abs(lat1 - lat2) < abs(lon1 - lon2)
then sqrt(1 - pow(lat1 - lat2, 2)/pow(lon1 - lon2, 2))
else sqrt(1 - pow(lon1 - lon2, 2)/pow(lat1 - lat2, 2)) end,
antipode_fraction = s12 / 20003931.459;
}}}
The test column `e` ("first eccentricity" of the bounds) is used to tell
very thin quadrilaterals (1.0) to roundish (0.0). The `antipode_fraction`
is near 0.0 for short geodesics, and 1.0 for ones that are 180°. These are
necessary to filter out some end-member tests for fair comparisons with
Vincenty's algorithms. The test data has a wide mixture of geodesics,
including a number of corner cases like near-antipodes. Please keep this
in mind for both interpretations for both timings and analysis of
precision distributions, as they are probably biased against Vincenty's.
= Direct geodesic calculation (`ST_Project`) - update =
The only change is to avoid casting `s12` from numeric, and just use
`distance` (double precision).
{{{
drop table if exists f;
select id, st_project(geog1, distance, azimuth) into unlogged f from
testgeod;
}}}
This boots the performance of each test to 1871 ms vs 2742 ms, or 1.46
times slower for Karney.
The precision results are the same as before, but to describe the %
accuracy at a few other arbitrary cutoffs:
{{{
select
avg(case when diff_distance <= 1e-2 then 100 else 0 end)::numeric(8,2)
|| '%' as cm,
avg(case when diff_distance <= 1e-3 then 100 else 0 end)::numeric(8,2)
|| '%' as mm,
avg(case when diff_distance <= 1e-6 then 100 else 0 end)::numeric(8,2)
|| '%' as µm,
avg(case when diff_distance <= 1e-8 then 100 else 0 end)::numeric(8,2)
|| '%' as ten_nm
from (
select st_distance(st_project, geog2) as diff_distance
from f join testgeod g on g.id=f.id
) as f;
}}}
||Method|| 1 cm || 1 mm || 1 µm || 10 nm ||
||Vincenty||100.00%||96.33%||28.05%||24.90%||
||Karney||100.00%||100.00%||100.00%||99.98%||
= Inverse geodesic calculations =
== Azimuth ==
Test performance of `ST_Azimuth`:
{{{
drop table if exists f;
select id, st_azimuth(geog1, geog2) into unlogged f from testgeod;
}}}
With these tests, Karney's is 3 times faster than Vincenty's, which is the
opposite of what I've found earlier, and I think it is due to the dataset
which has challenging geodesics that require more iteration with
Vincenty's. Compare 11384 ms to 3786 ms for Karney's.
Test accuracy of `ST_Azimuth`:
{{{
select f.id, s12, antipode_fraction, e, degrees(st_azimuth) as
st_azimuth_degrees, azi1,
(((degrees(st_azimuth) - azi1 + 180.0)::numeric % 360.0) - 180.0) as
diff_degrees
from f join testgeod g on g.id=f.id
-- where antipode_fraction between 0.004 and 0.996
order by abs(((degrees(st_azimuth) - azi1 + 90.0)::numeric % 180.0) -
90.0) desc
limit 10;
}}}
With the near-antipodes included (no filter), compare maximum absolute
errors 89.992633730349° to 0.031697236597° error. Removing the comment to
apply the filter to remove near-antipodes, compare 0.000006260972° to
0.000000000003°. And the % within precision table:
{{{
select
avg(case when diff_degrees <= 90 then 100 else 0 end)::numeric(8,3) ||
'%',
avg(case when diff_degrees <= 1e+0 then 100 else 0 end)::numeric(8,3) ||
'%',
avg(case when diff_degrees <= 1e-3 then 100 else 0 end)::numeric(8,3) ||
'%',
avg(case when diff_degrees <= 1e-6 then 100 else 0 end)::numeric(8,3) ||
'%',
avg(case when diff_degrees <= 1e-9 then 100 else 0 end)::numeric(8,3) ||
'%',
avg(case when diff_degrees <= 1e-12 then 100 else 0 end)::numeric(8,3)
|| '%',
avg(case when diff_degrees <= 1e-15 then 100 else 0 end)::numeric(8,3)
|| '%'
from (
select
abs(((degrees(st_azimuth) - azi1 + 180.0)::numeric % 360.0) - 180.0)
as diff_degrees
from f join testgeod g on g.id=f.id
) as f;
}}}
||Method||90°||1°||1e-3°||1e-6°||1e-9°||1e-12°||1e-15°||
||Vincenty||91.675%||89.842%||84.959%||78.527%||73.551%||32.790%||31.574%||
||Karney||100.000%||100.000%||99.995%||86.328%||82.551%||61.332%||57.107%||
== Distance ==
Test performance of `ST_Distance`:
{{{
drop table if exists f;
select id, st_distance(geog1, geog2) into unlogged f from testgeod;
}}}
The results are essentially the same for azimuth, since they use the same
underlying methods.
For the accuracy:
{{{
select f.id, s12, antipode_fraction, e, distance, st_distance,
(st_distance - distance) as diff_distance
from f join testgeod g on g.id=f.id
-- where antipode_fraction between 0.004 and 0.996
order by abs(st_distance - distance) desc
limit 10;
}}}
With the near-antipodes included (no filter), compare maximum absolute
errors 133912.525 m to 1.49e-8 m error. Removing the comment to apply the
filter to remove near-antipodes, the error of Vincenty's reduces to
0.01953874 m, and Karney's remains the same, at 14.9 nm. And the % within
precision table:
{{{
select
avg(case when diff_distance <= 1e+5 then 100 else 0 end)::numeric(8,3)
|| '%',
avg(case when diff_distance <= 1e+3 then 100 else 0 end)::numeric(8,3)
|| '%',
avg(case when diff_distance <= 1e+0 then 100 else 0 end)::numeric(8,3)
|| '%',
avg(case when diff_distance <= 1e-3 then 100 else 0 end)::numeric(8,3)
|| '%',
avg(case when diff_distance <= 1e-6 then 100 else 0 end)::numeric(8,3)
|| '%',
avg(case when diff_distance <= 1e-8 then 100 else 0 end)::numeric(8,3)
|| '%',
avg(case when diff_distance <= 1e-9 then 100 else 0 end)::numeric(8,3)
|| '%'
from (
select
abs(st_distance - distance) as diff_distance
from f join testgeod g on g.id=f.id
) as f;
}}}
||Method|| 500 km || 1 km || 1 m || 1 mm || 1 µm || 10 nm || 1 nm ||
||Vincenty||99.659%||97.404%||96.400%||65.387%||19.603%||6.425%||3.135%||
||Karney||100.000%||100.000%||100.000%||100.000%||100.000%||99.873%||58.909%||
== Area ==
The geodesic test data use a the quadrilateral AFHB (Fig. 1, Karney 2013),
which is area under the geodesic to the equator. Unfortunately, the
current algorithm in PostGIS cannot use points in the same hemisphere
("ERROR: ptarray_area_spheroid: cannot handle ptarray that crosses
equator"), so only 325088 tests or 65.02% can be used. Here are the
performance tests:
{{{
drop table if exists f;
select id, st_area(polyg) into unlogged f from testgeod
where (lat1 > 0 and lat2 < 0) or (lat1 < 0 and lat2 > 0);
}}}
The performance is identical, about 1950 ms. I think I need to dig around
a bit more to find out why, as I'm seeing really odd behaviour. This
Polygon has a NaN area for both unpatched and patched:
{{{
select ST_Area('POLYGON((0 78.703946026663,0 0,179.999997913235
0,179.999997913235 -33.0888306884702,0 78.703946026663))'::geography);
}}}
Looking at the debug messages, it seems to be happy to calculate the dot
product instead of the area
{{{
NOTICE: [lwgeom.c:lwgeom_set_srid:1499] entered with srid=4326
NOTICE: [lwgeom.c:lwgeom_is_empty:1277] lwgeom_is_empty: got type Polygon
NOTICE: [lwgeom.c:lwgeom_is_empty:1277] lwgeom_is_empty: got type Polygon
NOTICE: [lwgeodetic.c:edge_point_side:655] dot product 3.64209202e-08
NOTICE: [lwgeodetic.c:edge_point_side:655] dot product -3.97170884e-09
NOTICE: [lwgeom_transform.c:PROJ4SRSCacheDelete:139] deleting projection
object (0x7f35d11ba360) with MemoryContext key (0x7f35d11d8c10)
}}}
whereas with `select ST_Area('POLYGON((0 0, 1 1, 1 0, 0 0))'::geography)`
the messages in this location are:
{{{
NOTICE: [lwgeom.c:lwgeom_set_srid:1499] entered with srid=4326
NOTICE: [lwgeom.c:lwgeom_is_empty:1277] lwgeom_is_empty: got type Polygon
NOTICE: [lwgeom.c:lwgeom_is_empty:1277] lwgeom_is_empty: got type Polygon
NOTICE: [lwspheroid.c:ptarray_area_spheroid:143] geod_polygon_addpoint 1:
1 1
NOTICE: [lwspheroid.c:ptarray_area_spheroid:143] geod_polygon_addpoint 2:
0 1
NOTICE: [lwspheroid.c:ptarray_area_spheroid:143] geod_polygon_addpoint 3:
0 0
NOTICE: [lwspheroid.c:ptarray_area_spheroid:151] geod_polygon_compute
area: -6154854786.72
NOTICE: [lwgeom_transform.c:PROJ4SRSCacheDelete:139] deleting projection
object (0x7f35d11ba360) with MemoryContext key (0x7f35d11d8c10)
}}}
so why the dot product instead of ptarray_area_spheroid?
--
Ticket URL: <http://trac.osgeo.org/postgis/ticket/2918#comment:15>
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.
More information about the postgis-tickets
mailing list