5

We've run into an odd query result with PostGIS. We have a large polygon and several smaller ones, and we're interested in checking whether or not the smaller ones are contained in the large polygon.
A visual inspection reveals they are indeed within the polygon (zoomed-out here so that the bottom boundary of the large polygon is visible; the smaller ones are in blue):
Zoomed-out comparison of polygons

However, running ST_Within(boundaries, geometry) results in False for all these polygons. After some introspection, ST_Within(ST_Transform(boundaries, 4326), ST_Transform(geometry, 4326)) results in True.
Both geometries are in local CRS (EPSG:3067), and both tables' entries in geometry_columns have srid=3067, same geometry type (GEOMETRY), same coord_dimension (2), etc.

Since this is resolved by transforming the coordinates, I can only assume that postgis does some partial transformation behind the scenes, but I'm not sure how to resolve this without transforming all our geometries to 4326.

For reference, here's a EWKT excerpt with the geometries:

with bbox as (select st_geomfromewkt('SRID=3067;POLYGON((64562.78588868573 6649100.87912667,205761.97180753836 7794679.579243775,674306.9889334934 7782478.683824445,757812.7263923198 6632497.631782339,64562.78588868573 6649100.87912667))') geom),
     polys as (select st_geomfromewkt('SRID=3067;POLYGON((284346.675638346 6640055.333759279,284346.46427643043 6640053.187765737,284345.83831320365 6640051.1242415225,284344.8218040813 6640049.222486716,284343.45381293906 6640047.555584686,284341.7869109092 6640046.187593544,284339.885156102 6640045.171084422,284337.82163188816 6640044.545121195,284335.675638346 6640044.333759279,284333.52964480384 6640044.545121195,284331.46612059 6640045.171084422,284329.5643657828 6640046.187593544,284327.89746375295 6640047.555584686,284326.5294726107 6640049.222486716,284325.51296348835 6640051.1242415225,284324.8870002616 6640053.187765737,284324.675638346 6640055.333759279,284324.8870002616 6640057.479752821,284325.51296348835 6640059.5432770355,284326.5294726107 6640061.445031842,284327.89746375295 6640063.111933872,284329.5643657828 6640064.479925014,284331.46612059 6640065.496434136,284333.52964480384 6640066.122397363,284335.675638346 6640066.333759279,284337.82163188816 6640066.122397363,284339.885156102 6640065.496434136,284341.7869109092 6640064.479925014,284343.45381293906 6640063.111933872,284344.8218040813 6640061.445031842,284345.83831320365 6640059.5432770355,284346.46427643043 6640057.479752821,284346.675638346 6640055.333759279))') geom)
select st_within(polys.geom, bbox.geom), st_within(st_transform(polys.geom, 4326), st_transform(bbox.geom, 4326))
from bbox, polys;

st_within|st_within|
---------+---------+
false    |true     |

@geozelot pointed out the distortion when projecting from a geography to a geometry - I did not anticipate it cause such a large discrepancy (~14km at the point of interest here). The obvious answer is to modify our original bounding box to be larger.

Below is a visual inspection as suggested; in red - the bounding box (large polygon) in WGS84, in purple - the projected bounding box in EPSG:3067, and in pink dots - the polygons of interest (in EPSG:3067).

Difference between projections

Improve this question
7
  • What is in your defined default in the spatial_ref_sys ? see gis.stackexchange.com/a/23477/276
    – Mapperz ♦
    Commented Apr 3 at 16:19
  • 3
    In fact, your results are is correct: the projected geometries do not intersect at all. Use a client software capable of plotting geometries in the scope of their projections (like QGIS) to get a visual validation. EPSG:4326 is not a projection - rather, geometries are said to be deprojected if given in a geographic coordinate reference - and calculations/predicates using them on a planar scale will be subject to distortions.
    – geozelot
    Commented Apr 3 at 16:48
  • 1
    Is your big polygon on the order of half the earth?
    – Ian Turton
    Commented Apr 3 at 17:02
  • 1
    @geozelot You're right! I was visualizing via dbeaver and it did not capture this difference. Our bounding box (the large polygon) starts out "deprojected" and we then project it to the target CRS (3067, in this case). This causes a distortion like you said, about 14km in difference. That's quite a large distortion - is that expected? I've edited the original post to reflect this finding.
    – hyit
    Commented Apr 3 at 18:25
  • 2
    Densify/segmentize the huge geometry before doing any re-projections.
    – user30184
    Commented Apr 3 at 18:58

1 Answer 1

Reset to default
7

This is a common issue when re-projecting huge polygons between geographic and projected coordinate systems when the geometries have only few vertices. The projection algorithms re-project just the vertices which hides the fact that the correct result contains curved segments.

See what happens if the EPSG:3067 geometry is densified with 10 km intervals

select st_transform(
st_segmentize(
st_geomfromewkt('SRID=3067;POLYGON((64562.78588868573 6649100.87912667,205761.97180753836 7794679.579243775,674306.9889334934 7782478.683824445,757812.7263923198 6632497.631782339,64562.78588868573 6649100.87912667))'),10000)
,4326);

The curved line with densified vertices shows the bbox approximately right. The green area presents fake area that is formed when the bbox has vertices only at the corners.

enter image description here

The reference geometry falls inside this fake area. The blue box presents the place of the geometry.

enter image description here

I do not know what makes the error in your map image but probably it is a same kind projection error of a huge rectangle with few vertices. Just have a look at the north coordinates of the EPSG:3067 geometries and you can see that the polys is closer to the equator than the bbox.

Improve this answer
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.