The steps you need to perform to complete this recipe are as follows:
- First, perform a self-spatial join between your MultiLineString dataset and the PostGIS ST_Intersects function and find intersections in the join context with the ST_Intersection PostGIS function. The following is the basic query, resulting in 1,448 records being selected:
postgis_cookbook=# SELECT r1.gid AS gid1, r2.gid AS gid2,
ST_AsText(ST_Intersection(r1.the_geom, r2.the_geom)) AS the_geom FROM chp03.rivers r1 JOIN chp03.rivers r2 ON ST_Intersects(r1.the_geom, r2.the_geom) WHERE r1.gid != r2.gid;
- You may hastily assume that all of the intersections are single points, but this is not the case; if you check the geometry type of the geometric intersections using the ST_GeometryType function, you have three different cases of intersection, resulting in the following geometries:
- An ST_POINT geometry for a simple intersection between two linear geometries.
-
- An ST_MultiPoint geometry, if two linear geometries intersect each other at more points.
- An ST_GeometryCollection geometry in cases where the two MultiLineString objects intersect and share part of the line. In such a case, the geometry collection is composed of ST_Point and/or ST_Line geometries.
- You can check the different cases with a query, shown as follows:
postgis_cookbook=# SELECT COUNT(*),
ST_GeometryType(ST_Intersection(r1.the_geom, r2.the_geom))
AS geometry_type
FROM chp03.rivers r1
JOIN chp03.rivers r2
ON ST_Intersects(r1.the_geom, r2.the_geom)
WHERE r1.gid != r2.gid
GROUP BY geometry_type;

(3 rows)
- First, try to compute the intersection for just the first two cases (intersections composed of the ST_Point and ST_MultiPoint geometries). Just generate a table with the Point and MultiPoint geometries, excluding the records that have an intersection composed of a geometric collection. By executing the following commands, 1,444 of the 1,448 records are imported (the four records with geometry collections are ignored using the ST_GeometryType function):
postgis_cookbook=# CREATE TABLE chp03.intersections_simple AS
SELECT r1.gid AS gid1, r2.gid AS gid2,
ST_Multi(ST_Intersection(r1.the_geom,
r2.the_geom))::geometry(MultiPoint, 4326) AS the_geom
FROM chp03.rivers r1
JOIN chp03.rivers r2
ON ST_Intersects(r1.the_geom, r2.the_geom)
WHERE r1.gid != r2.gid
AND ST_GeometryType(ST_Intersection(r1.the_geom,
r2.the_geom)) != 'ST_GeometryCollection';
- In case you want to import the points from the geometry collection too (but just the points, ignoring the eventual linestrings), one way to go about it is by using the ST_CollectionExtract function in the context of a SELECTCASE PostgreSQL conditional statement; this way, you can import all the 1,448 intersections as follows:
postgis_cookbook=# CREATE TABLE chp03.intersections_all AS
SELECT gid1, gid2, the_geom::geometry(MultiPoint, 4326) FROM (
SELECT r1.gid AS gid1, r2.gid AS gid2,
CASE
WHEN ST_GeometryType(ST_Intersection(r1.the_geom,
r2.the_geom)) != 'ST_GeometryCollection' THEN
ST_Multi(ST_Intersection(r1.the_geom,
r2.the_geom))
ELSE ST_CollectionExtract(ST_Intersection(r1.the_geom,
r2.the_geom), 1)
END AS the_geom
FROM chp03.rivers r1
JOIN chp03.rivers r2
ON ST_Intersects(r1.the_geom, r2.the_geom)
WHERE r1.gid != r2.gid
) AS only_multipoints_geometries;
- You may see the difference between the two processes, counting the total number of points in each of the generated tables, as follows:
postgis_cookbook=# SELECT SUM(ST_NPoints(the_geom))
FROM chp03.intersections_simple; --2268 points per 1444 records postgis_cookbook=# SELECT SUM(ST_NPoints(the_geom))
FROM chp03.intersections_all; --2282 points per 1448 records
- In the following screenshot (taken from QGIS), you may notice the generated intersections generated by both approaches. In the case of the intersection_all layer, you will notice that some more intersections have been computed (in red):

Layers of river intersections visualized in QGIS