How to Download Roads and Rivers for a District

This guide has 2 routes. Pick one, or read both. Each route is complete on its own.

Roads and rivers are the lines that most maps need. Roads tell you where people can reach, and water lines tell you where water flows and where it floods. Both come from OpenStreetMap, and How to download OpenStreetMap data explains how OSM tags work. This guide uses the same ideas for lines, for a whole district, and ends with the lengths by type.

The district is Ludhiana, Punjab, and the boundary is ludhiana_district.geojson from the shared data of this section. Download the zip and unzip it.

The tags

Roads have the key highway. The value says how important the road is, from motorway down to track. For a district map, the main classes are motorway, trunk, primary, secondary and tertiary. residential and service roads make up most of the length in a city and most of the size of the download.

Water lines have the key waterway. The values you will use are river, stream, canal, drain and ditch. In Punjab many of the water lines are canals and drains built by people, and the rivers are fewer, such as the Sutlej. Look at all the values before you pick. In the result below, canals and drains are more than four times the length of the rivers, so a query for river alone would leave most of the water lines out.

Route A: QuickOSM in QGIS

Open Vector > QuickOSM > QuickOSM. If QuickOSM is not installed, the first steps of the OpenStreetMap guide show how.
On the Quick query tab, type the Key highway and the Value primary. For several values, click the plus button to add a row, and add secondary and tertiary.
Under In, choose Layer extent and pick the layer ludhiana_district. This keeps the search inside the box of the district.
Click Run query. QuickOSM adds the roads as a line layer.
For the water, run a second query with the Key waterway, and the Values river, canal and drain. Leave the value empty to see every waterway, and then decide.
QuickOSM takes everything in the box of the district, and the box is larger than the district. To cut the lines to the boundary, open Vector > Geoprocessing Tools > Clip. Set Input layer to the roads, Overlay layer to ludhiana_district, and click Run.
To measure length, open the attribute table, click Open field calculator, make a new decimal field, and enter $length / 1000. $length uses the ellipsoid of the project, and the project should be in a metric CRS or set to measure in metres.

The steps use the names of the QuickOSM form as they are in QGIS 3.44. They change a little from version to version.

Route B: Python

You need Python with requests, geopandas and shapely.

import geopandas as gpd
import requests
from shapely.geometry import LineString

# Step 1: the district boundary, and its box. Overpass wants south, west, north, east.
district = gpd.read_file("ludhiana_workflow_data/ludhiana_district.geojson")
west, south, east, north = district.total_bounds
print("box:", round(south, 2), round(west, 2), round(north, 2), round(east, 2))
print("district area:", round(float(district.to_crs(32643).area.iloc[0] / 1e6)), "km2")

OVERPASS = "https://overpass-api.de/api/interpreter"
HEADERS = {"User-Agent": "timeandspatial-guide/1.0 (write your email address here)"}

def overpass(query, tries=3):
    # Step 2: send the query. The service is shared, so retry when it answers with something else than JSON.
    for _ in range(tries):
        reply = requests.post(OVERPASS, data={"data": query}, headers=HEADERS, timeout=180)
        if reply.ok and reply.text.lstrip().startswith("{"):
            return reply.json()["elements"]
    reply.raise_for_status()
    raise RuntimeError("Overpass did not return JSON")

def lines_from(ways, tag):
    # Step 3: turn each way into a line. A way carries its vertices as a list of points.
    return gpd.GeoDataFrame(
        [{"osm_id": w["id"], tag: w["tags"][tag], "name": w["tags"].get("name")} for w in ways],
        geometry=[LineString([(p["lon"], p["lat"]) for p in w["geometry"]]) for w in ways],
        crs="EPSG:4326",
    )

box: 30.56 75.36 31.02 76.33
district area: 3697 km2

# Step 1: main roads in the box. "out geom" returns the coordinates of every vertex.
box = f"({south},{west},{north},{east})"
road_classes = "motorway|trunk|primary|secondary|tertiary"
ways = overpass(f'[out:json][timeout:120]; way["highway"~"^({road_classes})$"]{box}; out geom;')
roads = lines_from(ways, "highway")
print(len(roads), "road pieces in the box")

