使用 geopandas 查找地理边界 N 米内的所有点

Question

我需要找到美国边境 N 米内的所有点（名称、纬度、经度）。我的做法是：

1. 制作一个包含多多边形的地理数据框，其中每个多边形都是一个缓冲边界。

# Get the U.S. borders
border = gpd.read_file('cb_2018_us_nation_5m')
border = border.to_crs("EPSG:32636")
# Get just the boundary
just_border = border.boundary
# Put a 10km buffer around the line. We now have multiple polygons (AK, HI, CONUS), each that look like a donut.
border_buffered = just_border.buffer(10000, cap_style=3, join_style=2) 

然后将我的数据框转换为地理数据框。

gdf = gpd.GeoDataFrame(temp, geometry=gpd.points_from_xy(temp.Longitude, temp.Latitude))
gdf = gdf.set_crs("EPSG:32636")

然后裁剪多边形内的点。

ngdf = gpd.clip(gdf, border_buffered, keep_geom_type=False)

这没有找到点：

ngdf.shape
(0, 15)

所以我尝试使用面具：

gdf.reset_index(inplace=True)
gdf.shape
(25746326, 17)
pip_mask = gdf.within(border_buffered.loc[0, 'geometry'])

Traceback (most recent call last):
  File "/Users/kovar/work/a50_utils/foo.py", line 45, in <module>
    pip_mask = gdf.within(border_buffered.loc[0, 'geometry'])
  File "/Users/kovar/miniforge3/envs/a50-dev/lib/python3.9/site-packages/pandas/core/indexing.py", line 889, in __getitem__
    return self._getitem_tuple(key)
  File "/Users/kovar/miniforge3/envs/a50-dev/lib/python3.9/site-packages/pandas/core/indexing.py", line 1063, in _getitem_tuple
    self._has_valid_tuple(tup)
  File "/Users/kovar/miniforge3/envs/a50-dev/lib/python3.9/site-packages/pandas/core/indexing.py", line 720, in _has_valid_tuple
    self._validate_key_length(key)
  File "/Users/kovar/miniforge3/envs/a50-dev/lib/python3.9/site-packages/pandas/core/indexing.py", line 761, in _validate_key_length
    raise IndexingError("Too many indexers")
pandas.core.indexing.IndexingError: Too many indexers

我猜我的多边形太多了，或者缓冲的线太复杂了。

Answer 1

• 您的问题清楚地说明了 N 米。这需要使用 UTM CRS 几何。美国边界穿过多个 UTM CRS 区域
• 首先将部分边界与 UTM 区域（线段）相关联
• 对每个线段使用适当的 UTM CRS 并添加缓冲区。这个例子创建了一个 50KM 的缓冲区。存在一些问题，因此如果需要，可以减少缓冲区的大小
• 现在有带有多边形的地理数据框，其中包括缓冲区，以及用于显示实际线串和实际使用的缓冲区的列

import geopandas as gpd
import shapely.geometry
import numpy as np
import plotly.express as px
import requests, io
from pathlib import Path
from zipfile import ZipFile
import urllib
import pandas as pd

# fmt: off
# US geometry
url = "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_nation_5m.zip"
f = Path.cwd().joinpath(urllib.parse.urlparse(url).path.split("/")[-1])

if not f.exists():
    r = requests.get(url, stream=True, headers={"User-Agent": "XY"})
    with open(f, "wb") as fd:
        for chunk in r.iter_content(chunk_size=128):
            fd.write(chunk)
    zfile = ZipFile(f)
    zfile.extractall(f.stem)

gdf2 = gpd.read_file(list(f.parent.joinpath(f.stem).glob("*.shp"))[0])
gdf2["geometry"] = gdf2["geometry"].apply(lambda g: g.boundary)

# get UTM zones
utm = gpd.GeoDataFrame.from_features(requests.get("https://opendata.arcgis.com/datasets/b294795270aa4fb3bd25286bf09edc51_0.geojson").json()).set_crs(gdf2.crs)
utm = utm.loc[~utm["ZONE"].eq(0)]
# fmt: on

# join US boundary to UTM zomes using spatial join, i.e. segments of boundary that pass through multiple zone
us_utm = gpd.sjoin(utm, gdf2)
us_utm = us_utm.merge(
    gdf2["geometry"], left_on="index_right", right_index=True, suffixes=("", "_right")
)

# create linestrings for each zone / boundary segment
us_utm_sects = us_utm["geometry"].intersection(gpd.GeoSeries(us_utm["geometry_right"]))

# function to create a buffer on line seqment
def buffer(g, buffer=1):
    gdf_sect = gpd.GeoDataFrame(geometry=[g], crs="EPSG:4326")
    utm_crs = gdf_sect.estimate_utm_crs()
    while True:
        g2 = gdf_sect.to_crs(utm_crs).buffer(buffer, cap_style=2).to_crs("EPSG:4326")
        tb = g2.total_bounds[[0,2]]
        bounds = abs(tb[0] - tb[1])
        # make sure epsg:4326 geometry makes sense, if not reduce size of buffer on this segment
        if bounds < 60:
            break
        buffer = buffer / 2
    return {"geometry": g2.values[0], "utm_crs": str(utm_crs), "line_geometry":g, "buffer":buffer}


# create geo data frame of us border with 50km buffer
us_buf = gpd.GeoDataFrame(us_utm_sects.apply(buffer, buffer=5 * 10 ** 4).apply(pd.Series))
us_buf["color"] = pd.factorize(us_buf["utm_crs"])[0]

使用它

• 已经使用世界范围内的地震来证明
• 在点和缓冲美国边界之间进行空间连接的简单案例
• 曾在靠近美国边境的墨西哥地震上空盘旋，以证明它有效

gdf_eq = gpd.GeoDataFrame.from_features(
    requests.get(
        "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_week.geojson"
    ).json()
)

# earthquakes near US border...
us_eq = gpd.sjoin(gdf_eq, us_buf)

px.scatter_mapbox(us_eq, lat=us_eq.geometry.y, lon=us_eq.geometry.x, hover_data=["title"]).update_layout(
    mapbox={
        "style": "carto-positron",
        "zoom":1,
        "center":{"lat":sum(us_eq.total_bounds[[0,2]])/2,"lon":sum(us_eq.total_bounds[[1,3]])/2},
        "layers": [
            {
                "source": us_buf.geometry.__geo_interface__,
                "type": "fill",
                "color": "yellow",
                "opacity": 0.2,
            }
        ],
    },
    margin={"l":0,"t":0,"b":0,"r":0}
    
)

只过滤墨西哥边境

• 简单地连接到墨西哥多边形

world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres")).set_crs("EPSG:4326")
us_buf = gpd.sjoin(us_buf, world.loc[world["iso_a3"].eq("MEX")]).drop(columns=["index_right"])

Answer 2

将迈克尔的评论推广到答案。

“看起来border_buffered是一个GeoSeries，而不是一个GeoDataFrame，因为GeoDataFrame.boundary只是返回一个GeoSeries，所以就做border_buffered.loc[0]”

• https://stackoverflow.com/users/3888719/michael-delgado

这解决了问题。