[postgis-users] PostGIS transforms st_dwithin to && on geography column without preserving the custom spatial reference system

[Michael Loveridge]

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 && bounding box 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). But it seems that when using 
st_dwithin, the search is being done in srid 4326 instead of my custom one.

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?