03 PostGIS如何有效比较2D/3D两种维度下的几何是否相同

一 背景

近期遇到一些需求，使用PostGIS存储几何数据，做的过程中需求产生过一些变化，大概如下：

1. 服务端采用PostGIS存储几何数据，几何类型为Point、LineString、Polygon、MultiPoint、MultiLineString、MultiPolygon等常见的几种几何类型，维度为二维，不包含Z坐标；
2. 创建几何时采用几何信息的WKT格式进行创建，例如“POINT (x y)”；
3. 创建时如果几何存在则返回指定记录id（涉及到几何对比）；
4. 新需求，几何允许Z坐标的存储；

在要求存储Z坐标之前，解决的主要思路为：

1. 服务端接受几何WKT数据构建几何对象，实体中几何属性采用org.postgis.geometry或者org.locationtech.jts.geom；
2. 自定义配置MyBatisPlus对postgis的几何（postgis-Jdbc）与自定义类型几何（org.postgis.geometry或者org.locationtech.jts.geom)的转换；
3. 采用MyBatis/MyBatisPlus进行几何数据记录数据的CRUD，插入几何记录前调用在Mapper.xml层自定义的是否有指定几何存在的方法（PostGIS的ST_Equals(geo1,geo2)函数）来判断是否有存在的几何记录；

整个思路在未添加Z维度坐标之前是正确的，但是在添加Z坐标维度后，有些地方需要进行调整，主要为：

1. 实体中几何属性不建议采用GeoTools（org.locationtech.jts.geom）中组织的几何，因为其不具有Z坐标的支持，可使用postgis的几何（postgis-Jdbc）替换（此部分不是重点，仅提醒一下）；

以点类型对比，虽然GeoTools中的Point类中Z信息也存在于Coordinate中，但Point类并没有很好的描述，所以如果存在Z甚至M等多维度数据，建议使用postgis-jdbc中对几何的组织方式。

      	<!-- postgis数据库驱动 -->
        <dependency>
            <groupId>net.postgis</groupId>
            <artifactId>postgis-jdbc</artifactId>
        </dependency>

2. 数据库中建立几何表，采用3D几何类型创建
   1. Point -> PointZ
   2. Polygon -> PolygonZ
   3. LineString -> LineStringZ
   4. MultiPoint ->MultiPointZ
   5. MultiLineString -> MultiLineStringZ
   6. MultiPolygon -> MultiPolygonZ

-- ----------------------------
-- 添加PostGIS扩展
-- ----------------------------
-- CREATE EXTENSION postgis;


DROP TABLE IF EXISTS "public"."data_point";
CREATE TABLE data_point (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(Point, 4326)
);
CREATE INDEX data_point_geom_index ON data_point USING GIST (geometry);

DROP TABLE IF EXISTS "public"."data_pointz";
CREATE TABLE data_pointz (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(PointZ, 4326)
);
CREATE INDEX data_pointz_geom_index ON data_pointz USING GIST (geometry);



DROP TABLE IF EXISTS "public"."data_line";
CREATE TABLE data_line (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(LineString, 4326)
);
CREATE INDEX data_line_geom_index ON data_line USING GIST (geometry);

DROP TABLE IF EXISTS "public"."data_linez";
CREATE TABLE data_linez (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(LineStringZ, 4326)
);
CREATE INDEX data_linez_geom_index ON data_linez USING GIST (geometry);



DROP TABLE IF EXISTS "public"."data_polygon";
CREATE TABLE data_polygon (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(Polygon, 4326)
);
CREATE INDEX data_polygon_geom_index ON data_polygon USING GIST (geometry);

DROP TABLE IF EXISTS "public"."data_polygonz";
CREATE TABLE data_polygonz (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(PolygonZ, 4326)
);
CREATE INDEX data_polygonz_geom_index ON data_polygonz USING GIST (geometry);



DROP TABLE IF EXISTS "public"."data_point_multi";
CREATE TABLE data_point_multi (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(MultiPoint, 4326)
);
CREATE INDEX data_point_multi_geom_index ON data_point_multi USING GIST (geometry);

DROP TABLE IF EXISTS "public"."data_pointz_multi";
CREATE TABLE data_pointz_multi (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(MultiPointZ, 4326)
);
CREATE INDEX data_pointz_multi_geom_index ON data_pointz_multi USING GIST (geometry);



DROP TABLE IF EXISTS "public"."data_line_multi";
CREATE TABLE data_line_multi (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(MultiLineString, 4326)
);
CREATE INDEX data_line_multi_geom_index ON data_line_multi USING GIST (geometry);

DROP TABLE IF EXISTS "public"."data_linez_multi";
CREATE TABLE data_linez_multi (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(MultiLineStringZ, 4326)
);
CREATE INDEX data_linez_multi_geom_index ON data_linez_multi USING GIST (geometry);



DROP TABLE IF EXISTS "public"."data_polygon_multi";
CREATE TABLE data_polygon_multi (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(MultiPolygon, 4326)
);
CREATE INDEX data_polygon_multi_geom_index ON data_polygon_multi USING GIST (geometry);

