Skip to content

Untested

Zach Levitt edited this page Sep 15, 2022 · 1 revision

untested from Tim Wallace: get OSM data for raster boundary

Make bounding box geojson for clipping

tif_name="${1%.tif}" gdaltindex $tif_name".geojson" $1

Get projection from tif

tif_proj=$(gdalinfo -proj4 $1 | sed -n '/proj=/,/Origin/p' | awk -F'Origin' '{print $1}'); echo Image projection: $tif_proj

Get latlon bbox for OSM query

mapshaper -i $tif_name".geojson" -proj from="$tif_proj" "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" -o "wgs_"$tif_name".geojson" latlonbounds=$(mapshaper -i "wgs_"$tif_name".geojson" -rectangle -info 2>&1 | awk 'NR==7' | awk -F' ' '{print $2}') echo LatLon Bounds: $latlonbounds

Query OSM for that area

curl --location --globoff "http://www.overpass-api.de/api/xapi?way[building=*][bbox="$latlonbounds"]" -o "buildings"$tif_name".osm" curl --location --globoff "http://www.overpass-api.de/api/xapi?way[highway=*][bbox="$latlonbounds"]" -o "roads"$tif_name".osm"

Convert all road lines from the .osm file

ogr2ogr -f "ESRI Shapefile" $tif_name"roads.shp" "roads"$tif_name".osm" -skipfailures -sql "SELECT * FROM lines"

Clip. This assumes the projection of the tif is webmercator. Change the -proj flag if not.

mapshaper-xl -i $tif_name"_roads.shp" encoding=UTF-8 -proj "$tif_proj" -clip $tif_name.geojson -o $tif_name"_roads.shp" force

Convert all polygon buildings from the .osm file

ogr2ogr -f "ESRI Shapefile" $tif_name"buildings.shp" "buildings"$tif_name".osm" -skipfailures -sql "SELECT * FROM multipolygons"

Clip. This assumes the projection of the tif is webmercator. Change the -proj flag if not.

mapshaper-xl -i $tif_name"_buildings.shp" encoding=UTF-8 -proj "$tif_proj" -clip $tif_name.geojson -o $tif_name"_buildings.shp" force