0

I am trying to get the intersection points between the boundary of a city (stored in the database as a MultiPolygon) and the OSM roads (stored as Multiline string).

I am using ST_intersect(city.geom, roads.geom) as a where clause but I am not getting any results.

If I open the tables in QGis I see that there are intersections. Do I need to transform any of the geom ?

3

1 Answer 1

3

Do you really need to work with MultiPolygon and MultiLineString ? If not, first you need to check if you multi geometries really contains multi geometries - i.e if spatially disconnected entities are present. I often work with OSM data and I also get multi geometries whenever I import QGIS-created shapefiles in PostGIS.

Run the following query for each of your tables :

SELECT count(
    CASE WHEN ST_NumGeometries(t.geom) > 1
    THEN 1 END
    ) AS multi,
    count(t.geom) AS total
FROM the_table t; --Replace the_table by the table you want to test

If the first column returns 0, it means you don't actually have multi geometries, they're only stored as such. Then you can easily convert your geometry column to simple geometry using :

ALTER TABLE multipolygon_table
    ALTER COLUMN geom TYPE geometry(Polygon, 4326) --Replace 4326 by the SRID you are using
    USING ST_GeometryN(geom, 1);

ALTER TABLE multilinestring_table
    ALTER COLUMN geom TYPE geometry(LineString, 4326) 
    USING ST_GeometryN(geom, 1);

If your data really contains multi geometries (first column of the first query returns a value > 0), you should dump them into simple geometries (Take a look here https://stackoverflow.com/questions/21719941/postgis-convert-multipolygon-to-single-polygon)

To get the intersection between the boundary of your city and the roads, I suggest you first create a CTE (kind of a temporary table created with WITH) containing the border of your polygons as line(s) and then run the ST_Intersection function combined with ST_Intersects (to work only on roads that actually intersect the boundaries) to obtain points :

WITH boundary AS (
    SELECT id, ST_ExteriorRing(geom) as geom
    FROM city_table
    )
SELECT r.id as road_id, ST_Intersection(r.geom, b.geom) as geom
FROM road_table r, boundary b
WHERE ST_Intersects(r.geom, b.geom);

This is suppose to work because the intersection between two lines returns a point. Be aware that ST_ExteriorRing only works with simple polygons so you need to convert them in the first time.

Forgive me if some explanations are obvious but I tried to be as explicit as possible.

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.