DROP TABLE IF EXISTS "public"."data_polygonz_multi";
CREATE TABLE data_polygonz_multi (
  id SERIAL PRIMARY KEY,
  geometry GEOMETRY(MultiPolygonZ, 4326)
);
CREATE INDEX data_polygonz_multi_geom_index ON data_polygonz_multi USING GIST (geometry);

3. 原本采用PostGIS的空间分析函数ST_Equals(geo1,geo2)来比较几何相等，但具有Z坐标后无法正确比较，即采用ST_Equals比较"POINT (1 1 1)"和“POINT (1 1 0)”是同一个几何

二 如何有效比较2D/3D两种维度下的几何是否相同

我将我尝试过的方案划分为了两种：一种是使用PostGIS空间分析函数进行的比较；另外一种是不使用空间分析函数，也就是直接将几何转为WKT字符串，进行字符串的比较，这种方案首先说可以，但是效率低、开销大，不建议使用。

1 利用空间分析函数的比较方式

涉及到的空间分析函数参考下表：

函数	名称	描述
ST_OrderingEquals(geo1, geo2)	Exactly Equal	通过逐点比较两个几何体来确定精确的相等性，以确保它们在位置上相同。
ST_Equals(geo1, geo2)	Spatially Equal	精确相等（Exactly Equal）并没有考虑几何图形的空间性质。有一个函数，恰当地命名为ST_Equals，可用于测试几何图形的空间相等性或等效性。
geo1 ~= geo2	Bounds Equal	在最坏的情况下，精确相等需要对几何体中的每个顶点进行比较，以确定相等。这可能很慢，并且可能不适合比较大量的几何形状。为了更快地进行比较，提供了相等边界运算符~=。这只对边界框（矩形）进行操作，确保几何图形占据相同的二维范围，但不一定占据相同的空间。

这部分官方提供了名为“Equality”可比较性的专题说明，建议浏览参考。
下列测试以Point、PointZ、MultiPoint、MultiPointZ等进行，测试数据如下：

-- Point
-- id 1
INSERT INTO data_point (geometry)
VALUES (ST_SetSRID(ST_MakePoint(1, 1), 4326));
-- id 2
INSERT INTO data_point (geometry)
VALUES (ST_SetSRID(ST_MakePoint(1, 1), 4326));
-- id 3
INSERT INTO data_point (geometry)
VALUES (ST_SetSRID(ST_MakePoint(1, 2), 4326));

-- PointZ
-- id 1
INSERT INTO data_pointz (geometry)
VALUES (ST_SetSRID(ST_MakePoint(1, 1, 1), 4326));
-- id 2
INSERT INTO data_pointz (geometry)
VALUES (ST_SetSRID(ST_MakePoint(1, 1, 0), 4326));
-- id 3
INSERT INTO data_pointz (geometry)
VALUES (ST_SetSRID(ST_MakePoint(1, 2, 0), 4326));

-- MultiPoint
-- id 1
INSERT INTO data_point_multi (geometry)
VALUES (ST_GeomFromText('MULTIPOINT ((1 1),(2 2))',4326));
-- id 2
INSERT INTO data_point_multi (geometry)
VALUES (ST_GeomFromText('MULTIPOINT ((2 2),(1 1))',4326));
-- id 3
INSERT INTO data_point_multi (geometry)
VALUES (ST_GeomFromText('MULTIPOINT ((2 2),(1 2))',4326));

-- MultiPointZ
-- id 1
INSERT INTO data_pointz_multi (geometry)
VALUES (ST_GeomFromText('MULTIPOINT ((1 1 1),(2 2 2))',4326));
-- id 2
INSERT INTO data_pointz_multi (geometry)
VALUES (ST_GeomFromText('MULTIPOINT ((2 2 2),(1 1 1))',4326));
-- id 3
INSERT INTO data_pointz_multi (geometry)
VALUES (ST_GeomFromText('MULTIPOINT ((2 2 2),(1 1 0))',4326));
-- id 4
INSERT INTO data_pointz_multi (geometry)
VALUES (ST_GeomFromText('MULTIPOINT ((1 1 1),(3 3 3))',4326));

（1）ST_Equals(geo1, geo2)

考虑几何图形的空间性质，可用于测试几何图形的空间相等性或等效性，更符合我们对几何相等的直观理解（即包含相同的区域）。需要注意，无论是绘制多边形的方向、定义多边形的起点，还是使用的点数，在这里都不重要，重要的是几何包含相同的空间。
下列的示例将对一些特别的地方进行说明：

a.忽略Z维度（特别注意！特别注意！特别注意！）

SELECT id,ST_AsText(geometry) AS wkt
FROM data_pointz
WHERE 
ST_Equals(
	geometry, 
	ST_GeomFromText(
		'POINT (1 1 1)',
		 4326
	)
)

注意！虽然data_pointz.1和data_pointz.2的Z不相同，但仍然会判定相等，因为ST_Equals函数仅仅对二维比较

b.忽略几何构件顺序

