[GIS] Intersection postgis TopologyException

postgistopology

I'm trying to make a quite simple analysis in Postgis, but I'm finding some problems. I want to in the first place create buffers for then, intersects themselves but one, this will be differenced. With this product I want to intersects again with another layer. But, here is the problem Postgis complains about a non-existent error topology. I say non-existent because there is no feature in that coordinate.

> Error performing intersection: TopologyException: Input geom 1 is
> invalid: Self-intersection at or near point 723120.23091093975
> 4379700.195452339 at 723120.23091093975 4379700.195452339

Here is the SQL sentence that I'm using:

     select 1 as id, st_intersection(diff_1.st_difference, mapa_ruido_dia.geom)
    from
    (select 1 as id, st_difference(intersectionEMT_valenbisi_carrilBici.st_intersection, buf_ocio.st_union)
    from
    (select 1 as id, st_intersection(buf_carril_bici.st_union,inter_emt_valenbisi.st_intersection)
    from
    (select 1 as id, st_intersection(buf_emt.st_union,buf_valenbisi.st_union)
    from
    (SELECT 1 as id, ST_union(st_buffer(emt_paradas.geom, 200)) FROM emt_paradas )as buf_emt,
    (SELECT 1 as id, ST_union(st_buffer(valenbisi.geom, 200)) FROM valenbisi )as buf_valenbisi) as inter_emt_valenbisi,
    (SELECT 1 as id, ST_union(st_buffer(carril_bici.geom, 200)) FROM carril_bici )as buf_carril_bici) as intersectionEMT_valenbisi_carrilBici,
    (SELECT 1 as id, ST_union(st_buffer(ocio.geom, 250)) FROM ocio )as buf_ocio)as diff_1, mapa_ruido_dia 
     where mapa_ruido_dia.gridcode='1' or mapa_ruido_dia.gridcode='2'

Here is the database

This is the place where postgis says there's a problem,and the two layers I want to intersect

The Yellow star is where I have the topology exception.

EDIT:

What I try to do is the next:

First I want to create a buffer for emt_paradas, and the layer valenbisi, and then intersect them, the result I will call "intersect_1".
Second, I want to create a buffer for carril_bici, and intersect with "intersect_1". The result I will call "intersect_2".
Third, I want to create a buffer for layer ocio and make a difference with "intersect_2". The result I will call: "diff_1".
Finnaly I want to intersect "diff_1" with layer mapa_ruido_dia, ideally where the field gridcode equals one and two.

Best Answer

I just had a similar problem, and the solution is the following:

You have many queries, so the first question is which one of them fails. I'd suggest bisecting the entire SQL query to find where the problem happens exactly. Comment one half of the query (but make it still work) and see if the problem is still there, if so, bisect this half, if not, the error is in the other, bisect that half instead.

When you see which particular query fails, check which table has invalid geometries. My example:

 with resi as (select st_transform(way, 3857) way from residential),
      road_buff as (select st_buffer(way, 40) way from roads)
 select st_intersection(r.way, b.way) from resi r, road_buff b
 where st_intersects(r.way, b.way)

This query failed, and I just looked for which CTE contains invalid geometry:

 with resi as (select st_transform(way, 3857) way from residential),
      road_buff as (select st_buffer(way, 40) way from roads)
 select r.way, st_isvalid(r.way) val from resi r

Then the same for road_buff CTE. Turned out resi had invalid geometries. Then I filtered by val=false, looked at that table in QGIS, and saw that some polygons were drawn self-intersecting.

The solution was to add st_makevalid(way) on that CTE:

    with resi as (select st_transform(st_makevalid(way), 3857) way from residential

And this solved the issue.

Here's a good tool to fix different invalid geometries.

Related Solutions

PostGIS Polygons – How to Separate Polygons Based on Intersection

Since you said you get a group of intersecting polygons for each polygon you're interested in, you may want to create what is referred to as a "polygon overlay".

This isn't exactly what Adam's solution is doing. To see the difference, take a look at this picture of an ABC intersection:

ABC intersection

I believe Adam's solution will create an "AB" polygon that covers both the area of "AB!C" and "ABC", as well as an "AC" polygon that covers "AC!B" and "ABC", and a "BC" polygon that is "BC!A" and "ABC". So the "AB", "AC", and "BC" output polygons would all overlap the "ABC" area.

A polygon overlay produces non-overlapping polygons, so AB!C would be one polygon and ABC would be one polygon.

Creating a polygon overlay in PostGIS is actually pretty straightforward.

There are basically three steps.

Step 1 is extract the linework [Note that I'm using the exterior ring of the polygon, it does get a little more complicated if you want to correctly handle holes]:

SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines

Step 2 is to "node" the linework (produce a node at every intersection). Some libraries like JTS have "Noder" classes you can use to do this, but in PostGIS the ST_Union function does it for you:

SELECT ST_Union(the_geom) AS the_geom FROM (...your lines...) AS noded_lines

Step 3 is to create all the possible non-overlapping polygons that can come from all those lines, done by the ST_Polygonize function:

SELECT ST_Polygonize(the_geom) AS the_geom FROM (...your noded lines...)

You could save the output of each of those steps into a temp table, or you can combine them all into a single statement:

CREATE TABLE my_poly_overlay AS
SELECT geom FROM ST_Dump((
    SELECT ST_Polygonize(the_geom) AS the_geom FROM (
        SELECT ST_Union(the_geom) AS the_geom FROM (
            SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines
        ) AS noded_lines
    )
)

I'm using ST_Dump because the output of ST_Polygonize is a geometry collection, and it is (usually) more convenient to have a table where each row is one of the polygons that makes up the polygon overlay.

[GIS] PostGIS intersection function

Clip and Intersect are two similar, but different operations, but that's not your problem here. It looks to me like your Intersect worked fine, the problem is that the shape you're Intersecting with, the red frame, is a line object and not a polygon. So, the only road segments being output are the ones that directly touch that line. If you want to clip the roads to just those segments inside the polygon, make the red frame shape a polygon rather than a line feature and try again. I don't have PostGIS experience so I can't help with that part but if your red frame is a line like you said, making it a polygon should solve the problem.

Best Answer

Related Solutions

PostGIS Polygons – How to Separate Polygons Based on Intersection

[GIS] PostGIS intersection function

Related Question