组 lat/lon ，如果它们出现在边界框内

Posted 2023-03-27

【中文标题】组 lat/lon ，如果它们出现在边界框内

【英文标题】：Group lat/lon , if they occur witin a bounding box

【发布时间】：2021-12-26 18:16:45

【问题描述】：

我有一个如下的数据框：

我创建了一个 10 公里乘 10 公里的边界框。

def create_space(lat, lon, s=10):
    n = (180/math.pi)*(500/6378137)*10
    return lat - n, lon - n, lat + n, lon + n

现在，我想看看我的数据框中是否有两行或多行（纬度和经度）在边界框内。如果任何纬度和经度落在边界框内，我想添加该事件。 例如，如果 index[9] 落在 index[0] 的边界框上，则出现次数为 6495+23，并且 index[9] 将被删除。

我知道了：

我试过了

step =0.1
to_bin = lambda x: np.floor(x / step) * step
df["latbin"] = df.lat.map(to_bin)
df["lonbin"] = df.lon.map(to_bin)
#groups = df.groupby(("latbin", "lonbin"))

但它并没有解决我的问题，我不知道如何进一步。

【问题讨论】：

这能回答你的问题吗？ Is a lat lon within a bounding box?

是的，它部分是。但我想检查边界框内是否有 lon 和 lat，如果是，则添加事件并合并行。

嗨 Nicky - 欢迎来到 ***！ Please do not upload images of code/errors when asking a question. - 图像不能被索引/搜索，不能被屏幕阅读器读取，容易损坏等。相反，在显示输入/输出时使用格式化的代码块。谢谢！

另外，你能描述一下为什么它没有解决你的问题吗？如果导致错误，请提供traceback。 groupby 方法看起来很合适，因此它可以帮助我们更多地了解问题所在。

【参考方案1】：

我的解决方案是创建一个geopandas.geoseries.GeoSeries，它代表您在地图上的边界框。然后在 Python 中有现有的工具来测试一个点是否在它里面。

但是由于我没有你的数据，这里我只是用最简单的例子来向你展示我的代码是如何工作的。

# import packages
import geopandas as gpd
from shapely.ops import cascaded_union
from shapely import geometry

# create a function to build the "geopandas.geoseries.GeoSeries" for your bounding box
def create_boundingbox(p1,p2,p3,p4): 
    '''provide four corner points as (lon,lat), 
       the order is bottom-left(p1), bottom-right(p2), top-right(p3),top-left(p4) on a map'''
    p1 = geometry.Point(p1)
    p2 = geometry.Point(p2)
    p3 = geometry.Point(p3)
    p4 = geometry.Point(p4)
    pointList = [p1, p2, p3, p4, p1]
    boundingbox = geometry.Polygon([[p.x, p.y] for p in pointList])
    boundingbox = gpd.GeoSeries(cascaded_union(boundingbox))
    return boundingbox

# use some simple points as an example

# construct your box
p1 = (0,0)
p2 = (1,0)
p3 = (1,1)
p4 = (0,1)

box1 = create_boundingbox(p1,p2,p3,p4)

# now test if these points are inside or not
p5 = geometry.Point(0.5,0.5)
p6 = geometry.Point(15,15)

print(box1.contains(p5)) # this is True
print(box1.contains(p6)) # this is False

【参考方案2】：

请注意，当使用 lon-lat 时，您是在暗示某种地球，并且一组点描述了球体表面上两个可能的球形多边形。在这里，我假设您想要较小的。如果您以逆时针方向对点进行排序，则可以使用此处描述的方法使用线性代数。

这会考虑到您绕极点或逆子午线的情况，并查看连接点的大圆（因此连接 (000, 45) 和 (100, 45) 的线不是 遵循 45 度平行线。因此，请仔细考虑您在正方形中的含义，因为正方形（我们通常认为它们的方式）不能很好地覆盖在球体上。

根据您的申请，Jeremy 的回答可能就足够了，但在某些情况下可能不够。

我在下面建议的方法是基于对I asked here 的问题的回答，它解释了我在这里实施的数学。

首先，您需要将点转换为矢量（可以使用单位球体）

import numpy as np
def deg2rad(theta):
    return(theta * np.pi / 180)
def get_xyz(p):
    lon = deg2rad(p[0])
    lat = deg2rad(p[1])
        
    return(np.array([np.cos(lat) * np.cos(lon),
                     np.cos(lat) * np.sin(lon),
                     np.sin(lat)]))

因此，为您的四个点创建四个向量来描述笛卡尔坐标中的位置。例如，

p1 = [ 170, -10]
p2 = [-170, -10]
p3 = [-170, 10]
p4 = [ 170, 10]
v1 = get_xyz(p1)
v2 = get_xyz(p2)
v3 = get_xyz(p3)
v4 = get_xyz(p4)

多边形的每一边是great circle 的一部分，垂直于平面的向量是n1 = np.cross(v1, v2)，等等：

n1 = np.cross(v1, v2)
n2 = np.cross(v2, v3)
n3 = np.cross(v3, v4)
n4 = np.cross(v4, v1)

如果某个点 v5 在 v1、v2、v3 和 v4 所描述的多边形内，则 v5 与每个 n 向量的点积将大于 1。（如果它在一个边缘边数等于 1。）

p5 = [180, 0]
v5 = get_xyz(p5)

in_poly = (np.dot(v5, n1)) > 0 and \
    (np.dot(v5, n2)) > 0 and \
    (np.dot(v5, n3)) > 0 and \
    (np.dot(v5, n4)) > 0

print(in_poly) # True

p5 = [180, 20]
v5 = get_xyz(p5) 

in_poly = (np.dot(v5, n1)) > 0 and \
    (np.dot(v5, n2)) > 0 and \
    (np.dot(v5, n3)) > 0 and \
    (np.dot(v5, n4)) > 0

print(in_poly) # False

编辑：

今天早上我意识到我没有解释如何对数据进行分类。我在这里把它放在字典里，但你可以在数据框中创建一个列来做类似的事情。另外，并不是说我的网格框需要是凸的（没有 pacman 形状），并且网格点的角需要按逆时针顺序提供。

def in_poly(poly, point):
   v_vec = [get_xyz(p) for p in poly]
   n_vec = [np.cross(v_vec[i], v_vec[(i+1)%len(v_vec)]) for i in range(len(v_vec))]
   v_p = get_xyz(point)
   dot_prod = [np.dot(n, v_p) for n in n_vec]
   if all(d > 0 for d in dot_prod):
      return True
   else:
      return False

p1 = [ 170, -10]
p2 = [-170, -10]
p3 = [-170,  10]
p4 = [ 170,  10]
p5 = [ 150, -10]
p6 = [ 150,  10]
P = 1: [p1, p2, p3, p4],
     2: [p1, p4, p6, p5]
InBox = 1: 0, 2: 0

Npts = 100
lonlim = [150, 190]
latlim = [-20, 20]
points = np.stack([np.random.randint(lonlim[0], lonlim[1], Npts),
                   np.random.randint(latlim[0], latlim[1], Npts)])


for i in range(Npts):
   for key in InBox.keys():
      if in_poly(P[key], points[:,i]):
         InBox[key] += 1

print(InBox) 
#1: 25, 2: 22