1

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.

2
  • Overpass API doesn't have a concept of a selector. That's something one of your libraries you're using has invented (probably OSMPythonTools). I'd suggest to try an Overpass QL query in overpass-turbo.eu, and then figure out how to translate that into something your wrapper libraries can understand.
    – mmd
    Sep 20 at 15:49
  • Thank you, I didn't know that website before, but seems very useful!
    – zabop
    Sep 20 at 16:12
1

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.

3
  • Thank you! Works like a charm; I'm going to leave the question open for a while in case someone wants to add something, if not, then I'll accept this answer.
    – zabop
    Sep 20 at 16:43
  • I got stuck a bit with implementing this; related question: gis.stackexchange.com/questions/412116/…
    – zabop
    Sep 20 at 19:12
  • 1
    Solved! :)
    – zabop
    Sep 20 at 21:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.