We have to filter out the roads and ditches without removing streams that cross roads or follow them closely. I’m going to use PostGIS to find the intersection of the streams lines data with a buffered roads polygon. If the intersected line is less than 50% of the length of the stream line, then the line is considered a road and removed from the streams dataset.
Buffering the roads
I found a subset of osm_roads that covered the whole area. Using QGIS, I performed a simple 25 foot buffer on the whole dataset. Instead of dissolving the boundaries in QGIS, however, I did in in PostGIS below using ST_Union().
Finding the intersection
Once the shapefiles were imported to postgres, I made a table that contained the ‘streams’ that were actually roads. This being a large dataset, I wanted to make sure a few things were in order: (1) that the tables I was using were indexed with a spatial index (using
gist()), and (2) that I wasn’t working with a multi-polygon layer. I used ST_Dump() to ensure I was working with single polygons in the roads layer.
drop table if exists big_roads; --create table of single polygon create table big_roads as select (ST_Dump(ST_Union(geom))).geom as geom from odm_roads_all; -- primary key alter table big_roads add Column gid serial primary key; -- create index create index idx_big_roads_geom on big_roads using gist(geom); -- new table drop table if exists intersected_streams; create table intersected_streams ( stream_id integer, stream_length numeric, geom geometry ); -- create intersected streams insert into intersected_streams (stream_id, geom) ( with clipped_streams as ( select * from ( select DISTINCT ON (streams.id) streams.id as stream_id, ST_Length(streams.geom) as stream_length, ST_Intersection(streams.geom, roads.geom) as geom from big_roads as roads INNER JOIN streams ON ST_Intersects(roads.geom, streams.geom)) as clipped ), int_streams as ( select * from clipped_streams where ST_Length(geom) > 0.25 * stream_length -- AND clipped_streams.stream_id IN (select id from streams) ) select id, geom from streams where id in (select stream_id from int_streams) );
After about 45 minutes, here is the result of my efforts:
(The image above actually has a roads polygon that is more accurate than a buffered roads)
Does it capture all the streams being misclassified as roads? For the most part it seems, yes. Does it remove all the other errors, such as ditches, wide highways, and parking lots? No. The biggest takeaway for me doing this project has been that urban areas add so much complexity to all of my problems. Send me back to Vermont!!!
You can see above the stream enters a culvert and passes under the road, and that really confused the algorithm.
My next steps are to instead of creating a new table of intersected streams, add a new “isroad” column to the existing stream dataset and reclassify “roads”. I want to make sure I preserve the connectivity of my network for the next step in the process of brook trout habitat analysis.