postgresqlgispostgis
How to accurately determine ST_Intersects of point and geography near boundaries (ST_Intersects geography vs. geometry discrepancy)
I have a geography column of polygons. This column has been ST_Subdivided for point-in-polygon lookup performance.
When I perform an ST_Intersects between a point < 10m from the edge of a polygon with a geography type, I get false. But after casting to geometry types, I get true:
SELECT
ST_Intersects(ST_GeogFromText('POINT(-118.46131 34.31171)'), latlong_shape) AS intersects_geog,
ST_Intersects(ST_GeogFromText('POINT(-118.46131 34.31171)')::geometry, latlong_shape::geometry) AS intersects_geom
FROM region_ids_for_latlong_shapes
WHERE ST_Distance(ST_GeogFromText('POINT(-118.46131 34.31171)'), latlong_shape) < 10;
| intersects_geog | intersects_geom |
| --------------- | --------------- |
| false | true |
| false | false |
I'm aware that there are expected differences between ST_Intersects(geography, geography) and ST_Intersects(geometry, geometry), but I don't understand why they would be occurring in this instance:
- The geometries are valid
- Even if it were some great circle discrepancy, I'd expect the point to be intersecting one of the polygons, not outside both of them.
What is the most accurate way to determine which polygon this point intersects? I understand I can cast to geometry types, but I'm surprised that ST_Intersect's geography support isn't working here. Is there to do this accurately without casting?
Here is the point and the two closest shapes:
Some more diagnostics:
SELECT
ST_Intersects(ST_GeogFromText('POINT(-118.46131 34.31171)'), latlong_shape) AS geogs_intersect,
ST_Intersects(ST_GeogFromText('POINT(-118.46131 34.31171)')::geometry, latlong_shape::geometry) AS geoms_intersect,
ST_IsValid(latlong_shape::geometry) AS geom_isvalid,
ST_Relate(ST_GeogFromText('POINT(-118.46131 34.31171)')::geometry, latlong_shape::geometry) AS relate_geoms,
ST_Distance(ST_GeogFromText('POINT(-118.46131 34.31171)'), latlong_shape) AS geogs_distance,
ST_Distance(ST_GeogFromText('POINT(-118.46131 34.31171)'), latlong_shape, false) AS geogs_no_spheroid_distance,
ST_Distance(ST_GeogFromText('POINT(-118.46131 34.31171)')::geometry, latlong_shape::geometry) AS geoms_distance,
ST_3DDistance(ST_GeogFromText('POINT(-118.46131 34.31171)')::geometry, latlong_shape::geometry) AS geoms_3d_distance
FROM region_ids_for_latlong_shapes
WHERE ST_Distance(ST_GeogFromText('POINT(-118.46131 34.31171)'), latlong_shape) < 10
| geogs_intersect | geoms_intersect | geom_isvalid | relate_geoms | geogs_distance | geogs_no_spheroid_distance | geoms_distance | geoms_3d_distance |
| --------------- | --------------- | ------------ | ------------ | -------------- | -------------------------- | ----------------------- | ----------------------- |
| false | true | true | 0FFFFF212 | 3.4064663 | 3.4146671 | 0 | 0 |
| false | false | true | FF0FFF212 | 5.63478588 | 5.64835122 | 0.000054321999996886916 | 0.000054321999996886916 |
SELECT VERSION()
| version |
| ----------------------------------------------------------------------------------------------------------------- |
| PostgreSQL 13.9 on x86_64-apple-darwin19.6.0, compiled by Apple clang version 11.0.3 (clang-1103.0.32.62), 64-bit |
SELECT PostGIS_Version();
| postgis_version |
| ------------------------------------- |
| 3.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1 |
For convenience, here is the GeoJSON that I visualized above. It contains the lookup point and the two closest polygons (as returned by ST_AsGeoJSON):
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"properties": {},
"geometry": {
"type": "Point",
"coordinates": [
-118.46131,
34.31171
]
}
},
{
"type": "Feature",
"properties": {},
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
-118.590173754,
34.818203384
],
[
-118.554669477,
34.818248424
],
[
-118.519233938,
34.81833323
],
[
-118.500917157,
34.818263939
],
[
-118.484492998,
34.818255519
],
[
-118.466714818,
34.818356235
],
[
-118.431575103,
34.818707983
],
[
-118.413900556,
34.818924424
],
[
-118.396521556,
34.819207589
],
[
-118.375913258,
34.819322959
],
[
-118.326206757,
34.819675701
],
[
-118.298911338,
34.819887117
],
[
-118.298911338,
34.747356668
],
[
-118.30726467,
34.747324674
],
[
-118.311221449,
34.747259523
],
[
-118.315177432,
34.747194256
],
[
-118.316156934,
34.747178062
],
[
-118.31616929,
34.747315306
],
[
-118.316271488,
34.752065567
],
[
-118.312325229,
34.752140489
],
[
-118.312326535,
34.75219544
],
[
-118.312531114,
34.760874375
],
[
-118.311419271,
34.760887938
],
[
-118.311447254,
34.761794819
],
[
-118.307547133,
34.761876226
],
[
-118.307688072,
34.76914709
],
[
-118.312139537,
34.769051848
],
[
-118.312013231,
34.761782985
],
[
-118.320871559,
34.761596736
],
[
-118.320643006,
34.751982309
],
[
-118.32518382,
34.751892427
],
[
-118.325146684,
34.750645594
],
[
-118.320698938,
34.750725582
],
[
-118.320619613,
34.747238796
],
[
-118.320621472,
34.74711144
],
[
-118.319132498,
34.74712522
],
[
-118.31912928,
34.746989586
],
[
-118.318963759,
34.742294886
],
[
-118.324094459,
34.742214348
],
[
-118.324100254,
34.742351721
],
[
-118.324765599,
34.742340771
],
[
-118.324902701,
34.742340774
],
[
-118.324782787,
34.739791629
],
[
-118.324656399,
34.739792562
],
[
-118.321818215,
34.73983491
],
[
-118.321723743,
34.737421328
],
[
-118.318794936,
34.737463571
],
[
-118.318943061,
34.741690458
],
[
-118.315014788,
34.741751738
],
[
-118.314961133,
34.739937086
],
[
-118.311042888,
34.739996866
],
[
-118.311131734,
34.743628911
],
[
-118.307195019,
34.743692679
],
[
-118.307022176,
34.73292151
],
[
-118.307017236,
34.732701729
],
[
-118.306961869,
34.729197417
],
[
-118.311327393,
34.729174956
],
[
-118.311346965,
34.732617423
],
[
-118.3113478,
34.732644898
],
[
-118.311351901,
34.732727164
],
[
-118.315711282,
34.73267014
],
[
-118.315713308,
34.732587934
],
[
-118.315701602,
34.729152316
],
[
-118.320071884,
34.729129528
],
[
-118.320062563,
34.725505708
],
[
-118.315688277,
34.725527134
],
[
-118.315674941,
34.721907411
],
[
-118.320054027,
34.721894246
],
[
-118.320044883,
34.718366592
],
[
-118.320044698,
34.71828416
],
[
-118.318949128,
34.718285882
],
[
-118.317852726,
34.718287568
],
[
-118.316757189,
34.718289242
],
[
-118.315661616,
34.718290466
],
[
-118.315591765,
34.711052422
],
[
-118.31177053,
34.711045178
],
[
-118.306851348,
34.711035697
],
[
-118.306807821,
34.703851184
],
[
-118.306807578,
34.703768345
],
[
-118.311163994,
34.703793557
],
[
-118.311164549,
34.703710593
],
[
-118.311153205,
34.700151308
],
[
-118.315521476,
34.700168694
],
[
-118.315520836,
34.698345595
],
[
-118.315520196,
34.696522689
],
[
-118.315519752,
34.694760206
],
[
-118.315519343,
34.69294554
],
[
-118.315518911,
34.691133595
],
[
-118.315518483,
34.689349154
],
[
-118.315518511,
34.689321676
],
[
-118.315518613,
34.689225465
],
[
-118.324240593,
34.6892215
],
[
-118.324199826,
34.681929424
],
[
-118.324181765,
34.678286336
],
[
-118.324013286,
34.678286206
],
[
-118.321990626,
34.678294268
],
[
-118.319802221,
34.678302293
],
[
-118.319789286,
34.676482105
],
[
-118.319776231,
34.674800764
],
[
-118.319777068,
34.674668997
],
[
-118.324156113,
34.674640133
],
[
-118.324305833,
34.667353292
],
[
-118.324299819,
34.665534142
],
[
-118.324196404,
34.665535338
],
[
-118.319819256,
34.66554772
],
[
-118.319836703,
34.661900053
],
[
-118.321427846,
34.661592475
],
[
-118.322854394,
34.661224449
],
[
-118.324355708,
34.661602373
],
[
-118.324379145,
34.660067608
],
[
-118.322436434,
34.660073928
],
[
-118.320078162,
34.66007227
],
[
-118.319016614,
34.659630555
],
[
-118.318265028,
34.659228426
],
[
-118.317204016,
34.65892139
],
[
-118.316696353,
34.658343507
],
[
-118.316433733,
34.657886214
],
[
-118.316399174,
34.657374451
],
[
-118.316233563,
34.656701388
],
[
-118.315872914,
34.656190468
],
[
-118.315317537,
34.655922512
],
[
-118.314716272,
34.655843026
],
[
-118.314152005,
34.656042868
],
[
-118.313658715,
34.656335886
],
[
-118.31325996,
34.656463014
],
[
-118.312825923,
34.656501405
],
[
-118.312240655,
34.656360416
],
[
-118.311823758,
34.656137014
],
[
-118.311295459,
34.655640322
],
[
-118.310486705,
34.655074199
],
[
-118.309612031,
34.654525776
],
[
-118.308409437,
34.654003693
],
[
-118.307518342,
34.653679193
],
[
-118.306495295,
34.653471832
],
[
-118.306491282,
34.652731529
],
[
-118.306449169,
34.645426711
],
[
-118.304251649,
34.645432722
],
[
-118.302054126,
34.645438694
],
[
-118.299868246,
34.645445971
],
[
-118.298911338,
34.645448524
],
[
-118.298911338,
34.311655678
],
[
-118.590173754,
34.311655678
],
[
-118.590173754,
34.818203384
]
]
]
}
},
{
"type": "Feature",
"properties": {},
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
-118.625390854,
34.05061708
],
[
-118.625390854,
34.311655678
],
[
-118.456022326,
34.311655678
],
[
-118.456022326,
34.05061708
],
[
-118.625390854,
34.05061708
]
]
]
}
}
]
}
Solution
Even if it were some great circle discrepancy, I'd expect the point to be intersecting one of the polygons, not outside both of them.
This is not always a valid assumption. This would work if the polygons had the same edge with exactly the same end points. But this is not the case, so each follows its own great circle. Here is the exaggerated version of what probably happens here, using two Geographies that would be perfectly aligned as Geometries (I used BigQuery GeoViz to draw this):
select 1 x, st_geogfromtext('polygon((20 20, 40 20, 40 40, 20 40, 20 20))') g1
union all
select 2, st_geogfromtext('polygon((25 40, 50 40, 50 60, 25 60, 25 40))')
As you can see, there is both overlap and a hole between two "rectangles" when the edges follow great circle.
Basically, you need to use the type (Geometry vs Geography) that really represents the shape you have. PostGIS does the right calculations, but when you cast to a different type, you may get wrong results.
If the reason you use Geometry is the lack of PostGIS ST_Subdivide for Geography type, consider making your own. I've created a ST_Subdivide version for BigQuery that works with Geography (BigQuery lacks Geometry type) - it should be possible to port it to PostGIS https://mentin.medium.com/subdivide-and-conquer-any-geometry-ca4f0a4b8491.
- Call a function for each row in select - Postgres
- How to remove last element from array?
- Prioritizing one row over another, if it exists
- Golang postgresql query with transactions and squirrel
- Docker Compose: The "Hh" variable is not set. Defaulting to a blank string
- Why does this multiple "DO" block script work in DBeaver but fail in pgAdmin?
- postgresql "keepalives" parameter
- How to make supabase reflect whatever is on my payload cms
- Create Postgis Polygon from Geometry Points with adding X metres to polygon boundaries
- Convert string back to array of tuples
- How to connect GraphQL and PostgreSQL
- Why does updating a table in PostgreSQL need disk storage?
- PostgreSQL splitting time range into days
- How to concatenate strings of a string field in a PostgreSQL 'group by' query?
- Difference between CURRENT_TIMESTAMP and CURRENT_DATE
- How to return a table by rowtype in PL/pgSQL
- Allow docker container to connect to a local/host postgres database
- postgres. How to check if a field contains at least one substring from a list of substings?
- SQL ASC and DESC ordering on same column based on another column
- ANY and ALL against a small data set producing confusing results
- Generating a sequence number in a PostgreSQL database - concurrency & isolation levels
- How to get all the categories and sub categories in PostgreSQL?
- Only one child record is being created when trying to create a record with one to many relationship
- The entity type 'XmlSchemaCompilationSettings' requires a primary key to be defined. If you intended to use a keyless entity type
- How to use a PostgreSQL custom type as a nested POJO in jdbi?
- Store query result in a PL/pgSQL variable
- PostgreSQL drop constraint with unknown name
- How to match these 2 tables?
- PostgreSQL error : 22007 invalid_datetime_format
- How to use PGPASS file in Powershell to avoid password prompt?