地学数据处理提速
地理空间处理因其资源密集型需求而非常慢，尤其是在处理1G+的栅格数据时。基于向量的处理也会使系统资源紧张，在大型数据集的情况下。
Python 生态系统有几个关键地理空间库构建，包括 Geopandas、Shapely、Fiona 和 Pyproj。虽然 GeoPandas 可以高效地处理地理空间数据，但数据量大的时候，仍可能成为数据处理的瓶颈，尤其是在处理需要多次迭代的大型数据集时。
最近我在用SWOT卫星数据观察河道和湖泊时，一个湖泊的数据近 900 万个数据点，需要数万个空间查询并与其他图层进行关联。导致渲染特别慢。
这是一个很好的例子，说明这时候我们要使用GeoPandas 加速优化技术了
技巧1：使用pyogrio新引擎
在 GeoPandas 中，所选的底层引擎会极大地影响数据读取效率。从 0.11 版本（2022 年 7 月发布）开始，GeoPandas 引入了一个名为 pyogrio 的新引擎，用于读取矢量数据集。但是默认引擎仍然是 Fiona。
为了提高性能，我们可以通过在 read_file 命令中传递 engine=' pyogrio' 来切换到 pyogrio 引擎。如果我们选择 pyogrio，我们可以通过设置 use_arrow=True 来进一步增强从 GDAL 到 Python 的数据传输。这利用 Arrow 进行更高效的数据处理。
代码和测试如下：
import time
import matplotlib.pyplot as plt
import geopandas as gpd
file_path = '/home/mw/input/river3479/geoft_bho_2017_5k_curso_dagua.gpkg'
engines = {
'fiona': {'engine': 'fiona'},
'pyogrio': {'engine': 'pyogrio'},
'pyogrio+arrow': {'engine': 'pyogrio', 'use_arrow': True}
}
times = []
for engine in engines:
start_time = time.time()
df = gpd.read_file(file_path, **engines[engine])
end_time = time.time()
times.append(end_time - start_time)
# Plotting
plt.figure(figsize=(10, 6))
plt.bar(engines.keys(), times, color='skyblue')
plt.xlabel('Engine')
plt.ylabel('Time (seconds)')
plt.title('Performance of GeoPandas `read_file` function with Different Engines')
plt.show()
Fiona用了40秒，而pyogrio+arrow10秒内解决战斗，这下不是标题党了吧。。
技巧2：表格筛选
使用大型数据集时，读取整个数据集可能非常耗时。可以在读取数据之前根据特定条件或条件筛选数据子集
SQL 查询过滤：从 0.12 版本开始，GeoPandas 还支持基于 SQL 查询的过滤。使用此方法，仅加载满足查询条件的数据。这是在处理大型数据集时优化性能的有效方法。
times = []
# without tabular filtering
start_time = time.time()
df = gpd.read_file(file_path, **engines['pyogrio+arrow'])
end_time = time.time()
times.append(end_time - start_time)
# with tabular filtering
start_time = time.time()
df = gpd.read_file(file_path, **engines['pyogrio+arrow'], where="cocursodag='676'")
end_time = time.time()
times.append(end_time - start_time)
# Plotting
fig, ax = plt.subplots(figsize=(6, 3))
ax.bar(['without filtering', 'with filtering'], times, color='skyblue')
ax.set_ylabel('Time (seconds)')
ax.set_title('Different masking times')
没有经过过滤需要7秒，而SQL过滤后仅仅1秒。
技巧3：边界框筛选
利用 pyogrio 引擎，考虑边界框筛选。 它的工作原理如下：不是精确定义多边形，而是加载一个包含所需区域的更大区域（边界框）。 虽然这种方法读取的数据比绝对必要的要多，但由此产生的性能改进是非常值得的。
这可能有点反直觉，但是速度要更快。
# load the brazilian states
states = gpd.read_file("/home/mw/input/river3479/br_states.json")
se = states.query("PK_sigla == 'SE'")
# setup the scenarios
engines = {
'fiona+mask': {'engine': 'fiona', 'mask': se},
'pyogrio+arrow+bbox': {'engine': 'pyogrio', 'use_arrow': True, 'bbox': tuple(se.total_bounds)}
}
# loop through the scenarios
times = []
for engine in engines:
start_time = time.time()
df = gpd.read_file(file_path, **engines[engine])
end_time = time.time()
times.append(end_time - start_time)
# Plotting
fig, axs = plt.subplots(1, 2, figsize=(10, 5))
axs[0].bar(engines.keys(), times, color='skyblue')
axs[0].set_xlabel('Engine')
axs[0].set_ylabel('Time (seconds)')
axs[0].set_title('Different masking times')
se.plot(ax=axs[1], facecolor='none')
df.plot(ax=axs[1], zorder=-1, linewidth=0.5)
axs[1].set_title('Final dataset')
技巧4：坐标索引
高效的空间操作同样重要。
在第一种情况下，我们在河流数据集执行简单的空间连接
在第二种场景中，我们在使用坐标索引 CX 运算符执行空间操作之前，按状态的边界框过滤数据集。
times = []
reaches = states
# without CX indexing
start_time = time.time()
join = reaches.sjoin(se)
end_time = time.time()
times.append(end_time - start_time)
# Using CX indexing prior the spatial join
start_time = time.time()
xmin, ymin, xmax, ymax = se.total_bounds
reaches2 = reaches.cx[xmin:xmax, ymin:ymax]
join2 = reaches2.sjoin(se)
end_time = time.time()
times.append(end_time - start_time)
# Plotting
plt.figure(figsize=(10, 6))
plt.bar(['Without CX filtering', ' WIth CX filtering'], times, color='skyblue')
plt.xlabel('Engine')
plt.ylabel('Time (seconds)')
plt.title('Performance comparison with and without CX indexing')
技巧5：并行处理
并行处理（使用 concurrent.futures）：
我们创建一个 PoolExecutor 来管理并行进程。 通过将工作负载分配到多个内核，我们可以获得更好的性能。
虽然提升不显著（约 35%）在这种情况下，改进程度取决于硬件规格、计算核心数和具体的工作负载。并行处理可以显著提高效率，尤其是对于大规模地理空间任务！
技巧6：几何简化
当其他优化技术无法满足要求时，请考虑创建更简单的地理空间数据集表示形式。GeoPandas 中的 simplify（） 方法是实现此目的的强大工具。
from pathlib import Path
simplified = reaches.simplify(tolerance=0.5)
# measure clipping times for AL state
times = []
start_time = time.time()
reaches.clip(states.query("PK_sigla == 'AL'"))
end_time = time.time()
times.append(end_time - start_time)
start_time = time.time()
simplified.clip(states.query("PK_sigla == 'AL'"))
end_time = time.time()
times.append(end_time - start_time)
# save the bases to disk
reaches.to_file("/home/mw/temp/original_bho.gpkg", driver='GPKG')
simplified.to_file("/home/mw/temp/simplified_bho.gpkg", driver='GPKG')
size_original = Path("/home/mw/temp/original_bho.gpkg").stat().st_size
size_simplified = Path("/home/mw/temp/simplified_bho.gpkg").stat().st_size
# Plotting file sizes
_, axs = plt.subplots(1, 2 ,figsize=(10, 5))
axs[0].bar(['Original', 'Simplified'], [size_original / 1024**2, size_simplified / 1024**2], color='skyblue')
axs[0].set_ylabel('Size (MB)')
axs[0].set_title('File Sizes')
# Plotting clipping times
axs[1].bar(['Original', 'Simplified'], times, color='skyblue')
axs[1].set_ylabel('Time (s)')
axs[1].set_title('Clipping Times')
除了裁剪性能提高 30% 之外，简化版本占用的磁盘空间也明显减少
求求你点个在看吧，这对我真的很重要
