如何转换Overpass Turbo请求为原生请求并解决Python调用API异常？

问题背景 #

我需要根据区域名称获取其边界轮廓并导入GeoPandas。
以下请求在Overpass Turbo中可正常运行：

[out:json][timeout:25];
{{geocodeArea:Велижский район}}[type="boundary"]->.a;
(
    relation["boundary"="administrative"]["admin_level"="6"](area.a);
);

(
  way(r);
  node(w);
);
out skel qt;

另一个请求也可正常运行：

[out:json];
relation(3342303)[type="boundary"];
out geom;

我可将结果下载为GeoJSON并通过GeoPandas绘图：

import geopandas as gpd

df = gpd.read_file('export.geojson')
df.geometry = df.geometry.to_crs(epsg=32636)
df.plot()

已尝试的方案 #

import overpass

api = overpass.API()

query = """
relation(3342303)[type="boundary"];
out geom;
"""
res = api.Get(query)

预期结果 #

获取JSON格式的返回数据

实际出现的错误 #

---------------------------------------------------------------------------
UnknownOverpassError                      Traceback (most recent call last)
<ipython-input-10-2eef0e720784> in <cell line: 9>()
      7 out geom;
      8 """
----> 9 res = api.Get(query)

1 frames
/usr/local/lib/python3.10/dist-packages/overpass/api.py in _as_geojson(self, elements)
    212                             polygons.append([points])
    213                         else:
--> 214                             raise UnknownOverpassError("Received corrupt data from Overpass (incomplete polygon).")
    215                 # Then get the inner polygons
    216                 for member in elem.get("members", []):

UnknownOverpassError: Received corrupt data from Overpass (incomplete polygon).

更新内容1 #

我成功获取数据到Python，但这些线条总面积为0，未解决问题：

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

def extract_geometries(element, geometries):
    if 'members' in element:
        members = element['members']
        for member in members:
            if member['type'] == 'way' and 'geometry' in member:
                geometry = member['geometry']
                coordinates = [(point['lon'], point['lat']) for point in geometry]
                line = LineString(coordinates)
                geometries.append(line)
            elif member['type'] == 'relation':
                extract_geometries(member, geometries)

def close_polygon(geometry):
    if geometry.geom_type == 'Polygon':
        coordinates = list(geometry.exterior.coords)
        if coordinates[0] != coordinates[-1]:
            coordinates.append(coordinates[0])
            return Polygon(coordinates)
    return geometry

def get_relation_geometry(relation_id):
    # Define the Overpass API URL
    url = "https://overpass-api.de/api/interpreter"

    # Define the Overpass query for the relation ID
    query = f"""
    [out:json];
    rel({relation_id});
    out geom;
    """

    # Send the GET request to the Overpass API
    response = requests.get(url, params={"data": query})

    # Check if the request was successful
    if response.status_code == 200:
        # Convert the response to GeoPandas dataframe
        data = response.json()

        # Extract the 'elements' list
        elements = data['elements']

        # Create an empty list to store geometries
        geometries = []

        # Iterate over the elements and extract geometries
        for element in elements:
            extract_geometries(element, geometries)

        # Create a GeoPandas dataframe from the geometries
        gdf = gpd.GeoDataFrame(geometry=geometries)

        # Forcefully close the polygons
        gdf['geometry'] = gdf['geometry'].apply(close_polygon)

        return gdf
    else:
        print(f"Request failed with status code {response.status_code}")
        return None

# Example usage
relation_id = 3342303
gdf = get_relation_geometry(relation_id)
if gdf is not None:
    print(gdf.head())

内容的提问来源于stack exchange，提问作者Kirill Setdekov