如何加速千万级地理坐标数据集的地形剖面构建速度？

现有1200万行数据的CSV地形数据集，需通过坐标构建地形剖面图。原CSV单高度查询耗时约2-2.5秒，转Parquet后缩短至1-1.7秒，但构建500-2000个点的剖面仍耗时过长，且未来数据集可能扩容，处理速度会更慢。求进一步降低处理耗时的方案。

代码示例：

import dask.dataframe as dk
import numpy as np
import pandas as pd
import time

filename = 'n46_e032_1arc_v3.csv'

df = dk.read_csv(filename)
df.to_parquet('n46_e032_1arc_v3_parquet')

Latitude1y, Longitude1x = 46.6276, 32.5942
Latitude2y, Longitude2x = 46.6451, 32.6781
sec, steps, k = 0.00027778, 1, 11.73
Latitude, Longitude = [Latitude1y], [Longitude1x]
sin, cos = Latitude2y - Latitude1y, Longitude2x - Longitude1x
y, x = Latitude1y, Longitude1x
while Latitude[-1] < Latitude2y and Longitude[-1] < Longitude2x:
    y, x, steps = y + sec * k * sin, x + sec * k * cos, steps + 1
    Latitude.append(y)
    Longitude.append(x)

time_start = time.time()
long, elevation_data = [], []
df2 = dk.read_parquet('n46_e032_1arc_v3_parquet')
for i in range(steps + 1):
    elevation_line = df2[(Longitude[i] <= df2['x']) & (df2['x'] <= Longitude[i] + sec) &
                    (Latitude[i] <= df2['y']) & (df2['y'] <= Latitude[i] + sec)].compute()
    elevation = np.asarray(elevation_line.z.tolist())
    if elevation[-1] < 0:
        elevation_data.append(0)
    else:
        elevation_data.append(elevation[-1])
    long.append(30 * i)

plt.bar(long, elevation_data, width = 30)
plt.show()
print(time.time() - time_start)

1. 替换循环单查为批量处理 #

当前代码的核心瓶颈是循环中对每个点单独调用compute()，每次都会触发Dask的任务调度和数据分片读取，重复开销极大。可以先一次性筛选出所有目标点覆盖的区域数据，再在内存中完成匹配：

# 生成所有查询区间的元数据
query_df = pd.DataFrame({
    'min_x': Longitude,
    'max_x': [lon + sec for lon in Longitude],
    'min_y': Latitude,
    'max_y': [lat + sec for lat in Latitude],
    'point_idx': range(len(Longitude))
})

# 先缩小全局查询范围，减少加载的数据量
x_min, x_max = query_df['min_x'].min(), query_df['max_x'].max()
y_min, y_max = query_df['min_y'].min(), query_df['max_y'].max()
full_data = dk.read_parquet('n46_e032_1arc_v3_parquet')
filtered_data = full_data[(full_data['x'].between(x_min, x_max)) & 
                          (full_data['y'].between(y_min, y_max))].compute()

# 内存中批量匹配每个点的高程
def match_elevation(row):
    mask = (filtered_data['x'].between(row['min_x'], row['max_x'])) & \
           (filtered_data['y'].between(row['min_y'], row['max_y']))
    elevations = filtered_data.loc[mask, 'z'].values
    return 0 if len(elevations) == 0 or elevations[-1] < 0 else elevations[-1]

query_df['elevation'] = query_df.apply(match_elevation, axis=1)
elevation_data = query_df['elevation'].tolist()
long = query_df['point_idx'].apply(lambda x: 30 * x).tolist()

2. 给Parquet文件添加空间分区 #

针对经度x和纬度y字段对Parquet文件做分区，查询时Dask会自动跳过无关分区，大幅减少磁盘IO量：

# 写入Parquet时按经纬度分箱分区（可根据数据分布调整分箱粒度，比如0.1度）
raw_df = dk.read_csv(filename)
raw_df['x_bin'] = (raw_df['x'] // 0.1).astype(int)
raw_df['y_bin'] = (raw_df['y'] // 0.1).astype(int)
raw_df.to_parquet('n46_e032_1arc_v3_parquet_partitioned', partition_on=['x_bin', 'y_bin'])

# 后续查询时，Dask只会加载目标点所在分区的数据，无需遍历全量
df_part = dk.read_parquet('n46_e032_1arc_v3_parquet_partitioned')

3. 构建KD树实现快速空间搜索 #

如果服务器内存足够容纳全量数据，构建KD树可以将空间查询复杂度从O(n)降到O(logn)，适合高频坐标查询场景：

import scipy.spatial as spatial

# 加载全量数据到内存
full_data = dk.read_parquet('n46_e032_1arc_v3_parquet').compute()
# 用经纬度构建KD树
coords = full_data[['x', 'y']].values
kd_tree = spatial.KDTree(coords)

# 批量查询每个点的高程
elevation_data = []
for lon, lat in zip(Longitude, Latitude):
    # 搜索sec范围内的所有点
    indices = kd_tree.query_ball_point([lon, lat], r=sec)
    if not indices:
        elevation_data.append(0)
    else:
        elev = full_data.iloc[indices]['z'].values[-1]
        elevation_data.append(0 if elev < 0 else elev)

4. 转换为栅格格式（终极优化） #

你的数据是1arc分辨率的规则格网地形数据，转换为TIFF栅格后，读取高程的速度会比表格形式快一个数量级：

import rasterio
from rasterio.transform import from_origin

# 加载数据并提取栅格参数
full_data = dk.read_parquet('n46_e032_1arc_v3_parquet').compute()
x_min, x_max = full_data['x'].min(), full_data['x'].max()
y_min, y_max = full_data['y'].min(), full_data['y'].max()
res = sec  # 1arc对应的分辨率

# 创建栅格变换矩阵
transform = from_origin(x_min, y_max, res, res)
width = int((x_max - x_min) / res)
height = int((y_max - y_min) / res)
raster = np.zeros((height, width), dtype=np.float32)

# 填充栅格数据
for _, row in full_data.iterrows():
    x_idx = int((row['x'] - x_min) / res)
    y_idx = int((y_max - row['y']) / res)
    raster[y_idx, x_idx] = row['z']

# 保存为TIFF文件
with rasterio.open('terrain.tif', 'w', driver='GTiff',
                   height=height, width=width, count=1,
                   dtype=raster.dtype, crs='EPSG:4326',
                   transform=transform) as dst:
    dst.write(raster, 1)

# 后续读取高程直接通过坐标索引
with rasterio.open('terrain.tif') as src:
    elevation_data = []
    for lon, lat in zip(Longitude, Latitude):
        row, col = src.index(lon, lat)
        if 0 <= row < height and 0 <= col < width:
            elev = src.read(1)[row, col]
            elevation_data.append(0 if elev < 0 else elev)
        else:
            elevation_data.append(0)

内容的提问来源于stack exchange，提问作者Alejandro