使用值列从 GPS 数据创建插值多边形

Question

我有一些 GPS 数据，其中每个点都有一个值（就像空气质量一样）。我可以绘制这些点，例如用大叶，并将值映射到圆的大小，如下所示：

import pandas, numpy, folium
lat = numpy.random.uniform(45, 45.01, 250)
lon = numpy.random.uniform(3, 3.02, 250)
value = numpy.random.uniform(0,50,250)
df = pandas.DataFrame({'lat': lat, 'lon': lon, 'value': value})
mymap = folium.Map(location = [lat.mean(), lon.mean()], tiles="OpenStreetMap", zoom_start=14)
for elt in list(zip(df.lat, df.lon, df.value)):
    folium.Circle(elt[:2], color="blue", radius=elt[2]).add_to(mymap)
mymap.save('mymap.html')

我想处理这些数据，对其进行插值并创建一个 shapefile 作为输出，其中插值的多边形包含平均值并显示高值和低值区域（右侧的假图片）。当然，多边形的限制将通过插值自动生成。

我怎样才能实现这个目标？因为我尝试过在 folium 中使用 HeatMap 工具，但它的设计目的是插值点的密度，而不是与每个点关联的值！我希望它不会太复杂。谢谢，注意：我使用 folium，但我对其他 python 库没问题。

Answer 1

我们可以将这些值作为权重提供给 folium 中的 HeatMap（请参阅data文档中的参数）。

这是一个例子。我将方框下半部分的值除以 5，以证明尽管数据点密度相对均匀，但热图在图片的上半部分清楚地显示了更高的幅度：

import pandas as pd
import numpy as np
import folium
from folium.plugins import HeatMap
import matplotlib as mpl

# parameters
n = 250                     # number of points
lat0 = 40.7                 # coordinates will be generated uniformly with
lon0 = -73.9                # lat0 - eps <= lat < lat0 + eps
eps = 0.1                   # lon0 - eps <= lon < lon0 + eps
v_min, v_max = 0, 100       # min, max values

# generating values
lat = np.random.uniform(lat0 - eps, lat0 + eps, n)
lon = np.random.uniform(lon0 - eps, lon0 + eps, n)
value = numpy.random.uniform(v_min, v_max, n)
df = pandas.DataFrame({'lat': lat, 'lon': lon, 'value': value})

# to demonstrate the effect of weights on the heatmap,
# we'll divide values below the center of the box by K = 5
K = 5
df.loc[df['lat'] < lat0, 'value'] /= K

# plotting the map, both the points themselves and the heatmap
m = folium.Map(location = [lat0, lon0], tiles="OpenStreetMap",
               zoom_start=11, width=400, height=400)
for elt in list(zip(df.lat, df.lon, df.value)):
    folium.Circle(elt[:2], color="white", radius=elt[2]).add_to(m)

# df.values used here is a (250, 3) numpy.ndarray
# with (lat, lon, weight) for each data point
HeatMap(data=df.values, min_opacity=0.1).add_to(m)

m

输出：

更新：

这是一种稍微不同的方法，不使用内置的 HeatMap（正如 @agenis 提到的，鉴于https://github.com/python-visualization/folium/issues/1271 ，这可能是目前更好的选择）。

我们首先转换数据，制作方格网格，其中每个方格的值是该单元格内原始数据点值的平均值（因此它不取决于该单元格内的点数，只取决于它们的值）。

然后我们可以通过在地图上绘制 GeoJson 多边形来可视化这些正方形。这是一个例子：

# define the size of the square
step = 0.02

# calculate values for the grid
x = df.copy()
x['lat'] = np.floor(x['lat'] / step) * step
x['lon'] = np.floor(x['lon'] / step) * step
x = x.groupby(['lat', 'lon'])['value'].mean()
x /= x.max()
x = x.reset_index()

# geo_json returns a single square
def geo_json(lat, lon, value, step):
    cmap = mpl.cm.RdBu
    return {
      "type": "FeatureCollection",
      "features": [
        {
          "type": "Feature",
          "properties": {
            'color': 'white',
            'weight': 1,
            'fillColor': mpl.colors.to_hex(cmap(value)),
            'fillOpacity': 0.5,
          },
          "geometry": {
            "type": "Polygon",
            "coordinates": [[
                [lon, lat],
                [lon, lat + step],
                [lon + step, lat + step],
                [lon + step, lat],
                [lon, lat],
              ]]}}]}


# generating a map...
m = folium.Map(location=[lat0, lon0], zoom_start=11, width=400, height=400)

# ...with squares...
for _, xi in x.iterrows():
    folium.GeoJson(geo_json(xi['lat'], xi['lon'], xi['value'], step),
                   lambda x: x['properties']).add_to(m)

# ...and the original points
for elt in list(zip(df.lat, df.lon, df.value)):
    folium.Circle(elt[:2], color="white", radius=elt[2]).add_to(m)

m

输出：

更新（2）：

这是使用以下方法在网格上进行插值的版本griddata：

import pandas as pd
import numpy as np
import folium
from scipy.interpolate import griddata

# parameters
n = 250                     # number of points
lat0 = 40.7
lon0 = -73.9
eps = 0.1
v_min, v_max = 0, 100       # min, max values

# generating values
lat = np.random.normal(lat0, eps, n)
lon = np.random.normal(lon0, eps, n)
value = numpy.random.uniform(v_min, v_max, n)

# set up the grid
step = 0.02
xi, yi = np.meshgrid(
    np.arange(lat.min() - step/2, lat.max() + step/2, step),
    np.arange(lon.min() - step/2, lon.max() + step/2, step),
)

# interpolate and normalize values
zi = griddata((lat, lon), value, (xi, yi), method='linear')
zi /= np.nanmax(zi)
g = np.stack([
    xi.flatten(),
    yi.flatten(),
    zi.flatten(),
], axis=1)

# geo_json returns a single square
def geo_json(lat, lon, value, step):
    cmap = mpl.cm.RdBu
    return {
      "type": "FeatureCollection",
      "features": [
        {
          "type": "Feature",
          "properties": {
            'color': 'white',
            'weight': 1,
            'fillColor': mpl.colors.to_hex(cmap(value)),
            'fillOpacity': 0.5,
          },
          "geometry": {
            "type": "Polygon",
            "coordinates": [[
                [lon - step/2, lat - step/2],
                [lon - step/2, lat + step/2],
                [lon + step/2, lat + step/2],
                [lon + step/2, lat - step/2],
                [lon - step/2, lat - step/2],
              ]]}}]}


# generating a map...
m = folium.Map(location=[lat0, lon0], zoom_start=9, width=400, height=400)

# ...with squares...
for gi in g:
    if ~np.isnan(gi[2]):
        folium.GeoJson(geo_json(gi[0], gi[1], gi[2], step),
                       lambda x: x['properties']).add_to(m)

# ...and the original points
for elt in list(zip(lat, lon, value)):
    folium.Circle(elt[:2], color='white', radius=elt[2]).add_to(m)

m

输出：