绘制地理空间矢量场

用 Folium 绘制地理空间矢量场

地学和许多应用领域中，数据的视觉化非常重要。尤其是一些表示方向和速度的矢量数据，例如风速、海流、车速等，使用矢量图进行绘制能够更加直观地表达这些数据的特性。

示例数据集选择

为了便于说明矢量场的绘制方法，我们选择了一个实际的例子：纽约市的 Citi Bike 数据集。Citi Bike 是纽约市一个著名的自行车共享系统，该数据集涵盖了数百万次骑行，包含每次骑行的起点和终点位置。我们将利用这些信息计算每次骑行的速度向量。

数据准备与处理

首先，我们需要导入所需的Python库并准备工作环境。以下是我们将使用的主要库：

「Folium」：一个用于在Python中绘制交互式地图的库。
「H3」：一个地理空间索引库，支持六边形网格划分。
「Pandas」 和 「Numpy」：用于数据处理和数值计算。
「Matplotlib」：用于数据可视化。
「Pyproj」：用于地理计算。
「SciPy」：用于统计计算。

我们从本地的存储目录中加载Citi Bike数据集，读取每个CSV文件并将它们合并为一个DataFrame。

import jsonimport matplotlib.pyplot as plt
import numpy as npimport h3import folium
from geojson import Feature, FeatureCollectionimport math
import matplotlibfrom scipy import stats
import pandas as pdfrom pandarallel import pandarallelimport scipy.stats as scsimport os
import pandas as pd
import zipfileos.chdir('/home/mw/input/tripdata5123')zf = os.listdir('/home/mw/input/tripdata5123/202406-citibike-tripdata')
result_data = []
print("Reading:")for zfcsv in zf:print(zfcsv)df = pd.read_csv('/home/mw/input/tripdata5123/202406-citibike-tripdata/' + zfcsv, low_memory=False)result_data.append(df)df = pd.concat(result_data)
print("Done!")
result_data = None

数据清理

在进行分析之前，通常需要对数据进行清理。我们首先去掉所有缺失值（NaN）以确保后续计算的准确性。

df = df.dropna()
print(len(df))
#df = df.sample(10000)

计算骑行中点和近似航向

为了将每次骑行归纳到一个中心点，我们计算了起点和终点的中点。同时，为了简化计算，我们使用近似值来估算航向（从起点到终点的方向角度）。需要注意的是，这种近似方式并不完全精确。

# we choose the middle position instead of the start or end position to plot
df['latitude'] = (df['start_lat'] + df['end_lat'])/2.0
df['longitude'] = (df['start_lng'] + df['end_lng'])/2.0# we approximate bearing like this -- no, it is not entirely correct!
df['bearing_approx'] = 90.0 - 360 * (np.arctan2(df['end_lat'] - df['start_lat'], df['end_lng'] - df['start_lng']) / (2*math.pi))df['timestamp'] = pd.to_datetime(df['ended_at'])
df['duration'] = df['timestamp'] - pd.to_datetime(df['started_at'])
df['hour_of_day'] = df['timestamp'].dt.hour

计算距离和速度

接下来，我们使用H3库来计算每次骑行的直线距离，并利用计算出的距离和骑行时间来计算速度（km/h）。

pandarallel.initialize(progress_bar=True)distances = df.parallel_apply(lambda row: h3.point_dist((row.start_lat, row.start_lng), (row.end_lat, row.end_lng), unit='m'), axis = 1)df = df.assign(distance_in_meters=distances.values)
distances = None

speeds = df.parallel_apply(lambda row: 3.6 * row.distance_in_meters / row.duration.total_seconds(), axis = 1)df = df.assign(speed_kmh=speeds.values)
speeds = None

计算精确航向

为了得到更加精确的航向，我们使用了 pyproj 库中的 Geod 对象。这可以提供更精确的方位角计算。

import pyprojgeodesic = pyproj.Geod(ellps='WGS84')def bearing_from_lat_long_pair(start_lat, start_long, end_lat, end_long):bearing_angle, _, _ = geodesic.inv(start_long, start_lat, end_long, end_lat)if bearing_angle < 0:bearing_angle = bearing_angle + 360.0if bearing_angle > 360.0:bearing_angle = bearing_angle - 360.0return bearing_anglebearings_properly = df.parallel_apply(lambda row: bearing_from_lat_long_pair(row.start_lat, row.start_lng, row.end_lat, row.end_lng), axis = 1)df = df.assign(bearing=bearings_properly.values)bearings_properly = None

计算方向向量的x和y分量

为了更好地处理方向平均值，我们将每个方向的角度转换为x和y分量。这样，我们可以使用标准的聚合函数来计算方向的平均值。

df["direction_x"] = np.sin(2 * math.pi * df["bearing"] / 360.0)
df["direction_y"] = np.cos(2 * math.pi * df["bearing"] / 360.0)

将数据聚合到六边形网格中

