STIntersect 表示每个点都与每个区域相交 - 需要确认
STIntersect says every point is intersecting every zone - need to confirm
我正在为一个行业团体实施车辆安全计划。该小组提供了 3k+ 地理区域,覆盖了该小组成员的站点。提供的文件是一个 shapefile。我的同事把它提取出来,转换成4326,上传到一个新的数据库中。我使用 IsValidDetailed()
检查了区域。 7 是无效的,但当 MakeValid()
是 运行 时他们自己纠正了。地理 table SaferTogether 已创建空间索引,仅使用自动。
我最初使用以下方法进行测试。这将 运行 每个区域的每个点,我在一个多小时后得到了结果。
INSERT INTO @Intersects (LogId, ogr_fid, NodeId)
SELECT p.LogId, s.ogr_fid, p.NodeId
FROM LogPos p,
[gis]..safertogether as s
WHERE p.MyPoint.STIntersects(s.geog4326) = 1
and AssembledTime between @PeriodStart and @PeriodEnd
我现在正在测试 27 个日志点的单次行程,我注意到我有 81162 个结果这大致匹配 27 个点 x 3k+ 区域。我不相信每个区域都覆盖每个点。下面的代码是对一趟日志的测试。效率不高。
UPDATE #Logs
SET IsIntersect = 0
SET @Message = 'Detecting Log intersects with Zone ' + CONVERT(VARCHAR(10), @GeographyId) + ':' + @GeographyName
EXEC usp_Log @ProcessName, @Message, @SessionName, @LogInfo, @Section
UPDATE #Logs
SET IsIntersect = 1
WHERE @Geography.STIntersects(MyPoint) = 1
SELECT @rc = COUNT(1)
FROM #Logs WHERE IsIntersect = 1
SET @Message = 'Trip ' + CONVERT(VARCHAR(10), @TripLogId) + ' intersecting points with ' + CONVERT(VARCHAR(10), @GeographyId) + ':' + @GeographyName + ': ' + CONVERT(varchar(10), @rc)
EXEC usp_Log @ProcessName, @Message, @SessionName, @LogDebug, @Section
SET @Message = 'Save intersects with ' + CONVERT(VARCHAR(10), @GeographyId) + ':' + @GeographyName
EXEC usp_Log @ProcessName, @Message, @SessionName, @LogInfo, @Section
INSERT INTO @Intersects (LogId, NodeId, ogr_fid)
SELECT
LogId,
NodeId,
@GeographyId
FROM #Logs
WHERE IsIntersect = 1
我在这里寻找的是一种可视化区域与点的方法,这样我就可以判断它们是否损坏。我不知道如何比较内置空间结果选项卡中的字段。
问题可能出在将 shapefile(平面地图)数据转换为 SQL 服务器地理(球形)模型时。这些很容易出错并得到倒置的(免费的)多边形,即描述所有地球区域的多边形除了预期的形状。
另请参阅几个类似的案例,它们与 Google BigQuery 一起使用,它使用面向球面的多边形的相似语义:
您可以通过计算多边形的面积来检查多边形是否被正确摄取到数据库中 - 对于倒置的多边形来说,它会很大(大约 5e14 m^2),而对于正确摄取的多边形,您应该会得到合理的结果。
SQL 服务器有 ReorientObject()
反转多边形方向的方法。但是要小心,因为方向与有效性相互作用 - 所以最好在 MakeValid()
之前使用 ReorientObject()
。最重要的当然是修复摄取方法。
或者,您可以使用使用平面地图语义的 SQL 服务器 Geometry
(而不是 Geography
)类型。使用平面数据存在一些问题,但假设您没有交叉 anti-meridian,并且对哪种投影适合哪种计算有基本的了解,您通常应该没问题。
或者切换到使用球形模式但可以正确解释和转换平面数据的数据库。例如。 BigQuery 了解 GeoJson 几何图形是平面的,并且在将它们转换为内部球形地理模型时会忽略它们的方向,因此如果将几何图形加载为 GeoJson 字段,则永远不会遇到此问题。
以下是我使用的最终脚本,以防有人需要帮助。
@Shapes 是一个 table 变量。它被插入填充,然后由游标循环。
@Shape 的每个游标循环都经过有效性测试。如果无效,它有 MakeValid()
运行。然后测试它是否是 MULTIPOLYGON 。如果是,则将其分解为@SubShapes。如果是单个形状,则将该形状分配给第一个@SubShape。
每个@SubShape 都经过验证和面积测试。如果需要,MakeValid()
和 ReorientObject()
是 运行。
在@Shape 的所有@SubShapes 都是运行 之后,如果@Shape 是一个MULTIPOLYGON,它会被放回一起。这是通过 WellKnownText 的字符串连接和操作完成的。这是我最不满意的一点。感觉像是地理上的落差API.
一旦处理完所有的@Shapes,它们就会被写回真实的table
usp_Log_Error 是一个 in-house 记录过程。欢迎评论出来
SET NOCOUNT ON
DECLARE @ProcessName varchar(50) = 'ST_Geography_Validation'
DECLARE @SessionName varchar(50) = 'ST_GV_' + CONVERT(VARCHAR(19), GETDATE(), 126)
DECLARE @Message varchar(1000)
DECLARE @Id int
DECLARE @Name varchar(254)
DECLARE @Shape geography
DECLARE @ValidationDetail varchar(200)
DECLARE @ShapeCount int
DECLARE @SubValidationDetail varchar(200)
DECLARE @SubShape geography
Declare @Area float
DECLARE @AreaStr varchar(20)
DECLARE @i int
DECLARE @IsMulti int = 0
DECLARE @MultiWKT varchar(MAX)
DECLARE @SingleWKT varchar(MAX)
DECLARE @Shapes TABLE
(
Id int,
Name varchar(254),
Shape geography
)
INSERT INTO @Shapes (Id, Name, Shape)
SELECT
ogr_fid,
name,
geog4326
FROM safertogether2
DECLARE @SubShapes TABLE
(
SubId int PRIMARY KEY,
SubShape geography,
Processed bit default 0
)
DECLARE _c CURSOR FOR
SELECT
s.Id,
s.name,
s.Shape,
0,
'',
''
FROM @Shapes s
OPEN _c
FETCH NEXT FROM _c INTO @Id, @Name, @Shape, @IsMulti, @MultiWKT, @SingleWKT
WHILE @@FETCH_STATUS = 0
BEGIN
DELETE FROM @SubShapes
SET @ValidationDetail = @Shape.IsValidDetailed() -- This is for the base shape
IF @ValidationDetail <> '24400: Valid'
BEGIN
SET @Shape = @Shape.MakeValid()
END
SET @ShapeCount = @Shape.STNumGeometries()
IF @ShapeCount > 1
BEGIN
RAISERROR('%i:%s: is MULTIPOLYGON. Break into subshapes.', 0, 1, @Id, @Name) WITH NOWAIT
SET @IsMulti = 2
INSERT INTO @SubShapes(SubId, SubShape)
SELECT
smg.Id,
smg.Geog
FROM [dbo].uf_SplitMultiGeography(@Shape) smg
END
ELSE IF @ShapeCount = 1
BEGIN
RAISERROR('%i:%s: is SINGLE. Insert Shape whole as Id 1.', 0, 1, @Id, @Name) WITH NOWAIT
SET @IsMulti = 1
INSERT INTO @SubShapes(SubId, SubShape)
VALUES (1, @Shape)
END
ELSE
BEGIN
RAISERROR('%i:%s: is NONE.', 0, 1, @Id, @Name) WITH NOWAIT
SET @IsMulti = 0
END
SET @i = 1
WHILE @i <= @ShapeCount
BEGIN
SELECT @SubShape = s.SubShape
FROM @SubShapes s
WHERE SubId = @i
SET @SubValidationDetail = @SubShape.IsValidDetailed()
RAISERROR('%i:%s:%i: .IsValidDetailed(): %s', 0, 1, @Id, @Name, @i, @SubValidationDetail) WITH NOWAIT
IF @SubValidationDetail <> '24400: Valid'
BEGIN
RAISERROR('%i:%s:%i: .MakeValid()', 0, 1, @Id, @Name, @i) WITH NOWAIT
SET @SubShape = @SubShape.MakeValid() -- This is for the subshape
END
BEGIN TRY
SET @Area = @SubShape.STArea()
END TRY
BEGIN CATCH
SET @Message = 'While getting the Area of a subshape'
EXEC usp_Log_Error @ProcessName, @SessionName, @Message
END CATCH
SET @AreaStr = LTRIM(STR(@Area, 20, 5))
RAISERROR('%i:%s:%i: .STArea(): Pre: %s', 0, 1, @Id, @Name, @i, @AreaStr) WITH NOWAIT
IF @Area > 510000000000000
BEGIN
RAISERROR('%i:%s:%i: .ReorientObject()', 0, 1, @Id, @Name, @i) WITH NOWAIT
SET @SubShape = @SubShape.ReorientObject()
SET @Area = @SubShape.STArea()
SET @AreaStr = LTRIM(STR(@Area, 20, 5))
RAISERROR('%i:%s:%i: .STArea(): Post: %s', 0, 1, @Id, @Name, @i, @AreaStr) WITH NOWAIT
END
UPDATE @SubShapes
SET SubShape = @SubShape
WHERE SubId = @Id
IF @IsMulti = 2
BEGIN
RAISERROR('%i:%s:%i: MULTIPOLYGON member. Process WKT for recombination', 0, 1, @Id, @Name, @i) WITH NOWAIT
BEGIN TRY
SET @SingleWKT = @SubShape.STAsText()
END TRY
BEGIN CATCH
SET @Message = 'While getting the WKT of a subshape'
EXEC usp_Log_Error @ProcessName, @SessionName, @Message
END CATCH
SET @SingleWKT = REPLACE(@SingleWKT, 'POLYGON', '')
IF (@i <> @ShapeCount)
BEGIN
SET @SingleWKT = REPLACE(@SingleWKT, '))', ')),')
END
SET @MultiWKT = @MultiWKT + @SingleWKT
END
SET @i = @i + 1
END
IF @IsMulti = 2
BEGIN
RAISERROR('%i:%s: Is a MULTIPOLYGON. Recreate from updated WKT.', 0, 1, @Id, @Name, @i) WITH NOWAIT
SET @MultiWKT = 'MULTIPOLYGON(' + @MultiWKT + ')'
BEGIN TRY
SET @Shape = geography::STGeomFromText(@MultiWKT, 4326)
END TRY
BEGIN CATCH
SET @Message = 'While creating a multipolygon'
EXEC usp_Log_Error @ProcessName, @SessionName, @Message
END CATCH
END
ELSE IF @IsMulti = 1
BEGIN
RAISERROR('%i:%s: Is SINGLE. Set Shape to subshape', 0, 1, @Id, @Name, @i) WITH NOWAIT
SET @Shape = @SubShape
END
ELSE
BEGIN
RAISERROR('%i:%s: Is NONE. Set Shape to none', 0, 1, @Id, @Name, @i) WITH NOWAIT
END
RAISERROR('%i:%s: Update temp table', 0, 1, @Id, @Name) WITH NOWAIT
UPDATE @Shapes
SET Shape = @Shape
WHERE Id = @Id
FETCH NEXT FROM _c INTO @Id, @Name, @Shape, @IsMulti, @MultiWKT, @SingleWKT
END
CLOSE _c
DEALLOCATE _c
select
st.ogr_fid,
st.name,
st.geog4326 as ShapePre,
s.Shape as ShapePost,
CASE
WHEN st.geog4326.IsValidDetailed() = '24400: Valid' then st.geog4326.STArea()
ELSE st.geog4326.MakeValid().STArea()
END as AreaPre,
CASE
WHEN s.Shape.IsValidDetailed() = '24400: Valid' then s.Shape.STArea()
ELSE s.Shape.MakeValid().STArea()
END as AreaPost
from safertogether2 st
join @Shapes s ON s.Id = st.ogr_fid
RAISERROR('Update back to safertogether table', 0, 1) WITH NOWAIT
UPDATE st
SET geog4326 = s.Shape
FROM @Shapes s
JOIN safertogether2 st ON s.Id = st.ogr_fid
RAISERROR('End of Script', 0, 1, @Id, @Name) WITH NOWAIT
我正在为一个行业团体实施车辆安全计划。该小组提供了 3k+ 地理区域,覆盖了该小组成员的站点。提供的文件是一个 shapefile。我的同事把它提取出来,转换成4326,上传到一个新的数据库中。我使用 IsValidDetailed()
检查了区域。 7 是无效的,但当 MakeValid()
是 运行 时他们自己纠正了。地理 table SaferTogether 已创建空间索引,仅使用自动。
我最初使用以下方法进行测试。这将 运行 每个区域的每个点,我在一个多小时后得到了结果。
INSERT INTO @Intersects (LogId, ogr_fid, NodeId)
SELECT p.LogId, s.ogr_fid, p.NodeId
FROM LogPos p,
[gis]..safertogether as s
WHERE p.MyPoint.STIntersects(s.geog4326) = 1
and AssembledTime between @PeriodStart and @PeriodEnd
我现在正在测试 27 个日志点的单次行程,我注意到我有 81162 个结果这大致匹配 27 个点 x 3k+ 区域。我不相信每个区域都覆盖每个点。下面的代码是对一趟日志的测试。效率不高。
UPDATE #Logs
SET IsIntersect = 0
SET @Message = 'Detecting Log intersects with Zone ' + CONVERT(VARCHAR(10), @GeographyId) + ':' + @GeographyName
EXEC usp_Log @ProcessName, @Message, @SessionName, @LogInfo, @Section
UPDATE #Logs
SET IsIntersect = 1
WHERE @Geography.STIntersects(MyPoint) = 1
SELECT @rc = COUNT(1)
FROM #Logs WHERE IsIntersect = 1
SET @Message = 'Trip ' + CONVERT(VARCHAR(10), @TripLogId) + ' intersecting points with ' + CONVERT(VARCHAR(10), @GeographyId) + ':' + @GeographyName + ': ' + CONVERT(varchar(10), @rc)
EXEC usp_Log @ProcessName, @Message, @SessionName, @LogDebug, @Section
SET @Message = 'Save intersects with ' + CONVERT(VARCHAR(10), @GeographyId) + ':' + @GeographyName
EXEC usp_Log @ProcessName, @Message, @SessionName, @LogInfo, @Section
INSERT INTO @Intersects (LogId, NodeId, ogr_fid)
SELECT
LogId,
NodeId,
@GeographyId
FROM #Logs
WHERE IsIntersect = 1
我在这里寻找的是一种可视化区域与点的方法,这样我就可以判断它们是否损坏。我不知道如何比较内置空间结果选项卡中的字段。
问题可能出在将 shapefile(平面地图)数据转换为 SQL 服务器地理(球形)模型时。这些很容易出错并得到倒置的(免费的)多边形,即描述所有地球区域的多边形除了预期的形状。
另请参阅几个类似的案例,它们与 Google BigQuery 一起使用,它使用面向球面的多边形的相似语义:
您可以通过计算多边形的面积来检查多边形是否被正确摄取到数据库中 - 对于倒置的多边形来说,它会很大(大约 5e14 m^2),而对于正确摄取的多边形,您应该会得到合理的结果。
SQL 服务器有 ReorientObject()
反转多边形方向的方法。但是要小心,因为方向与有效性相互作用 - 所以最好在 MakeValid()
之前使用 ReorientObject()
。最重要的当然是修复摄取方法。
或者,您可以使用使用平面地图语义的 SQL 服务器 Geometry
(而不是 Geography
)类型。使用平面数据存在一些问题,但假设您没有交叉 anti-meridian,并且对哪种投影适合哪种计算有基本的了解,您通常应该没问题。
或者切换到使用球形模式但可以正确解释和转换平面数据的数据库。例如。 BigQuery 了解 GeoJson 几何图形是平面的,并且在将它们转换为内部球形地理模型时会忽略它们的方向,因此如果将几何图形加载为 GeoJson 字段,则永远不会遇到此问题。
以下是我使用的最终脚本,以防有人需要帮助。
@Shapes 是一个 table 变量。它被插入填充,然后由游标循环。
@Shape 的每个游标循环都经过有效性测试。如果无效,它有 MakeValid()
运行。然后测试它是否是 MULTIPOLYGON 。如果是,则将其分解为@SubShapes。如果是单个形状,则将该形状分配给第一个@SubShape。
每个@SubShape 都经过验证和面积测试。如果需要,MakeValid()
和 ReorientObject()
是 运行。
在@Shape 的所有@SubShapes 都是运行 之后,如果@Shape 是一个MULTIPOLYGON,它会被放回一起。这是通过 WellKnownText 的字符串连接和操作完成的。这是我最不满意的一点。感觉像是地理上的落差API.
一旦处理完所有的@Shapes,它们就会被写回真实的table
usp_Log_Error 是一个 in-house 记录过程。欢迎评论出来
SET NOCOUNT ON
DECLARE @ProcessName varchar(50) = 'ST_Geography_Validation'
DECLARE @SessionName varchar(50) = 'ST_GV_' + CONVERT(VARCHAR(19), GETDATE(), 126)
DECLARE @Message varchar(1000)
DECLARE @Id int
DECLARE @Name varchar(254)
DECLARE @Shape geography
DECLARE @ValidationDetail varchar(200)
DECLARE @ShapeCount int
DECLARE @SubValidationDetail varchar(200)
DECLARE @SubShape geography
Declare @Area float
DECLARE @AreaStr varchar(20)
DECLARE @i int
DECLARE @IsMulti int = 0
DECLARE @MultiWKT varchar(MAX)
DECLARE @SingleWKT varchar(MAX)
DECLARE @Shapes TABLE
(
Id int,
Name varchar(254),
Shape geography
)
INSERT INTO @Shapes (Id, Name, Shape)
SELECT
ogr_fid,
name,
geog4326
FROM safertogether2
DECLARE @SubShapes TABLE
(
SubId int PRIMARY KEY,
SubShape geography,
Processed bit default 0
)
DECLARE _c CURSOR FOR
SELECT
s.Id,
s.name,
s.Shape,
0,
'',
''
FROM @Shapes s
OPEN _c
FETCH NEXT FROM _c INTO @Id, @Name, @Shape, @IsMulti, @MultiWKT, @SingleWKT
WHILE @@FETCH_STATUS = 0
BEGIN
DELETE FROM @SubShapes
SET @ValidationDetail = @Shape.IsValidDetailed() -- This is for the base shape
IF @ValidationDetail <> '24400: Valid'
BEGIN
SET @Shape = @Shape.MakeValid()
END
SET @ShapeCount = @Shape.STNumGeometries()
IF @ShapeCount > 1
BEGIN
RAISERROR('%i:%s: is MULTIPOLYGON. Break into subshapes.', 0, 1, @Id, @Name) WITH NOWAIT
SET @IsMulti = 2
INSERT INTO @SubShapes(SubId, SubShape)
SELECT
smg.Id,
smg.Geog
FROM [dbo].uf_SplitMultiGeography(@Shape) smg
END
ELSE IF @ShapeCount = 1
BEGIN
RAISERROR('%i:%s: is SINGLE. Insert Shape whole as Id 1.', 0, 1, @Id, @Name) WITH NOWAIT
SET @IsMulti = 1
INSERT INTO @SubShapes(SubId, SubShape)
VALUES (1, @Shape)
END
ELSE
BEGIN
RAISERROR('%i:%s: is NONE.', 0, 1, @Id, @Name) WITH NOWAIT
SET @IsMulti = 0
END
SET @i = 1
WHILE @i <= @ShapeCount
BEGIN
SELECT @SubShape = s.SubShape
FROM @SubShapes s
WHERE SubId = @i
SET @SubValidationDetail = @SubShape.IsValidDetailed()
RAISERROR('%i:%s:%i: .IsValidDetailed(): %s', 0, 1, @Id, @Name, @i, @SubValidationDetail) WITH NOWAIT
IF @SubValidationDetail <> '24400: Valid'
BEGIN
RAISERROR('%i:%s:%i: .MakeValid()', 0, 1, @Id, @Name, @i) WITH NOWAIT
SET @SubShape = @SubShape.MakeValid() -- This is for the subshape
END
BEGIN TRY
SET @Area = @SubShape.STArea()
END TRY
BEGIN CATCH
SET @Message = 'While getting the Area of a subshape'
EXEC usp_Log_Error @ProcessName, @SessionName, @Message
END CATCH
SET @AreaStr = LTRIM(STR(@Area, 20, 5))
RAISERROR('%i:%s:%i: .STArea(): Pre: %s', 0, 1, @Id, @Name, @i, @AreaStr) WITH NOWAIT
IF @Area > 510000000000000
BEGIN
RAISERROR('%i:%s:%i: .ReorientObject()', 0, 1, @Id, @Name, @i) WITH NOWAIT
SET @SubShape = @SubShape.ReorientObject()
SET @Area = @SubShape.STArea()
SET @AreaStr = LTRIM(STR(@Area, 20, 5))
RAISERROR('%i:%s:%i: .STArea(): Post: %s', 0, 1, @Id, @Name, @i, @AreaStr) WITH NOWAIT
END
UPDATE @SubShapes
SET SubShape = @SubShape
WHERE SubId = @Id
IF @IsMulti = 2
BEGIN
RAISERROR('%i:%s:%i: MULTIPOLYGON member. Process WKT for recombination', 0, 1, @Id, @Name, @i) WITH NOWAIT
BEGIN TRY
SET @SingleWKT = @SubShape.STAsText()
END TRY
BEGIN CATCH
SET @Message = 'While getting the WKT of a subshape'
EXEC usp_Log_Error @ProcessName, @SessionName, @Message
END CATCH
SET @SingleWKT = REPLACE(@SingleWKT, 'POLYGON', '')
IF (@i <> @ShapeCount)
BEGIN
SET @SingleWKT = REPLACE(@SingleWKT, '))', ')),')
END
SET @MultiWKT = @MultiWKT + @SingleWKT
END
SET @i = @i + 1
END
IF @IsMulti = 2
BEGIN
RAISERROR('%i:%s: Is a MULTIPOLYGON. Recreate from updated WKT.', 0, 1, @Id, @Name, @i) WITH NOWAIT
SET @MultiWKT = 'MULTIPOLYGON(' + @MultiWKT + ')'
BEGIN TRY
SET @Shape = geography::STGeomFromText(@MultiWKT, 4326)
END TRY
BEGIN CATCH
SET @Message = 'While creating a multipolygon'
EXEC usp_Log_Error @ProcessName, @SessionName, @Message
END CATCH
END
ELSE IF @IsMulti = 1
BEGIN
RAISERROR('%i:%s: Is SINGLE. Set Shape to subshape', 0, 1, @Id, @Name, @i) WITH NOWAIT
SET @Shape = @SubShape
END
ELSE
BEGIN
RAISERROR('%i:%s: Is NONE. Set Shape to none', 0, 1, @Id, @Name, @i) WITH NOWAIT
END
RAISERROR('%i:%s: Update temp table', 0, 1, @Id, @Name) WITH NOWAIT
UPDATE @Shapes
SET Shape = @Shape
WHERE Id = @Id
FETCH NEXT FROM _c INTO @Id, @Name, @Shape, @IsMulti, @MultiWKT, @SingleWKT
END
CLOSE _c
DEALLOCATE _c
select
st.ogr_fid,
st.name,
st.geog4326 as ShapePre,
s.Shape as ShapePost,
CASE
WHEN st.geog4326.IsValidDetailed() = '24400: Valid' then st.geog4326.STArea()
ELSE st.geog4326.MakeValid().STArea()
END as AreaPre,
CASE
WHEN s.Shape.IsValidDetailed() = '24400: Valid' then s.Shape.STArea()
ELSE s.Shape.MakeValid().STArea()
END as AreaPost
from safertogether2 st
join @Shapes s ON s.Id = st.ogr_fid
RAISERROR('Update back to safertogether table', 0, 1) WITH NOWAIT
UPDATE st
SET geog4326 = s.Shape
FROM @Shapes s
JOIN safertogether2 st ON s.Id = st.ogr_fid
RAISERROR('End of Script', 0, 1, @Id, @Name) WITH NOWAIT