python - 地理分箱 - 地理边界内的平均值

Posted 2023-02-23

【中文标题】python - 地理分箱 - 地理边界内的平均值

【英文标题】：python - geo binning - averaging values within a geo boundary

【发布时间】：2020-02-17 12:53:21

【问题描述】：

使用如下数据， - 捕获各个近距离位置的测量结果

Lat Long    val
35.611053   139.628525  -72.82
35.61105336 139.6285236 -78.04
35.61105373 139.6285223 -72.99
35.61105409 139.6285209 -69.04
35.61105445 139.6285195 -65.4
35.61105482 139.6285182 -66.68
35.61105518 139.6285168 -65.82
35.61105555 139.6285155 -64.47
35.61105591 139.6285141 -71.26
35.61105627 139.6285127 -68.36
35.61105664 139.6285114 -74.48
35.611057   139.62851   -74.27
35.61105736 139.62851   -77.97
35.61105773 139.62851   -68.66
35.61105809 139.62851   -70.21
35.61105845 139.62851   -76.05
35.61105882 139.62851   -88.83
35.61105918 139.62851   -73.17
35.61105955 139.62851   -67.63
35.61105991 139.62851   -71.85
35.61106027 139.62851   -77.42
35.61106064 139.62851   -71.08
35.611061   139.62851   -79.27

需要对该数据执行分箱操作——即每 0.1x0.1 米获取val 中所有值的平均值。一种方法可能是找到边缘（如 NW、SW、NE 和 SE）并将其划分为一组 0.1x0.1 米的网格，并在每个网格内查找值并计算平均值并归因于中心的纬度/经度网格，以便我们得到如下结果。

Lat Long    Mean_val    Sample_count

虽然提议的方法可能很幼稚，但也想知道是否有基于 pandas 的方法

【问题讨论】：

我熟悉 2D 直方图，但不熟悉 Lat 和 Long。每列的最大值和最小值是多少？另外：我们可以假设心脏在每个正方形中都是平坦的，对吗？

是的，为此我们可以认为地球是平的 :) 东北角将由 (max Lat, min Long) 定义，而东南角将是 (min Lat, max Long)

我明白了，但让我们概括一下：我们要做一个二维直方图。 x 将是 Lat 和 y 将是 Long。我们将从min(x) 到max(x) 和同上为y 制作垃圾箱。可能不需要它，但我不想给你一个错误的代码：min(x), min(y), max(x), max(y) 的值是什么？这就是二维网格所在的位置。

在地理坐标中，x 将是 Long，y 将是 Lat。所以NW 角落将是(min(x),max(y)。 NE 角落将是 max(x), max(y)。 SE 角落将是 min(x), min(y)。 SW 角落将是 max(x), min(y)。

真正的问题是我认为你如何计算你的平均价值。假设您在该区域中只有一个点，在 (35.6 139.6) 处的值为 -77。它的热量是如何散发的？我的意思是必须有一个功能。像高斯分布或其他东西。

【参考方案1】：

按 0.1 m * 0.1 m 面积平均数据的简单解决方案

为此，您必须将纬度、经度坐标转换为 x,y 坐标。

这里我使用utm 模块：

x,y,_,_ = utm.from_latlon(latitude, longitude) 

之后，您可以创建一个新列，以分米表示您的 x,y 坐标：

def apply_fun (raw):
    x,y,_,_ = utm.from_latlon(raw['Lat'],raw['Long']) 
    return str(np.round(x*10))+"|"+str(np.round(y*10))

然后将其添加到您的数据框：

x = df.apply(lambda row : apply_fun(row),axis=1)
df.insert(3,'Group',x)

然后你应用 groupby 函数：

gdf = df.groupby(['Group']).agg("Lat":["mean"],"Long":["mean","count"],"val":["mean"])
gdf = gdf.reset_index().drop(columns=['Group'],level=0)
gdf.columns = [' '.join(col) for col in gdf.columns]

我们完成了！ :)

上一个解决方案的推广

要按 k1 米 * k2 米的面积对数据进行分组，只需修改此函数即可：

def apply_fun (raw):
    x,y,_,_ = utm.from_latlon(raw['Lat'],raw['Long']) 
    return str(np.round(x/k1))+"|"+str(np.round(y/k2)) 

对之前解决方案的批评

正如我之前指出的，要解决这个问题，我们必须将 lat、long 转换为 x、y 坐标。

在之前的解决方案中，我将 lat,long 转换为 utm 坐标。 utm 系统是一种地图投影，将地球划分为 120 个区域：北纬 60 和南纬 60。所以当我们这样做时：

x,y,area_number,NS = utm.from_latlon(raw['Lat'],raw['Long'])

(x,y) 是我们在(area_number,NS) 区域的位置。我们可以得出结论，当且仅当我们的传感器位于同一 UTM 区域时，我们的解决方案才有效。

我们也可以使用直接将 lat,long 转换为 x,y 坐标的 ECEF 转换来进行这种转换。我不知道这些方法的精度，因为我们被要求精确到十分之一米，所以我更喜欢选择看起来更准确的 utm 转换。

如果您想使用这样的 ECEF 方法：

import pyproj
def gps_to_ecef_pyproj(lat, lon, alt):
    ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
    lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
    x, y, z = pyproj.transform(lla, ecef, lon, lat, alt, radians=False)

    return x, y, z

x,y,z = gps_to_ecef_pyproj(raw['Lat'],raw['Long'],0)

（我从这里获取代码：https://gis.stackexchange.com/questions/230160/converting-wgs84-to-ecef-in-python）

【讨论】：

能否提供分米换算的参考。并且更概括它如何通过n meters网格完成任何n meters - 例如0.44 * 0.44米或1.21 * 1.21米

还有一种方法可以得到df.groupby(['Group']).mean()为组平均的样本数。

0.1 米 = 1 分米 = 10 厘米。要获得平均样本数，您可以执行 len(df) - len(df.groupBy(['Group']).mean())。如果该答案对您有所帮助，请考虑验证答案^^

哦，好吧，我去编辑那个。我没看到这部分

如果没有关于数据的更多地理背景，这个答案通常是不正确的。从utm 的documentation 注意到函数utm.from_latlon 返回一个4 元组：(EASTING, NORTHING, ZONE NUMBER, ZONE LETTER)。如果数据跨越多个utm zones，那么以这种方式使用它们的 (x,y) 是错误的。