Behavior of "out geom"

[uswoods]

When I don't use out geom I'm getting an error (see below), but when I use it, the json output is duplicated. Am I utilizing the option in a wrong matter?

import json
import overpass

endpoint = "https://overpass-api.de/api/interpreter"
api = overpass.API(timeout=100, endpoint=endpoint)
way = '45667934'
query = 'way('+str(way)+');out geom;'
result = api.get(query, responseformat="geojson")
print(result)

Error:

UnboundLocalError: local variable 'geometry' referenced before assignment

[mvexel]

This example seems to run fine. Make sure you have the latest version of the module.

Just a few pointers if I may:

• If you want to set the output verbosity to something else than the default body, use the verbosity parameter of the get call. See the relevant documentation.
• You don't need to set the interpreter if you're sticking with the default.

[uswoods]

But sorry, it runs fine but has duplicate geometry:

Geometry I

{
  "features": [
    {
      "geometry": {
        "coordinates": [
          [
            4.2711979,
            51.3542058
          ],
          ...
        ],
        "type": "LineString"
      },
      "id": 45667934,
      "properties": {
        "addr:housenumber": "600",
        "addr:postcode": "2040",
        "addr:street": "Scheldelaan",
        "alt_name": "BASF",
        "landuse": "industrial",
        "man_made": "works",
        "name": "BASF Antwerpen nv",
        "website": "http://www.basf.be",
        "wikipedia": "nl:BASF (chemieconcern)"
      },
      "type": "Feature"
    },

Geometry II

  {
     "geometry": {
       "coordinates": [
         [
           4.2711979,
           51.3542058
         ],
         ...
       ],
       "type": "LineString"
     },
     "id": 45667934,
     "properties": {
       "addr:housenumber": "600",
       "addr:postcode": "2040",
       "addr:street": "Scheldelaan",
       "alt_name": "BASF",
       "landuse": "industrial",
       "man_made": "works",
       "name": "BASF Antwerpen nv",
       "website": "http://www.basf.be",
       "wikipedia": "nl:BASF (chemieconcern)"
     },
     "type": "Feature"
   }
 ],
 "type": "FeatureCollection"
}

[mvexel]

That is because you use out geom in your query instead of using verbosity='geom' in the get() call:

import overpass
api = overpass.API()
query = 'way(45667934);'
result = api.get(query, verbosity='geom')
print(result)

gets you

{"features": [{"geometry": {"coordinates": [[4.2711979, 51.3542058], [4.2690866, 51.3543414],... [4.287715, 51.353279], [4.2870876, 51.3533091], [4.2804398, 51.3536934], [4.2798607, 51.353716], [4.2747813, 51.3540099], [4.2711979, 51.3542058]], "type": "LineString"}, "id": 45667934, "properties": {"addr:housenumber": "600", "addr:postcode": "2040", "addr:street": "Scheldelaan", "alt_name": "BASF", "landuse": "industrial", "man_made": "works", "name": "BASF Antwerpen nv", "website": "http://www.basf.be", "wikipedia": "nl:BASF (chemieconcern)"}, "type": "Feature"}], "type": "FeatureCollection"}

I updated the code in the repl as well.

[mvexel]

Does it work for you now @uswoods ?

[uswoods]

Actually I would like to use the output for calculating the area of the way:

import overpass
import json
from functools import partial
import pyproj    
from shapely.ops import transform
from shapely.geometry import mapping,shape
from shapely.geometry.polygon import Polygon

equal_area_projection = "aea"

api = overpass.API()
query = 'way(45667934);'
result = api.get(query, verbosity='geom')
# will fail
geom = shape(result)

# Create partially applied transformation
transformation = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:4326'),
    pyproj.Proj(
        proj=equal_area_projection,
        lat1=geom.bounds[1],  # TODO: Different parameters?
        lat2=geom.bounds[3])
)

geom_area = transform(transformation, geom)

print(round(geom_area.area/10000))

But this doesn't work due to the FeatureCollection. Do you any trick for unwinding this in python?

[mvexel]

overpass will output a FeatureCollection, so if you want to access its individual geometries you will have to access its features dict:

import overpass
from shapely.geometry import shape, Polygon
api = overpass.API()
query = 'way(45667934);'
result = api.get(query, verbosity='geom')
linestring = shape(result.features[0].geometry)
poly = Polygon(linestring)
print(poly.area)

[uswoods]

Thanks for the solution. But there is something weird:

According to JOSMs measurement plugin BASF Antwerp has a size of 606 hectares. Shapely gives 609 after projecting. poly.area gives 0.0007859869400399928
Why ?

[uswoods]

@mvexel When I query for a multipolygon, the output FeatureCollection is empty:

import overpass
api = overpass.API()
query = 'relation(2688333);'
result = api.get(query, verbosity='geom')
print(result)

[mvexel]

poly.area will yield the area in unprojected units (square degrees).

Relations are not supported in overpass yet. See #48.

[uswoods]

@mvexel I've now patched my local overpass version that relations do indeed work thanks to @t-g-williams work.

The polygon itself could be read by Shapely, but it's complaining due to the FeatureCollection. How can I remove the FeatureCollection from the GeoJSON's "DOM tree"? (Where should I better ask this question? On StackOverflow it was deleted)