如何不使用单一经纬度获取指定城市周边10km范围内的城镇？

核心思路 #

放弃单一经纬度的点半径搜索，转而先获取城市完整的多边形边界，生成向外扩展10公里的缓冲区域，再搜索该区域内的城镇，同时排除原城市范围内的结果。

具体步骤 #

1. 获取目标城市的多边形边界数据

主流地图/地理数据API都能返回城市的GeoJSON格式边界，关键是开启几何信息返回：

• 以Google Maps Geocoding API为例，请求时指定result_type为对应城市的行政层级（比如administrative_area_level_1或administrative_area_level_2），并添加include_geometry=true参数，就能拿到包含Polygon或MultiPolygon的几何数据。
• 用OpenStreetMap Nominatim的话，请求format=geojson并设置polygon_geojson=1，同样能获取边界。

2. 生成10公里缓冲区域

经纬度坐标系（WGS84）受地球曲率影响，不适合直接做距离缓冲，需先转换为平面投影（如UTM，根据城市所在纬度选对应带号）再生成缓冲：

• Python示例（依赖shapely和pyproj库）：

from shapely.geometry import shape
from pyproj import Transformer
import json

# 假设已获取城市的GeoJSON边界字符串
city_geo = json.loads("城市边界GeoJSON内容")
city_poly = shape(city_geo["geometry"])

# 转换为UTM投影（示例为EPSG:32633，对应北纬30-36度东半球）
proj_to_utm = Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)
city_poly_utm = city_poly.transform(proj_to_utm.transform)

# 生成10公里缓冲（单位：米）
buffer_poly_utm = city_poly_utm.buffer(10000)

# 转换回WGS84经纬度
proj_to_wgs84 = Transformer.from_crs("EPSG:32633", "EPSG:4326", always_xy=True)
buffer_poly = buffer_poly_utm.transform(proj_to_wgs84.transform)

• JavaScript示例（依赖@turf/turf库）：

import * as turf from '@turf/turf';

// 城市GeoJSON边界对象
const cityGeo = 城市边界GeoJSON对象;
// 生成10公里缓冲（turf内部自动处理投影转换）
const bufferGeo = turf.buffer(cityGeo, 10, {units: 'kilometers'});

3. 查询缓冲区域内的城镇

使用支持多边形过滤的API搜索目标要素：

• Google Maps Places API：调用textsearch接口，指定type=locality（对应城镇级地点），并通过polygon参数传入缓冲区域的经纬度坐标列表。
• OpenStreetMap Overpass API：编写查询语句搜索缓冲区域内的城镇，示例：

[out:json][timeout:25];
(
  node["place"="town"](poly:"缓冲区域经纬度串");
  node["place"="village"](poly:"缓冲区域经纬度串");
);
out body;
>;
out skel qt;

4. 过滤结果

最后排除原城市边界内的城镇，只保留缓冲区域内、原城市范围外的结果：

• 用地理空间库（如shapely）判断每个城镇的坐标点是否满足：不在原城市多边形内，但在缓冲多边形内。

注意事项 #

• 免费API可能存在请求次数限制或不支持复杂多边形搜索，优先考虑OpenStreetMap相关工具，数据开源且灵活。
• UTM投影带号需根据城市位置选择，避免跨带误差；也可使用turf这类自适应投影的库简化操作。

内容的提问来源于stack exchange，提问作者J.D.