CODEWITHSUNDEEP
postgresqlpostgis
loveridge
loveridge

Reputation: 23

PostGIS transforms st_dwithin to && on geography column without preserving the custom spatial reference system

When I use st_dwithin on a geography column with a GIST index, the explain plan shows that the narrowing condition on the index scan is changed into a && condition. With this different condition, the result set is missing rows that are included if a st_dwithin is performed using a full table scan when there is no index.

The condition in the explain plan without an index is

st_dwithin(
  geog,
  '010100002010A4000000000000006066400000000000000000'::geography,
  '3'::double precision,
  false
)

but the explain plan indicates that it is changed into the following when the column has a GIST index

geog && _st_expand(
  '010100002010A4000000000000006066400000000000000000'::geography,
  '3'::double precision
)

The geography column has a custom spatial reference system that is not being preserved in the transformation into a && condition. The custom spatial reference system is a sphere where the radius is such that one degree is equal to one meter of great-circle distance. So that a point of (0 0) is one distance away from (0 1).

The setup code is below

-- Insert the custom spatial reference system into the db.
insert into spatial_ref_sys values (42000, 'customsrs', 1,
     'GEOGCS[
       "Normal Sphere (r=57.2957795)",
       DATUM["unknown",
         SPHEROID["Sphere",57.29577951308232087679,0]
       ],
       PRIMEM["Greenwich",0],
       CS[ellipsoidal,2],
       AXIS["latitude",north],
       AXIS["longitude",east],
       UNIT["degree",0.0174532925199433]
     ]', '+proj=longlat +ellps=sphere +R=57.29577951308232087679 +no_defs');

-- Create a new table with a geography column in the new spatial reference system.
CREATE TABLE geographytest(gid serial PRIMARY KEY, geog geography(POINT, 42000));

-- Insert some data around the dateline
insert into geographytest (gid, geog) values
(1, 'srid=42000;POINT(179 0)'),
(2, 'srid=42000;POINT(178 0)'),
(3, 'srid=42000;POINT(-179 0)'),
(4, 'srid=42000;POINT(-179 90)'),
(5, 'srid=42000;POINT(0 0)');


-- Select all points within a distance of 3 from POINT(179 0).
--   The expected 3 points are returned, with st_distances of 0, 1 and 2.
select
 gid,
 st_distance(
   geog,
   'srid=42000;POINT(179 0)',
   false
 ),
 st_dwithin(
   geog,
   'srid=42000;POINT(179 0)',
   3,
   false
 )
from
 geographytest
where
 st_dwithin(
   geog,
   st_geogfromtext('srid=42000;POINT(179 0)'),
   3,
   false
 );

-- Create a GIST index on our geography column
CREATE INDEX geographytestindex ON geographytest USING gist (geog);
VACUUM analyze geographytest (geog);

-- Now select again using the same query.
--   Now only one result is returned, the row that has POINT(179 0).
--   The explain plan indicates that the index scan is using 
--   Index Cond: (geog && _st_expand('010100002010A4000000000000006066400000000000000000'::geography, '3'::double precision))
--   when performing the select.
select
 gid,
 st_distance(
   geog,
   'srid=42000;POINT(179 0)',
   false
 ),
 st_dwithin(
   geog,
   'srid=42000;POINT(179 0)',
   3,
   false
 )
from
 geographytest
where
 st_dwithin(
   geog,
   st_geogfromtext('srid=42000;POINT(179 0)'),
   3,
   false
 );

Increasing the st_dwithin distance to 230000 instead of 3 returns all three of the expected rows, since that is their distance in srid 4326.

How can I get PostGIS to use my custom spatial reference system in the query when an index is present on the column?

Upvotes: 1

Views: 316

Answers (1)

JGH
JGH

Reputation: 17906

The issue is with _st_expand. This function expands the bounding box of the given geography. Looking at the code, we can see that the desired distance is modified using WGS84 radius, so the custom CRS is indeed ignored. You can report the bug if you wish.

/* Read our distance value and normalize to unit-sphere. */
distance = PG_GETARG_FLOAT8(1);
/* Magic 1% expansion is to bridge difference between potential */
/* spheroidal input distance and fact that expanded box filter is */
/* calculated on sphere */
unit_distance = 1.01 * distance / WGS84_RADIUS;

Upvotes: 2

Related Questions