Skip to content

osm waterbodies

Jump to bottom
Tim Tisler edited this page Nov 1, 2016 · 10 revisions

Setting Up

Rasterizing OpenStreetMaps Water-body

  • Run ogr2ogr to create CSV files of rivers and lakes
# water lines
ogr2ogr -f "CSV" waterlines.csv planet-latest.osm.pbf -sql "SELECT OGR_GEOM_WKT from lines WHERE waterway='river' OR waterway='riverbank' OR waterway='canal' or natural='water'" -dialect OGRSQL

# water multilines
ogr2ogr -f "CSV" watermultilines.csv planet-latest.osm.pbf -sql "SELECT OGR_GEOM_WKT from multilinestrings WHERE waterway='river' OR waterway='riverbank' OR waterway='canal' or natural='water'" -dialect OGRSQL

# water polys
ogr2ogr -f "CSV" waterpolys.csv planet-latest.osm.pbf -sql "SELECT OGR_GEOM_WKT from multipolygons WHERE waterway='river' OR waterway='riverbank' OR waterway='canal' or natural='water'" -dialect OGRSQL

ogr2ogr -f "CSV" watercoastlines.csv coastline/water_polygons.shp -sql "select OGR_GEOM_WKT from water_polygons" -dialect OGRSQL

  • Merge the data files into one
tail -q -n +2 watercoastlines.csv waterlines.csv watermultilines.csv waterpolys.csv > osm-water.csv
  • You may want to randomize the vectors in the file, which will make the rasterization much more efficient. Basically, all the big ocean polygons are in one place in the file and it that one task will take a long time in rasterization
shuf -o osm-water-randomized.csv osm-water.csv
mv osm-water-randomized.csv osm-water.csv

Now use can use that file to rasterize water (remember to move it somewhere the mrgeo cluster can access it)

Remember MrGeo also needs a columns file (osm-water.csv.columns)

<?xml version="1.0" encoding="UTF-8"?>
<AllColumns firstLineHeader="false">
  <Column name="GEOMETRY" type="Nominal"/>
</AllColumns>

Run mrgeo and rasterize the vectors.

This command will rasterize the water at zoom level 12 (~30m)

mrgeo mapalgebra -o osm-waterbodies -e "rasterizevector([/mrgeo/vectors/osm-water.csv], 'MASK', '12z')"

Rasterizing OpenStreetMaps Roads

  • Run ogr2ogr to create CSV files of rivers and lakes
# roads lines
ogr2ogr -f "CSV" roads.csv planet-latest.osm.pbf -sql "SELECT OGR_GEOM_WKT,highway,maxspeed from lines WHERE highway='motorway' OR highway='trunk' OR highway='primary' OR highway='secondary' OR highway='tertiary' OR highway='unclassified' OR highway='residential' OR highway='service' OR highway='motorway_link' OR highway='trunk_link' OR highway='primary_link' OR highway='secondary_link' OR highway='tertiary_link' OR highway='living_street' OR highway='pedestrian' OR highway='track' OR highway='road'" -dialect OGRSQL

# roads muiltilines
ogr2ogr -f "CSV" roads-multilines.csv planet-latest.osm.pbf -sql "SELECT OGR_GEOM_WKT,highway,maxspeed from multilinestrings WHERE highway='motorway' OR highway='trunk' OR highway='primary' OR highway='secondary' OR highway='tertiary' OR highway='unclassified' OR highway='residential' OR highway='service' OR highway='motorway_link' OR highway='trunk_link' OR highway='primary_link' OR highway='secondary_link' OR highway='tertiary_link' OR highway='living_street' OR highway='pedestrian' OR highway='track' OR highway='road'" -dialect OGRSQL

# roads polys
ogr2ogr -f "CSV" roads-polys.csv roads-latest.osm.pbf -sql "SELECT OGR_GEOM_WKT,highway,maxspeed from multipolygons WHERE highway='motorway' OR highway='trunk' OR highway='primary' OR highway='secondary' OR highway='tertiary' OR highway='unclassified' OR highway='residential' OR highway='service' OR highway='motorway_link' OR highway='trunk_link' OR highway='primary_link' OR highway='secondary_link' OR highway='tertiary_link' OR highway='living_street' OR highway='pedestrian' OR highway='track' OR highway='road'" -dialect OGRSQL
  • Merge the data files into one
tail -q -n +2 roads-lines.csv roads-multilines.csv roads-polys.csv > osm-roads.csv
  • Make sure each road has a correct maxspeed value. Do this by running this python script.
#!/usr/bin/python
import csv


default_speeds = {
'motorway':100,
'trunk':70,
'primary':65,
'secondary':60,
'tertiary':50,
'unclassified':30,
'residential':30,
'service':20,
'motorway_link':70,
'trunk_link':65,
'primary_link':60,
'secondary_link':50,
'tertiary_link':40,
'living_street':10,
'pedestrian':6,
'track':20,
'road':35}

