PostGIS: Wege ausserhalb von Landuse-Polygonen ermitteln

Rubrik: 

Für meine Velocarte wollte ich einen Overlay-Layer erstellen, in dem Straßen, die sich wegen ihres geringen Verkehrsaufkommens besonders für Alltagsradler eignen, farblich hervorgehoben werden. In einem ersten Schritt habe ich dafür alle Unclassified- und Service-Wege selektiert. Es hat sich aber gezeigt, dass darunter auch viele stark befahrene Straßen im innerstädtischen Bereich und in Industrie- und Gewerbegebieten waren. Um diese Straßen auszuschließen, habe ich zunächst die PostGIS-Funktion ST_WITHIN verwendet:

select osm_id,name,ref,highway,tracktype,way 
      from planet_osm_line as l 
      where highway in ('unclassified','service') 
            and not exists (select * from planet_osm_polygon as p
                            where p.landuse in ('residential','industrial','commercial') 
                                  and st_isvalid(p.way) and st_within(l.way,p.way));

Da ST_Within auf vollständiges Enthaltensein prüft, werden von dieser Abfrage auch Straßen erfasst, die aus dem Umland in ein Wohn- oder Gewerbegebiet hineinführen. An dieser Stelle bin ich auf die Funktion ST_Difference gestoßen, welche, angewendet auf einen Weg und ein Polygon, den Teil des Wegs liefert, der nicht im Polygon liegt. Der erste Ansatz war:

select l.osm_id,l.name,l.ref,l.highway,l.tracktype,
       st_difference(l.way,p.way) 
       from (select * from planet_osm_line 
                      where highway in ('unclassified','track')) as l, 
            (select way from planet_osm_polygon 
                    where landuse in ('residential','industrial','commercial')) as p 
       where st_isvalid(p.way) and st_intersects(l.way,p.way);

Das funktioniert wie erwartet bei Wegen, die sich mit genau einem der betroffenen Polygone schneiden. Schneidet ein Weg mehrere Polygone liefert die Abfrage mehrere Schnitte, wie man an folgendem Beispiel sehen kann:

Bild

Der Weg 5172636 (blau) verläuft durch ein Residential- und ein Commercial-Gebiet. ST_Difference in der obigen SQL-Abfrage liefert zwei Wege (rot): Weg 216 als Ergebnis des Schnitts mit dem Commercial-Gebiet und Weg 217 als Ergebnis des Schnitts mit dem Residential-Gebiet. Das Teilstück rechts oben, das in keinem der beiden Gebiete liegt, ist in beiden Wegen enthalten.

Meine Lösung ist etwas aufwändig, aber sie funktioniert. Hier der SQL-Code:

create table small_roads (id serial PRIMARY KEY,osm_id bigint,name text,ref text,highway text,tracktype text,way geometry);
drop table small_roads2;
create table small_roads2 (id serial PRIMARY KEY,osm_id bigint,name text,ref text,highway text,tracktype text,way geometry);

insert into small_roads (osm_id, name, ref,highway,tracktype,way) 
 (select l.osm_id,l.name,l.ref,l.highway,l.tracktype,st_difference(l.way,p.way) 
 from (select * from planet_osm_line where highway in ('unclassified','track')) as l, 
 (select way from planet_osm_polygon where landuse in ('commercial')) as p 
 where st_isvalid(p.way) and st_intersects(l.way,p.way));
 
insert into small_roads (osm_id, name, ref,highway,tracktype,way) 
   (select l.osm_id,l.name,l.ref,l.highway,l.tracktype,l.way 
      from planet_osm_line as l 
      where highway in ('unclassified','track') 
            and not exists (select * from planet_osm_polygon as p where p.landuse in ('commercial') 
                                     and st_isvalid(p.way) and st_intersects(l.way,p.way)));

insert into small_roads2 (osm_id, name, ref,highway,tracktype,way) 
 (select l.osm_id,l.name,l.ref,l.highway,l.tracktype,st_difference(l.way,p.way) 
 from (select * from small_roads) as l, 
 (select way from planet_osm_polygon where landuse in ('residential')) as p 
 where st_isvalid(p.way) and st_intersects(l.way,p.way));

insert into small_roads2 (osm_id, name, ref,highway,tracktype,way) 
   (select l.osm_id,l.name,l.ref,l.highway,l.tracktype,l.way 
      from small_roads as l 
      where not exists (select * from planet_osm_polygon as p where p.landuse in ('residential') 
            and st_isvalid(p.way) and st_intersects(l.way,p.way)));

delete from small_roads;

insert into small_roads (osm_id, name, ref,highway,tracktype,way) 
 (select l.osm_id,l.name,l.ref,l.highway,l.tracktype,st_difference(l.way,p.way) 
 from (select * from small_roads2) as l, 
 (select way from planet_osm_polygon where landuse in ('industrial')) as p 
 where st_isvalid(p.way) and st_intersects(l.way,p.way));

insert into small_roads (osm_id, name, ref,highway,tracktype,way) 
   (select l.osm_id,l.name,l.ref,l.highway,l.tracktype,l.way 
      from small_roads2 as l 
      where not exists (select * from planet_osm_polygon as p where p.landuse in ('industrial') 
            and st_isvalid(p.way) and st_intersects(l.way,p.way)));

Hier das Ergebnis für obiges Beispiel:

Bild

Die Tabelle small_roads enthält am Ende nur noch das rote Teilstück des Wegs, welches in keinem der beiden Gebiete liegt. Das ganze kann man entsprechend für weitere Landuse-Typen durchführen.