shp2pgsql -c -D -s 4269 -W "latin1" -I tl_2013_48201_roads.shp tlharris | psql -U researcher -d vcdb3
As well as the US national file for the coastlines: shp2pgsql -c -D -s SRID SHAPEFILE4269 -W "latin1" -I tl_2017_us_coastline.shp SCHEMA.TABLE tlcoastline | psql -h HOST U researcher -d DATABASE vcdb3 Unfortunately, this doesn't show lakes... You can get all lines from https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2019&layergroup=All+Lines But you have to do them county by county. shp2pgsql -c -D -s 4269 -W "latin1" -I tl_2019_50007_edges.shp tlchittenden | psql -U USERresearcher -d vcdb3 And there's so many features...
One problem with this method is that there are partial intersections. We could use ST_Difference to return the part that doesn't intersect... a bigger problem is that we are restricted to pairs of geometries. Using a cross product we could test all rows against all other rows. But then we'd need to aggregate the intersections... One method is to use a recursive CTE [https://gis.stackexchange.com/questions/269875/aggregate-version-of-st-intersection]. ST_Union is truly an aggregate function but not what we want in this context.