SELECT id,ST_AsText(geometry) AS wkt
FROM data_point_multi
WHERE 
ST_Equals(
	geometry, 
	ST_GeomFromText('MULTIPOINT ((1 1),(2 2))',4326)
);

注意，data_point_multi.1和data_point_multi.2两个MultiPoint都是由Point(1 1)和Point(2 2)组成的，只不过顺序不同，但依然被判定相等，同理在LineString、Polygon中也是一样的。

（2）ST_OrderingEquals(geo1, geo2)

ST_OrderingEquals函数则比ST_Equals函数更加“苛刻”，无论是绘制多边形的方向、定义多边形的起点，还是使用的点数，比较时都会被考虑。

a.Z维度的有效比较

-- id 1
INSERT INTO data_pointz (geometry)
VALUES (ST_SetSRID(ST_MakePoint(1, 1, 1), 4326));
-- id 2
INSERT INTO data_pointz (geometry)
VALUES (ST_SetSRID(ST_MakePoint(1, 1, 0), 4326));
-- id 3
INSERT INTO data_pointz (geometry)
VALUES (ST_SetSRID(ST_MakePoint(1, 2, 0), 4326));

SELECT id,ST_AsText(geometry) AS wkt
FROM data_pointz
WHERE 
ST_OrderingEquals(
	geometry, 
	ST_GeomFromText(
		'POINT (1 1 1)',
		 4326
	)
)

b.几何构建顺序的有效比较

-- id 1
INSERT INTO data_pointz_multi (geometry)
VALUES (ST_GeomFromText('MULTIPOINT ((1 1 1),(2 2 2))',4326));
-- id 2
INSERT INTO data_pointz_multi (geometry)
VALUES (ST_GeomFromText('MULTIPOINT ((2 2 2),(1 1 1))',4326));
-- id 3
INSERT INTO data_pointz_multi (geometry)
VALUES (ST_GeomFromText('MULTIPOINT ((2 2 2),(1 1 0))',4326));
-- id 4
INSERT INTO data_pointz_multi (geometry)
VALUES (ST_GeomFromText('MULTIPOINT ((1 1 1),(3 3 3))',4326));

SELECT id,ST_AsText(geometry) AS wkt
FROM data_pointz_multi
WHERE 
ST_OrderingEquals(
	geometry, 
	ST_GeomFromText('MULTIPOINT ((2 2 2),(1 1 1))',4326)
);

2 不太常规的比较方式

（1）转WKT进行纯字符串层面的比较

-- 比较方式4：字符串比较
SELECT id,ST_AsText(geometry) AS wkt
FROM data_point
WHERE 
ST_AsText(geometry) = ST_AsText(ST_GeomFromText('POINT(1 1)'));

-- 比较方式4：字符串比较
SELECT id,ST_AsText(geometry) AS wkt
FROM data_pointz
WHERE 
ST_AsText(geometry) = ST_AsText(ST_GeomFromText('POINT(1 1 1)'));

（2）利用PostGIS的距离求算函数判断是否有距离为0的几何存在，如果有则代表已经存在；

SELECT id,ST_AsText(geometry) AS wkt
FROM data_pointz
WHERE 
(ST_Distance(geometry,ST_GeomFromText('POINT (1 1 1)',4326))) = 0;

SELECT id,ST_AsText(geometry) AS wkt
FROM data_pointz
WHERE 
(ST_3DDistance(geometry,ST_GeomFromText('POINT (1 1 1)',4326))) = 0;

三 总结

类似插入记录时，如果几何存在则返回已存在记录id的需求还是很常见的，但是需要结合具体的业务场景以及数据情况（2D/3D/…)来判断采用什么类型的空间分析函数。
ST_OrderingEquals足够精确，无论是绘制多边形的方向、定义多边形的起点，还是使用的点数，比较时都会被考虑。我用GeoJson举一个例子：

{
    "type": "FeatureCollection",
    "name": "水位监测",
    "crs": {},
    "features": [
        {
            "type": "Feature",
            "properties": {
                "time": "time 1",
                "waterLevel": [
                    1.1,
                    2.2,
                    3.3
                ]
            },
            "geometry": {
                "type": "MultiPoint",
                "coordinates": [
                    [1,1],
                    [2,2],
                    [3,3]
                ]
            }
        },
        {
            "type": "Feature",
            "properties": {
                "time": "time 2",
                "waterLevel": [
                    1.11,
                    2.22,
                    3.33
                ]
            },
            "geometry": {
                "type": "MultiPoint",
                "coordinates": [
                    [1,1],
                    [2,2],
                    [3,3]
                ]
            }
        }
    ]
}

数据包含两个要素，要素为MultiPoint类型，包含三个坐标点，每个坐标点为水位监测仪器的的空间位置，属性“waterLevel”则对应三个点的具体水位值。类似这样的场景，我们就不能仅仅考虑几何是否包含（占据）相同的空间。
假设某些场景中并不需要考虑Z坐标以及几何的构建顺序，则ST_Equals就是最佳选择，因为它的效率优于前者，