现在，我们将数据聚合到六边形网格中。在本文的例子中，我们使用H3库将数据转换为分辨率为9的六边形网格单元。

df.groupby("rideable_type", as_index=False).agg({"ride_id": "count", "speed_kmh": "mean" })
# CONSTANTSH3_RESOLUTION = 9pandarallel.initialize(progress_bar=True)hex_ids = df.parallel_apply(lambda row: h3.geo_to_h3(row.latitude, row.longitude, H3_RESOLUTION), axis = 1)df = df.assign(hex_id=hex_ids.values)
hex_ids = None

现在我们有了数百万个速度矢量，我们需要将它们聚合到网格单元中，本文的例子中是六边形。这里有个细节值得一提。取「速度」的平均值（或中位数，或其他值）很简单，你可以简单地用*“speed”: “mean”执行**groupby('hexid', …) **。取方位的平均值不能以同样的方式进行，因为方位具有360 度的圆形范围（弧度，则是 2π）。

例如，2 度（几乎向北；略微向东）和 356 度（几乎向北；略微向西）的平均值肯定不是「179」度（几乎向南！）。相反，在方向的世界中，合理的平均值是 359 度，它是一个 2 度角和一个 356 度角的圆平均值。

import numpy as np
import scipy.stats# compass bearings in degrees:
data = [2, 356]# average bearing:
print(np.mean(data)) # == 179.0, aka "wrong"
print(scipy.stats.circmean(data, high=360, low=0)) # == 359.0, aka "correct"print(np.rad2deg(np.angle(np.mean(np.cos(np.deg2rad(data)) + np.sin(np.deg2rad(data)) * 1j))%(2*math.pi)))
# == 359.0, aka also correct, but achieving the result in an overly complex way

当然，有现成的函数可以计算圆周均值 -例如scipy.stats.circmean - 但不知道如何将它与groupby一起使用，所以我采用了同样简单但不太优雅的方法，即独立平均速度矢量的 x 和 y 分量，并将arctan2函数应用于结果。

计算六边形网格单元的平均方向和速度

计算平均速度相对简单，可以直接使用 groupby 聚合函数。然而，对于方向的平均值，我们需要使用向量平均的方法。通过将方向角转换为x和y分量，然后计算这些分量的平均值，我们可以得出每个六边形单元的平均方向。

positions_df_byhex = df_subset.groupby("hex_id", as_index=False).agg({"longitude": "mean","latitude": "mean","direction_x": "mean","direction_y": "mean","bearing": "count","speed_kmh": "mean"})

另一个容易犯的错误是混淆方位（通常以北为零，顺时针测量）和数学角度（以正 x 轴为零，逆时针测量）。还有一点：如果您的矢量指的是风，风向通常是指风吹来的方向，而不是空气流动的方向。

❝
将这些错误加在一起，很容易就会得到本来应该是逆时针的方向，但却是顺时针方向，「和/或」偏移 +90、-90 或 +180 度，「而且」从一开始就是错误聚合的。甚至没有提到弧度与度数。因此，要小心！
❞

将自行车行程速度向量聚合到六边形单元格中，就能绘制它们了。使用六边形热图是不够的，因为对于每个单元格，除了「幅度（速度）之外，我们还有一个向量方向」（方位）。因此，我们在代码中手动设计了一个小箭头符号，并用这些箭头替换六边形，颜色可以反映速度。

arrow_tip_x = size_of_arrow_x * math.cos(bearing_rad) + center_lat
arrow_tip_y = size_of_arrow_y * math.sin(bearing_rad) + center_longarrow_back_left_x = size_of_arrow_x * math.cos(bearing_rad + 2*math.pi * 135.0 / 360.0 ) + center_lat
arrow_back_left_y = size_of_arrow_y * math.sin(bearing_rad + 2*math.pi * 135.0 / 360.0 ) + center_longarrow_back_x = size_of_arrow_x * 0.2 * math.cos(bearing_rad + 2*math.pi * 180.0 / 360.0 ) + center_lat
arrow_back_y = size_of_arrow_y * 0.2 * math.sin(bearing_rad + 2*math.pi * 180.0 / 360.0 ) + center_longarrow_back_right_x = size_of_arrow_x * math.cos(bearing_rad + 2*math.pi * 225.0 / 360.0 ) + center_lat
arrow_back_right_y = size_of_arrow_y * math.sin(bearing_rad + 2*math.pi * 225.0 / 360.0 ) + center_longcoordinates = ( (arrow_tip_x, arrow_tip_y), (arrow_back_left_x, arrow_back_left_y), (arrow_back_x, arrow_back_y), (arrow_back_right_x, arrow_back_right_y), (arrow_tip_x, arrow_tip_y) )geometry_for_row = { "type" : "Polygon", "coordinates": [ coordinates ] }

我们选择速度百分位数来利用全色彩范围，但对于某些应用程序，您可能希望使用原始值。