Strange difference with ST_DWithin - Geography vs Geometry

user_15371779 0
Viewed 323

(EDIT: Sorry. Forgot to mention. Postgresql 13, PostGIS: 3.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1. Sorry about that).

The following 2 queries, in my opinion, should return "true" in both cases. The point is clearly within the polygon (square). Displaying both certainly indicates so. I am starting to wonder what fondamental thing I am missing here.

To be honest, I have been circling around this issue for years, as I am never able to find something "final" to understand what is happening. In this case, the point is somewhat near the edge of the polygon, but from time to time I had the same issues with points that were well within the polygon.

The only difference between the 2 queries is that one uses geometry the other uses geography. As you can see, we use 0.0 for the distance value. This is like an intersects.

OH, and by the way, using intersects yields the same results. I am very confused about this.

select 
ST_Dwithin(
St_GeometryFromText(
  'POLYGON((-95.979486 41.676556,-81.432269 41.676556,-81.432269 30.196043,-95.979486 30.196043,-95.979486 41.676556))', 4326),
ST_GeometryFromText('POINT (-87.70551 30.24481)', 4326),
  0
)

select 
ST_Dwithin(
St_GeographyFromText(
  'POLYGON((-95.979486 41.676556,-81.432269 41.676556,-81.432269 30.196043,-95.979486 30.196043,-95.979486 41.676556))'),
St_GeographyFromText('POINT (-87.70551 30.24481)'),
  0
)

EDIT. Sorry, not sure of the right etiquette. Following the great responses, I am adding some more context, in case somebody finds it interesting.

I am working on a "world" dataset containing country polygons. Since some of these are huge, I have elected to st_subdivide them.

![st_subdivide usa] https://imgur.com/IQRz8vg

And when zooming we can see EXACTLY what is explained in the answer below:

![zoom st_subdivide usa] https://imgur.com/fgmpBx0

So wow, perfect explanation, thank you very much.

I know that I should not use geography, well, "should not". However, I also need to work a lot with distances and expansion, and geom(4326) is quite annoying with units I find.

postgis postgresql
2 Answers

First of all, you need to understand the difference between geometry and geography.

  1. Geometry, where it assumes all of your data lives on a Cartesian plane.
  2. Geography, where it assumes that your data is made up of points on the earth's surface.

The problem is the following. When you are projecting the geometry, PostGIS is assuming that the edges of the geometry are not planar (since the earth's surface is not planar) and transforming them to spherical (this process is called tessellation). Therefore, when you try to intersect them the point is not inside de polygon.

Planar vs Spherical Edges

If the geographies are small enough, the difference between spherical and planar edges can be ignored, but in your case, the polygon is a huge one.

You are projecting your geometry on the EPSG 4326 in order to project it on the earth's surface, however although the point is inside the square in a planar coordinate system, it is not in a spherical one.

To sum up, you need to know what type, geometry or geography, works depending on your use case.

You can solve this issue by forcing the edges to be planar (in some implementations of PostGIS as BigQuery):

St_GeogFromText(
    'POLYGON((-95.979486 41.676556,-81.432269 41.676556,-81.432269 30.196043,-95.979486 30.196043,-95.979486 41.676556))',
     planar => True)
user_125365400

In addition to the answer from @sergiomahi - which should be accepted - I would like to add a few examples to make his point clearer. The issue isn't the type geography itself, but in which surface the features are projected.

Given your polygon and point, if we ST_Distance using geometry the distance is zero, since it will use a Cartesian plane, but using geography it returns 16524.35577956 metres due to the spherical edges and the size of your polygon, as pointed out by @sergiomahi. So, they no longer overlap.

WITH j (pol,poi) AS (
  VALUES
   ('SRID=4326;POLYGON((-95.979486 41.676556,-81.432269 41.676556,-81.432269 30.196043,-95.979486 30.196043,-95.979486 41.676556))',
    'SRID=4326;POINT (-87.70551 30.24481)')
)
SELECT 
  ST_Distance(poi::geometry,pol::geometry) AS geom_dist, 
  ST_Distance(poi::geography,pol::geography) AS geog_dist 
FROM j;



geom_dist |   geog_dist    
-----------+----------------
         0 | 16524.35577956

However, if you stick to geometry and tell PostGIS that the distance is supposed to be calculated using a sphere you get similar results as with geography:

WITH j (pol,poi) AS (
 VALUES
   ('SRID=4326;POLYGON((-95.979486 41.676556,-81.432269 41.676556,-81.432269 30.196043,-95.979486 30.196043,-95.979486 41.676556))',
    'SRID=4326;POINT(-87.70551 30.24481)'))
SELECT 
  ST_Distance(poi::geography,pol::geography) AS geog_dist, 
  ST_DistanceSphere(poi::geometry,pol::geometry) AS geom_dist_sphere
FROM j;

   geog_dist    | geom_dist_sphere 
----------------+------------------
 16524.35577956 |   16574.61755546
(1 Zeile)
user_22753880
Related