Overpass API – How to Use Selector Keyword to Get All Water Bodies

openstreetmapoverpass-apipython

After installing non-default libraries in Google Colab:

%%shell
pip install OSMPythonTools
pip install geopandas

and importing useful libraries:

from OSMPythonTools.overpass import overpassQueryBuilder, Overpass
import numpy as np
from shapely.geometry import Polygon
import folium
import geopandas as gpd

I define a bounding box, within which I would like to find all water bodies. The bounding box is:

bbox = {'latLower':59.9,'lonLower':10.7,'latHigher':60.0,'lonHigher': 10.8}

I then use the overpass API to query data:

overpass = Overpass()
query = overpassQueryBuilder(bbox=[bbox['latLower'],bbox['lonLower'],bbox['latHigher'],bbox['lonHigher']],
                             elementType=['way', 'relation'], selector='"natural"="water"', includeGeometry=True)
result = overpass.query(query)

I also want to plot the bounding box, to check that I am on the right path, so I define a shapely polygon:

bbox_polygon_obj = Polygon(np.array([
                        [bbox['lonLower'],bbox['latLower']],
                        [bbox['lonLower'],bbox['latHigher']],
                        [bbox['lonHigher'],bbox['latHigher']],
                        [bbox['lonHigher'],bbox['latLower']]
                        ]))

Then I plot the results:

m = folium.Map(location=[0.5*(bbox['latLower']+bbox['latHigher']),
                         0.5*(bbox['lonLower']+bbox['lonHigher'])], zoom_start=12, tiles='CartoDB positron')

for each in result.ways():

    coord_pairs = np.array(each.geometry()['coordinates'][0])
    polygon_obj = Polygon(coord_pairs)

    geo_j = folium.GeoJson(data=gpd.GeoSeries(polygon_obj).to_json(),style_function=lambda x: {'fillColor': 'orange'})
    geo_j.add_to(m)

geo_j_bbox = folium.GeoJson(data=gpd.GeoSeries(bbox_polygon_obj).to_json(),style_function=lambda x: {'fillColor': 'red'})
geo_j_bbox.add_to(m)

m

I get this map:

enter image description here

As visible, quite a lot of water bodies are marked in the bounding box, but not all (ie see the big lake towards the north of the red area).

I would like to obtain all water bodies within the bouning box. I think what I did wrong is the selector='"natural"="water"' above, when I defined query.

How should I modify the defintition of query line above so I get all water bodies within the
bounding box?

Ready to run code is here.

Best Answer

I think you are experiencing one of the main issues with OSM and that is inconsistent tagging of data. I don't use the notebook style of coding you are using so not particularly knowledgeable about the syntax you are using, but a quick play on https://overpass-turbo.eu/ and I created the expression:

/*
This has been generated by the overpass-turbo wizard.
The original search was:
“(natural=water and type=multipolygon) or natural=water”
*/
[out:json][timeout:25];
// gather results
(
  // query part for: “natural=water and type=multipolygon”
  node["natural"="water"]["type"="multipolygon"]({{bbox}});
  way["natural"="water"]["type"="multipolygon"]({{bbox}});
  relation["natural"="water"]["type"="multipolygon"]({{bbox}});
  // query part for: “natural=water”
  node["natural"="water"]({{bbox}});
  way["natural"="water"]({{bbox}});
  relation["natural"="water"]({{bbox}});
);
// print results
out body;
>;
out skel qt;

This extracted out:

Results

Note the two examples I point to with red arrows, both lakes but one comes out because it has a multipolygon relationship, the other is an area tagged water.

Related Question