python - 地理装箱 - 地理边界内的平均值
python - geo binning - averaging values within a geo boundary
使用如下数据,- 捕获各个邻近位置的测量值
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网格的中心,以便我们得到如下结果。
Lat Long Mean_val Sample_count
虽然提议的方法可能很幼稚,但也想知道是否可以有一种基于 pandas
的方法
0.1米*0.1米面积平均数据的简单解法
为此,您必须将纬度、经度坐标转换为 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))
对之前解决方案的批评
正如我之前指出的那样,为了解决这个问题,我们必须将纬度、经度转换为 x、y 坐标。
在前面的解决方案中,我将经纬度转换为 utm 坐标。 utm 系统是一种制图投影,它将地球划分为 120 个区域:60 个北区和 60 个南区。所以当我们这样做时:
x,y,area_number,NS = utm.from_latlon(raw['Lat'],raw['Long'])
(x,y)
是我们在(area_number,NS)
区域的位置。我们可以得出结论,当且仅当我们的传感器位于同一 UTM 区域时,我们的解决方案才有效。
我们也可以使用直接将纬度、经度转换为 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)
使用如下数据,- 捕获各个邻近位置的测量值
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网格的中心,以便我们得到如下结果。
Lat Long Mean_val Sample_count
虽然提议的方法可能很幼稚,但也想知道是否可以有一种基于 pandas
0.1米*0.1米面积平均数据的简单解法
为此,您必须将纬度、经度坐标转换为 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))
对之前解决方案的批评
正如我之前指出的那样,为了解决这个问题,我们必须将纬度、经度转换为 x、y 坐标。
在前面的解决方案中,我将经纬度转换为 utm 坐标。 utm 系统是一种制图投影,它将地球划分为 120 个区域:60 个北区和 60 个南区。所以当我们这样做时:
x,y,area_number,NS = utm.from_latlon(raw['Lat'],raw['Long'])
(x,y)
是我们在(area_number,NS)
区域的位置。我们可以得出结论,当且仅当我们的传感器位于同一 UTM 区域时,我们的解决方案才有效。
我们也可以使用直接将纬度、经度转换为 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)