# Step 2: water lines. Ask for river, canal, drain and stream.
ways = overpass(f'[out:json][timeout:120]; way["waterway"~"^(river|canal|drain|stream)$"]{box}; out geom;')
water = lines_from(ways, "waterway")
print(len(water), "water pieces in the box")

4298 road pieces in the box
241 water pieces in the box

# Step 1: cut the lines to the district. The box is bigger than the district, so some lines lie outside it.
roads_in = gpd.clip(roads, district)
water_in = gpd.clip(water, district)
print("roads:", len(roads), "pieces in the box,", len(roads_in), "after the clip")
print("water:", len(water), "pieces in the box,", len(water_in), "after the clip")

# Step 2: lengths in kilometres. Lengths need metres, so go to the UTM zone first.
for layer in (roads_in, water_in):
    layer["length_km"] = layer.to_crs(32643).length / 1000

# Step 3: length by type
print(roads_in.groupby("highway")["length_km"].sum().round(0).sort_values(ascending=False).to_string())
print(water_in.groupby("waterway")["length_km"].sum().round(0).sort_values(ascending=False).to_string())

roads: 4298 pieces in the box, 3739 after the clip
water: 241 pieces in the box, 179 after the clip
highway
tertiary     2772.0
secondary     409.0
trunk         363.0
primary       215.0
waterway
canal     288.0
drain     191.0
river     114.0
stream     78.0

# Step 1: density, kilometres of line for each square kilometre of the district
area_km2 = float(district.to_crs(32643).area.iloc[0] / 1e6)
print("road density:", round(roads_in["length_km"].sum() / area_km2, 2), "km per km2")
print("water line density:", round(water_in["length_km"].sum() / area_km2, 2), "km per km2")

# Step 2: how many water lines have a name?
named = water_in["name"].notna()
print("water lines with a name:", int(named.sum()), "of", len(water_in))
print(water_in.loc[named, "name"].value_counts().head(5).to_string())

# Step 3: save both layers into one GeoPackage
roads_in.to_file("ludhiana_osm_lines.gpkg", layer="roads", driver="GPKG")
water_in.to_file("ludhiana_osm_lines.gpkg", layer="water", driver="GPKG")

road density: 1.02 km per km2
water line density: 0.18 km per km2
water lines with a name: 52 of 179
name
Sua               23
Sutlej             7
Sem                5
Sirhind Canal      4
Sidhwan Branch     4

For Ludhiana, the box holds 4298 road pieces and 241 water pieces, and the clip leaves 3739 and 179. The district has about 3,759 km of main roads: 2,772 km of tertiary, 409 km of secondary, 363 km of trunk and 215 km of primary. The water lines add up to 671 km, with 288 km of canal, 191 km of drain, 114 km of river and 78 km of stream. That is 1.02 km of main road and 0.18 km of mapped water line for each square kilometre of the district. Of the 179 water pieces, 52 have a name, and the commonest are the Sua, the Sutlej, the Sem and the Sirhind Canal.

Overpass answers for the box and not for the district, which is why the clip comes after. The clip keeps the part of each line inside the boundary, and cuts a line in two when it crosses. For this reason the number of pieces after the clip can differ from the number before it.

What OSM lines can and cannot tell you

A road line is a centre line. It has no width, and it does not tell you how many lanes the road has or what it is paved with.
The highway class is chosen by the mapper. Two similar roads can have different classes.
Water lines are drawn where mappers saw them. In Punjab the main canals are well mapped, and the small field channels are mostly missing. The length of water lines here is the length of what was mapped and not of what exists.
Seasonal streams and drains may be mapped as a line even when they are dry.
A river may be a line in one place and a polygon in another. Add natural=water and waterway=riverbank to your query when you need the wide river surface.

Common mistakes

Querying waterway=river only, and concluding that the district has no water. Look at all the values first.
Asking for every highway in a large area. The answer is very large, and a request may time out. Start with the main classes.
Measuring length in degrees. A length from a latitude and longitude layer is in degrees. Reproject to metres first.
Forgetting to clip. The box is larger than the district.

The tags

Route A: QuickOSM in QGIS

Route B: Python

What OSM lines can and cannot tell you

Common mistakes

Key takeaways