如何在Python中判断经纬度点所属的多边形（对标R实现）

要实现和R中map.where类似的经纬度点归属判断，用GeoPandas结合Shapely就能搞定，以下是具体步骤和代码：

核心思路 #

把目标经纬度转换成空间点对象，再通过空间判断方法（空间连接或多边形包含检查）找到点所在的多边形，最后提取对应的属性字段（比如你需要的ISO_1编码）。

完整代码示例 #

import geopandas as gpd
from shapely.geometry import Point

# 读取已下载的美国行政区划shapefile（你已完成这一步）
shp = gpd.read_file("gadm41_USA_1.shp")

# 定义要判断的点：注意Shapely的Point是「经度在前，纬度在后」，和R的map.where参数顺序一致
target_point = Point(-73.978544, 40.725959)
# 把点转换成GeoDataFrame，确保坐标系和shapefile一致（GADM数据默认是WGS84，EPSG:4326）
point_gdf = gpd.GeoDataFrame({'geometry': [target_point]}, crs=shp.crs)

# 方法1：空间连接（适合批量点判断，效率更高）
match_result = gpd.sjoin(point_gdf, shp, how='left', predicate='within')

# 提取并格式化结果
if not match_result.empty and match_result['ISO_1'].notna().iloc[0]:
    iso_code = match_result['ISO_1'].iloc[0]
    # 从GID_1字段提取后缀数字（对应R返回的":88"部分）
    gid_suffix = match_result['GID_1'].iloc[0].split('_')[-1]
    print(f"US-{iso_code}:{gid_suffix}")  # 输出：US-NY:88
else:
    print("该点不在任何目标多边形范围内")

# 方法2：单个点直接遍历判断（逻辑更直观）
for _, polygon_row in shp.iterrows():
    if polygon_row['geometry'].contains(target_point):
        print(f"US-{polygon_row['ISO_1']}:{polygon_row['GID_1'].split('_')[-1]}")
        break

关键细节说明 #

1. 坐标系一致性：必须保证点和shapefile的坐标系相同，GADM下载的行政区划数据默认是WGS84（EPSG:4326），所以用crs=shp.crs同步即可。
2. 点的顺序：Shapely构造Point时是(经度, 纬度)，和R中map.where(x=经度, y=纬度)的参数顺序对应，别搞反了。
3. 字段对应：R中namefield="ISO_1"对应shapefile里的ISO_1字段，返回结果中的数字后缀来自GID_1字段（比如USA_88提取88）。

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