2

I have two tables A and B, both have valid geometries when I check with ST_IsValid. In my query, I intersect the two and the query stops in the middle because of this error: ERROR: lwgeom_intersection: GEOS Error: TopologyException: Input geom 0 is invalid: Self-intersection at or near point

I read ST_Union fails with TopologyException despite valid polygons and using ST_SnapToGrid which seems very similar, and tried the ST_CollectionExtract as mentioned in the accepted answer, but I still the same error.

When I run the query on a subset of the data I notice some records are ST_MultiPolygon while others are ST_Polygon, is this causing the issue and, if so, how can I fix it?

The original query has other fields which I've left out here but it generally looks like this:

SELECT 
   ST_CollectionExtract(ST_Intersection(ST_Transform(r.geom, 5070), 
                                        ST_Transform(c.geom, 5070)),3) as geom,
   a.place,
   b.name

FROM table1 as a 
JOIN table2 as b
WHERE a.geom && b.geom
AND ST_Intersects(ST_Transform(r.geom, 5070), ST_Transform(c.geom, 5070))

I am projecting it because I also calculate areas in the full query.

  • 2
    Reprojecting geometries can harm self-integrity, even without intersection. – Vince Nov 8 '19 at 0:58
  • 1
    I'm not saying to avoid projecting the data, just that some projections corrupt some geometries. I am also suggesting that focusing the question on the place where the error occurs will make it easier to solve. BTW: You can compute a geodesic area by casting to geography (if not using a buggy PostGIS) -- the units will be meters. – Vince Nov 8 '19 at 2:24
  • 3
    just because the two input geometries are valid doesn't imply that the intersection is a valid geometry. – Ian Turton Nov 8 '19 at 8:57
  • 2
    note that you wouldn't need to transform the filter geometries (ST_Intersects). you can also create the ST_Intersection without projection, and, as @Vince pointed at, cast (::) the input to ST_Area to GEOGRAPHY to receive precise measurements with meter as unit. both would increase performance significantly, as PostgreSQL can lookup the spatial index on the initial tables (transformed geometries are not covered by those indexes, so lookup will be sequential and slow). note that, technically, this still doesn't guarantee a valid intersection geometry, but reduces the chance for them. – ThingumaBob Nov 8 '19 at 9:08
  • 1
    As noted above reprojection can cause a valid polygonal geometry to become invalid. In this case ST_MakeValid could be used to ensure the transformed polygons are valid. – dr_jts Nov 18 '19 at 0:55
2

To resume:

The reprojection seems to be the problem here. If your need is to have an area in meters, the simplest way, as you have 4326 geom, would be to use the geography type:

Add an indexes on your geoms casted as a geography (to speed up the intersection search):

CREATE INDEX idx_table_1_geog ON table_1 USING gist (CAST(geom AS geography));
CREATE INDEX idx_table_2_geog ON table_2 USING gist (CAST(geom AS geography));
ANALYZE table_1; ANALYZE table_2;

Change your request to use geography (also you don't need bounding box comparison, ST_Intersects automatically use it):

SELECT 
        ST_CollectionExtract(ST_Intersection(a.geom::geography, b.geom::geography),3) as geom,
        a.place,
        b.name
    FROM table1 as a, table2 as b
    WHERE ST_Intersects(a.geom::geography, b.geom::geography);

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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