在大型 GIS 数据库中按属性高效溶解相邻多边形-海口c网

您是否拥有一个大型 GIS 数据集，并希望高效地融合所有具有相同属性的相邻多边形？在本文中，我将分享如何使用 PostGIS 处理包含超过 75 万行数据的土地利用数据集来实现这一目标。

我将以维多利亚州土地利用数据集为例。该数据可从Data VIC免费获取。

动机

geopandas首先，让我们使用Python绘制一小组数据样本，这样我就可以解释我想要溶解多边形的动机。

由于数据保存在 PostGIS 表中，我们首先要创建与数据库的连接。

import psycopg
from sqlalchemy import create_engine# database configuration
user = '<username>'
host = 'localhost'
database = '<dbname>'
driver = 'postgresql+psycopg'
connection_str = f'{driver}://{user}@{host}/{database}'engine = create_engine(connection_str)

我只想对数据子集进行采样，因此我只会过滤随机位置 (-36.8, 144.1) 0.025 度范围内的多边形。

import geopandas as gpdsql = """
SELECT alumv8, geometry
FROM vic_landuse
WHERE ST_DWithin(geometry,ST_GeomFromText('POINT(144.1 -36.8)', 7844), 0.025)AND ST_geometrytype(geometry) = 'ST_Polygon'
"""sample = gpd.read_postgis(sql, engine, geom_col="geometry")
sample['alumv8'] = sample['alumv8'].astype('category')ax = sample.plot('alumv8',legend=True,legend_kwds=dict(loc='upper left',bbox_to_anchor=(1.05, 1),title='ALUM codes')
)
ax.tick_params('x', rotation=45)

没有任何

在上面的地图中，我们有许多相邻的地块具有相同的属性（在本例中为alumv8）。我想根据土地用途创建自定义图块；因此我对单个地块不感兴趣。出于处理目的，我可以通过合并具有相同属性的相邻多边形来减少几何图形的数量。

使用数据子集进行演示

首先，我将逐步演示如何使用上面创建的小子集数据框来实现这一点。然后，我们可以创建一个脚本来自动化整个表的流程。

%load_ext sql
%sql postgresql+psycopg://<username>@localhost/<dbname

我将使用示例数据创建一个新表。

%%sql
CREATE TABLEsample_landuse AS
SELECT*
FROMvic_landuse
WHEREST_DWithin (geometry,ST_GeomFromText ('POINT(144.1 -36.8)', 7844),0.025)AND ST_geometrytype (geometry) = 'ST_Polygon'
%%sql
SELECTalumv8,COUNT(alumv8)
FROMsample_landuse
GROUP BYalumv8;

没有任何

每个 ALUM 代码中的多边形数量。

简要流程

我遵循的流程是：

为每个 ALUM 代码创建一个视图。
在每个代码中，如果几何图形相邻，则将其溶解，从而生成一个单独的表。
连接在步骤 2 中创建的所有表。
清理所有中间视图和表。

我将通过下面的示例演示这些步骤。

为每个 ALUM 代码创建一个视图

这很简单。在本例中，我们使用 ALUM 代码320。

%%sql
CREATE VIEWlanduse_320 AS
SELECT*
FROMsample_landuse
WHEREalumv8 = 320;

溶解每个代码中的相邻几何图形

ST_ClusterDBSCAN是一个窗口函数，它返回每个输入几何体的聚类编号。相邻的几何体会被赋予相同的聚类编号。

我们将使用ST_ClusterDBSCAN(geometry, 0, 2)，因为我们希望将彼此相距 0 个单位（度）以内的几何图形进行聚类，并且我们希望聚类中至少有两个几何图形。

%%sql
CREATE TABLElanduse_320_clustered AS
SELECTalumv8,ST_ClusterDBSCAN (geometry, 0, 2) OVER () AS cluster_id,geometry
FROMlanduse_320;

我们可以看到已经创建了多个至少包含 2 个多边形的集群。

%%sql
SELECTcluster_id,COUNT(cluster_id)
FROMlanduse_320_clustered
GROUP BYcluster_id
ORDER BYCOUNT desc
LIMIT5;

没有任何

每个簇内的多边形数量。

接下来，如果行具有相同的cluster_id ，我们将使用 PostGIS 的 ST_Union 方法分解这些行。

%%sql
CREATE TABLElanduse_320_dissolved AS
SELECTalumv8,ST_Union (geometry) AS geometry
FROMlanduse_320_clustered
GROUP BYalumv8,cluster_id;

现在，我们只有 10 个多边形，而不是最初开始的 30 个。

%sql SELECT COUNT(*) FROM landuse_320_dissolved;

没有任何

重复所有 ALUM 代码

现在，我们可以循环遍历表中的所有 ALUM 代码，并对每个代码重复上述过程。