country_speeds = {
'at:urban':50,
'at:rural':100,
'at:trunk':100,
'at:motorway':130,
'ch:urban':50,
'ch:rural':80,
'ch:trunk':100,
'ch:motorway':120,
'cz:urban':50,
'cz:rural':90,
'cz:trunk':130,
'cz:motorway':130,
'dk:urban':50,
'dk:rural':80,
'dk:motorway':130,
'de:living_street':7,
'de:urban':50,
'de:walk':7,
'de:rural':100,
'fi:urban':50,
'fi:rural':80,
'fi:trunk':100,
'fi:motorway':120,
'fr:urban':50,
'fr:rural':90,
'fr:trunk':110,
'fr:motorway':130,
'hu:urban':50,
'hu:rural':90,
'hu:trunk':110,
'hu:motorway':130,
'it:urban':50,
'it:rural':90,
'it:trunk':110,
'it:motorway':130,
'jp:national':60,
'jp:motorway':100,
'ro:urban':50,
'ro:rural':90,
'ro:trunk':100,
'ro:motorway':130,
'ru:living_street':20,
'ru:rural':90,
'ru:urban':60,
'ru:motorway':110,
'sk:urban':50,
'sk:rural':90,
'sk:trunk':130,
'sk:motorway':130,
'si:urban':50,
'si:rural':90,
'si:trunk':110,
'si:motorway':130,
'se:urban':50,
'se:rural':70,
'se:trunk':90,
'se:motorway':110,
'gb:nsl_single':96.54,
'gb:nsl_dual':112.63,
'gb:motorway':112.63,
'ua:urban':60,
'ua:rural':90,
'ua:trunk':110,
'ua:motorway':130,
'living_street':6}

words = ['maxspeed', 'ambiguous', 'signals', 'none', 'walk', 'variable', 'national', 'fixme', 'unposted', 'implicit']
def is_number(s):
    try:
        float(s)
        return True
    except ValueError:
        return False

def default_speed(highway):
    if not highway in default_speeds:
        return default_speeds['road']
    else:
        return default_speeds[highway]

def get_maxspeed(speed, units, highway):
    speeds = speed.split(';|,-')
    maxspeed = 0
    for sp in speeds:
        sp = sp.replace(units, '')
        if (is_number(sp)):
            if units == 'mph':
                sp = 1.609 * float(sp) 
            elif units == 'knots':
                sp = 1.852 * float(knots)
            else:
                sp = float(sp)
                
            if sp > maxspeed:
                maxspeed = sp
    if maxspeed > 0:
        speed = maxspeed
    else:
        speed = default_speed(highway)

    return speed

file = open('osm-roads.csv')
roads = csv.reader(file)

output = open('osm-roads-corrected.csv', 'w')
roads_corrected = csv.writer(output)

row = 0
for road in roads:
    if row > 0:
        if (row % 100000) == 0:
            print row
            
        geom = road[0]
        highway = road[1]
        speed = road[2].lower().strip()
        
        # if we don't have a speed, give it a default
        if len(speed) == 0:
            speed = default_speed(highway)
        elif not is_number(speed):
            if 'kph' in speed:
                speed = get_maxspeed(speed, 'kph', highway)
            elif 'km/h' in speed:
                speed = get_maxspeed(speed, 'km/h', highway)
            elif 'kmh' in speed:
                speed = get_maxspeed(speed, 'kmh', highway)
            elif 'mph' in speed:
                speed = get_maxspeed(speed, 'mph', highway)
            elif 'knots' in speed:
                speed = get_maxspeed(speed, 'knots', highway)
            elif speed in country_speeds:
                speed = country_speeds[speed]
            elif speed in words:
                speed = default_speed(highway)
            else:
                speed = get_maxspeed(speed, '', highway)
        
        road[2] = 3.6 / float(speed)  # kph to spm
    
    roads_corrected.writerow(road)
    row = row + 1

print row

Now use can use that file to rasterize water (remember to move it somewhere the mrgeo cluster can access it)

Remember MrGeo also needs a columns file (osm-roads.csv.columns)

<?xml version="1.0" encoding="UTF-8"?>
<AllColumns firstLineHeader="false">
  <Column name="GEOMETRY" type="Nominal"/>
  <Column name="highway" type="Nominal"/>
  <Column name="maxspeed" type="Numeric"/>
</AllColumns>

Run mrgeo and rasterize the vectors.

This command will rasterize the roads at zoom level 12 (~30m)

mrgeo mapalgebra -o osm-roads -e "rasterizevector([/mrgeo/vectors/osm-roads.csv], 'MASK', '12z')"
Clone this wiki locally