BEGIN
-- the new table
CREATE TABLE landuse_dissolved (alumv8 int, geometry geometry);FOR alum_code INSELECT DISTINCT alumv8 FROM sample_landuse
LOOPEXECUTE format('DROP VIEW IF EXISTS landuse_%s', alum_code);-- create view for the codeEXECUTE format('CREATE VIEW landuse_%s ASSELECT * FROM sample_landuseWHERE alumv8 = %L',alum_code, alum_code);-- create clustered tableEXECUTE format('CREATE TABLE landuse_%s_clustered ASSELECTalumv8,ST_ClusterDBSCAN(geometry, 0, 2) OVER() AS cluster_id,geometryFROM landuse_%s', alum_code, alum_code);-- dissolve adjacent geometriesEXECUTE format('CREATE TABLE landuse_%s_dissolved ASSELECT alumv8, ST_Union(geometry) as geometryFROM landuse_%s_clusteredGROUP BY alumv8, cluster_id',alum_code, alum_code);-- add dissolved geometries to new tableEXECUTE format('INSERT INTO landuse_dissolved (alumv8, geometry)SELECT alumv8, geometryFROM landuse_%s_dissolved', alum_code);-- cleanup intermediary views and tablesEXECUTE format('DROP TABLE landuse_%s_clustered', alum_code);EXECUTE format('DROP VIEW landuse_%s', alum_code);EXECUTE format('DROP TABLE landuse_%s_dissolved', alum_code);COMMIT;
END LOOP;
END $$
original_rows = %sql SELECT COUNT(*) FROM sample_landuse;
new_rows = %sql SELECT COUNT(*) FROM landuse_dissolved;print(f"original rows = {original_rows[0].count}")
print(f"rows after dissolving = {new_rows[0].count}")

原始行数 = 94 溶解后行数 = 29

因此，在溶解相邻多边形之后，多边形数量从原来的 94 个减少到 29 个。现在让我们将溶解后的表格绘制landuse_dissolved在原始表格旁边sample，并直观地检查上述过程的影响。

import matplotlib.pyplot as pltsql = "SELECT * FROM landuse_dissolved"
sample_dis = gpd.read_postgis(sql, engine, geom_col="geometry")
sample_dis['alumv8'] = sample_dis['alumv8'].astype('category')fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True)sample.plot('alumv8', ax=ax1)
sample_dis.plot('alumv8', ax=ax2)for ax in [ax1, ax2]:ax.set_xticks([])ax.set_yticks([])

没有任何

左图：原始数据集。右图：使用 ALUM 代码溶解相邻多边形后的结果。

我们可以看到，除了具有相同 ALUM 代码的相邻几何体中单个地块的边界外，没有任何信息丢失。因此，我们可以确认获得了所需的结果，现在可以对整个数据集进行相同的操作。

面积小了？

这里需要注意的是，总面积会小于原始面积。这主要是由于融合过程中几何形状被简化，或者消除了多边形之间的重叠区域。

让我们计算每个 ALUM 代码的总面积，并比较两个表的结果。注意：由于几何体的 CRS 是地理坐标系（EPSG:7844，单位为度），因此我们需要将几何体转换为相应的平面 CRS（EPSG:3111，单位为米）来计算面积。

我们可以创建一个数据框，记录面积差异（以平方米为单位）。结果显示，123个ALUM代码中，有77个的总面积存在差异。

%%sql
original_area <<
SELECTalumv8,SUM(ST_Area (ST_Transform (geometry, 3111))) AS area
FROMsample_landuse
GROUP BYalumv8;
%%sql
new_area <<
SELECTalumv8,SUM(ST_Area (ST_Transform (geometry, 3111))) AS area
FROMlanduse_dissolved
GROUP BYalumv8;

我们可以创建一个数据框，记录面积差异（以平方米为单位）。结果显示，123个ALUM代码中，有77个的总面积存在差异。

area= original_area.DataFrame().set_index('alumv8').join(new_area.DataFrame().set_index('alumv8'), how='inner', rsuffix='_new')
area['diff'] = area['area'] - area['area_new']
area['diff_pct'] = 100 * (area['area'] - area['area_new']) / area['area']area

没有任何

原始数据集和溶解数据集之间每个 ALUM 代码的总面积（和差异）。

我们可以看到，5/8 ALUM 代码的面积有所不同。虽然对于此示例数据集而言，差异很小，但当数据集较大时，差异可能会累积起来。

因此，在决定这种溶解是否适合当前任务时，牢记这一点很重要。

溶解整个数据集中的相邻多边形

当然，对于超过 75 万行的整个数据集，使用相同的 ALUM8 代码溶解所有相邻多边形是一个昂贵且耗时的过程。或许可以休息一会儿喝杯咖啡，或者舒展一下筋骨？

但结果相当令人印象深刻。最终只有约 112k 行，减少了 85%。

original_rows = %sql SELECT COUNT(*) FROM vic_landuse;
new_rows = %sql SELECT COUNT(*) FROM vic_landuse_dissolved;print(f"original rows = {original_rows[0].count}")
print(f"rows after dissolving = {new_rows[0].count

原始行数 = 751295 溶解后行数 = 112603

在某些任务中，例如使用减少的行数进行绘图，我们可以立即看到差异。在我的电脑上，使用geopandas绘制原始数据集大约需要 2 分钟，而处理合并后的数据集只需要大约 50 秒。

sql = "SELECT * FROM vic_landuse"
vic_landuse = gpd.read_postgis(sql, engine, geom_col="geometry")
vic_landuse["alumv8"] = vic_landuse["alumv8"].astype("category")%time vic_landuse.plot('alumv8')

没有任何