Skip to content

Commit 87b3fd2

Browse files
berttbertt
committed
filter boundingboxes based on original projection (not epsg:4326)
1 parent 45c749d commit 87b3fd2

8 files changed

Lines changed: 127 additions & 82 deletions

File tree

‎src/b3dm.tileset.tests/GeometryRepositoryWhereTests.cs‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
using NUnit.Framework;
2+
using Wkx;
3+
4+
namespace B3dm.Tileset.Tests;
5+
6+
public class GeometryRepositoryWhereTests
7+
{
8+
[Test]
9+
public void GetWhere_2D_UseSourceEpsgDirectly()
10+
{
11+
var from = new Point(100000.0, 400000.0);
12+
var to = new Point(200000.0, 500000.0);
13+
14+
var sql = GeometryRepository.GetWhere("geom", from, to, "", 5698);
15+
16+
Assert.That(sql, Does.Contain("5698"));
17+
Assert.That(sql, Does.Not.Contain("4326"));
18+
Assert.That(sql, Does.Not.Contain("ST_Transform"));
19+
Assert.That(sql, Does.Contain("ST_MakeEnvelope"));
20+
}
21+
22+
[Test]
23+
public void GetWhere_2D_Epsg4326_UseSourceEpsgDirectly()
24+
{
25+
var from = new Point(-75.8, 38.4);
26+
var to = new Point(-75.0, 39.8);
27+
28+
var sql = GeometryRepository.GetWhere("geom", from, to, "", 4326);
29+
30+
Assert.That(sql, Does.Contain("4326"));
31+
Assert.That(sql, Does.Not.Contain("ST_Transform"));
32+
Assert.That(sql, Does.Contain("ST_MakeEnvelope"));
33+
}
34+
35+
[Test]
36+
public void GetWhere_3D_UseSourceEpsgDirectly()
37+
{
38+
var from = new Point(100000.0, 400000.0, 200.0);
39+
var to = new Point(200000.0, 500000.0, 300.0);
40+
41+
var sql = GeometryRepository.GetWhere("geom", from, to, "", 5698);
42+
43+
Assert.That(sql, Does.Contain("5698"));
44+
Assert.That(sql, Does.Not.Contain("4979"));
45+
Assert.That(sql, Does.Not.Contain("4326"));
46+
Assert.That(sql, Does.Not.Contain("ST_Transform"));
47+
Assert.That(sql, Does.Contain("ST_3DMakeBox"));
48+
}
49+
50+
[Test]
51+
public void GetWhere_3D_ContainsCentroidAndSrid()
52+
{
53+
var from = new Point(100000.0, 400000.0, 200.0);
54+
var to = new Point(200000.0, 500000.0, 300.0);
55+
56+
var sql = GeometryRepository.GetWhere("geom", from, to, "", 5698);
57+
58+
Assert.That(sql, Does.Contain("st_setsrid"));
59+
Assert.That(sql, Does.Contain("ST_3DIntersects"));
60+
}
61+
}

‎src/b3dm.tileset/BoundingBoxRepository.cs‎

Lines changed: 27 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
1-
using System.Data;
1+
using System;
2+
using System.Data;
23
using Wkx;
34

45
namespace B3dm.Tileset;
@@ -16,6 +17,31 @@ public static (BoundingBox bbox, double zmin, double zmax) GetBoundingBoxForTabl
1617
return bbox3d;
1718
}
1819

20+
public static (BoundingBox bbox, double zmin, double zmax) GetBoundingBoxAs4979(IDbConnection conn,(BoundingBox bbox, double zmin, double zmax) bboxTable, int sourceEpsg)
21+
{
22+
if (sourceEpsg == 4979 || sourceEpsg == 4326) {
23+
return bboxTable;
24+
}
25+
var bbox = bboxTable.bbox;
26+
var sqlBounds = FormattableString.Invariant($@"SELECT st_xmin(geom1),st_ymin(geom1), st_xmax(geom1), st_ymax(geom1), st_zmin(geom1), st_zmax(geom1)
27+
FROM (
28+
SELECT ST_3DExtent(
29+
ST_Transform(
30+
ST_SetSRID(
31+
ST_3DMakeBox(
32+
ST_MakePoint({bbox.XMin}, {bbox.YMin}, {bboxTable.zmin}),
33+
ST_MakePoint({bbox.XMax}, {bbox.YMax}, {bboxTable.zmax})
34+
),
35+
{sourceEpsg}
36+
),
37+
4979
38+
)
39+
) AS geom1
40+
) AS t");
41+
return GetBounds(conn, sqlBounds);
42+
}
43+
44+
1945
private static (BoundingBox, double, double) GetBounds(IDbConnection conn, string sql)
2046
{
2147
conn.Open();

‎src/b3dm.tileset/FeatureCountRepository.cs‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,10 +7,10 @@ namespace B3dm.Tileset;
77

88
public static class FeatureCountRepository
99
{
10-
public static int CountFeaturesInBox(NpgsqlConnection conn, string geometry_table, string geometry_column, Point from, Point to, string query, int source_epsg, bool keepProjection = false, HashSet<string> excludeHashes = null)
10+
public static int CountFeaturesInBox(NpgsqlConnection conn, string geometry_table, string geometry_column, Point from, Point to, string query, int source_epsg, HashSet<string> excludeHashes = null)
1111
{
1212
var select = $"COUNT({geometry_column})";
13-
var where = GeometryRepository.GetWhere(geometry_column, from, to, query, source_epsg, keepProjection);
13+
var where = GeometryRepository.GetWhere(geometry_column, from, to, query, source_epsg);
1414

1515
// Add hash exclusion filter using parameterized query
1616
if (excludeHashes != null && excludeHashes.Count > 0) {

‎src/b3dm.tileset/GeometryRepository.cs‎

Lines changed: 7 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -13,48 +13,6 @@ namespace B3dm.Tileset;
1313

1414
public static class GeometryRepository
1515
{
16-
public static HashSet<string> FilterHashesByEnvelope(NpgsqlConnection conn, string tableName, string geometryColumn, BoundingBox bbox, int source_epsg, HashSet<string> geometryHashes, bool keepProjection)
17-
{
18-
if (geometryHashes.Count == 0) {
19-
return new HashSet<string>();
20-
}
21-
22-
var filteredHashes = new HashSet<string>();
23-
24-
conn.Open();
25-
try {
26-
var envelope = keepProjection ?
27-
$"ST_MakeEnvelope(@xmin, @ymin, @xmax, @ymax, {source_epsg})" :
28-
$"ST_Transform(ST_MakeEnvelope(@xmin, @ymin, @xmax, @ymax, 4326), {source_epsg})";
29-
30-
var query = $@"
31-
SELECT MD5(ST_AsBinary({geometryColumn})::text) as geom_hash
32-
FROM {tableName}
33-
WHERE MD5(ST_AsBinary({geometryColumn})::text) = ANY(@hashes)
34-
AND ST_Within(
35-
ST_Centroid(ST_Envelope({geometryColumn})),
36-
{envelope}
37-
)";
38-
39-
using var cmd = new NpgsqlCommand(query, conn);
40-
cmd.Parameters.AddWithValue("hashes", geometryHashes.ToArray());
41-
cmd.Parameters.AddWithValue("xmin", bbox.XMin);
42-
cmd.Parameters.AddWithValue("ymin", bbox.YMin);
43-
cmd.Parameters.AddWithValue("xmax", bbox.XMax);
44-
cmd.Parameters.AddWithValue("ymax", bbox.YMax);
45-
46-
using var reader = cmd.ExecuteReader();
47-
while (reader.Read()) {
48-
var hash = reader.GetString(0);
49-
filteredHashes.Add(hash);
50-
}
51-
52-
return filteredHashes;
53-
}
54-
finally {
55-
conn.Close();
56-
}
57-
}
5816

5917
/// <summary>
6018
/// Returns double array with 6 bounding box coordinates, xmin, ymin, xmax, ymax, zmin, zmax
@@ -85,13 +43,13 @@ public static double[] GetGeometriesBoundingBox(NpgsqlConnection conn, string ge
8543
}
8644
}
8745

88-
public static List<GeometryRecord> GetGeometrySubset(NpgsqlConnection conn, string geometry_table, string geometry_column, double[] bbox, int source_epsg, int target_srs, string shaderColumn = "", string attributesColumns = "", string query = "", string radiusColumn = "", bool keepProjection = false, HashSet<string> excludeHashes = null, int? maxFeatures = null, SortBy sortBy = SortBy.AREA)
46+
public static List<GeometryRecord> GetGeometrySubset(NpgsqlConnection conn, string geometry_table, string geometry_column, double[] bbox, int source_epsg, int target_srs, string shaderColumn = "", string attributesColumns = "", string query = "", string radiusColumn = "", HashSet<string> excludeHashes = null, int? maxFeatures = null, SortBy sortBy = SortBy.AREA)
8947
{
9048
var sqlselect = GetSqlSelect(geometry_column, shaderColumn, attributesColumns, radiusColumn, target_srs);
9149
var sqlFrom = "FROM " + geometry_table;
9250
var points = GetPoints(bbox);
9351

94-
var sqlWhere = GetWhere(geometry_column, points.fromPoint, points.toPoint, query, source_epsg, keepProjection);
52+
var sqlWhere = GetWhere(geometry_column, points.fromPoint, points.toPoint, query, source_epsg);
9553

9654
// Add hash exclusion filter using parameterized query
9755
if (excludeHashes != null && excludeHashes.Count > 0) {
@@ -117,7 +75,7 @@ public static List<GeometryRecord> GetGeometrySubset(NpgsqlConnection conn, stri
11775
}
11876
}
11977

120-
public static string GetWhere(string geometry_column, Point from, Point to, string query, int source_epsg, bool keepProjection)
78+
public static string GetWhere(string geometry_column, Point from, Point to, string query, int source_epsg)
12179
{
12280
var fromX = from.X.Value.ToString(CultureInfo.InvariantCulture);
12381
var fromY = from.Y.Value.ToString(CultureInfo.InvariantCulture);
@@ -128,22 +86,14 @@ public static string GetWhere(string geometry_column, Point from, Point to, stri
12886
var where = "";
12987

13088
if (!hasZ) {
131-
where = keepProjection ?
132-
$"ST_Centroid(ST_Envelope({geometry_column})) && ST_MakeEnvelope({fromX}, {fromY}, {toX}, {toY}, {source_epsg}) {query}" :
133-
$"ST_Centroid(ST_Envelope({geometry_column})) && st_transform(ST_MakeEnvelope({fromX}, {fromY}, {toX}, {toY}, 4326), {source_epsg}) {query}";
89+
where = $"ST_Centroid(ST_Envelope({geometry_column})) && ST_MakeEnvelope({fromX}, {fromY}, {toX}, {toY}, {source_epsg}) {query}";
13490
}
13591
else {
136-
var fromBox = keepProjection ?
137-
$"st_setsrid(ST_MakePoint({fromX}, {fromY}, {from.Z.Value.ToString(CultureInfo.InvariantCulture)}), {source_epsg})" :
138-
$"st_setsrid(ST_MakePoint({fromX}, {fromY}, {from.Z.Value.ToString(CultureInfo.InvariantCulture)}), 4979)";
139-
var toBox = keepProjection ?
140-
$"st_setsrid(ST_MakePoint({toX}, {toY}, {to.Z.Value.ToString(CultureInfo.InvariantCulture)}), {source_epsg})" :
141-
$"st_setsrid(ST_MakePoint({toX}, {toY}, {to.Z.Value.ToString(CultureInfo.InvariantCulture)}), 4979)";
92+
var fromBox = $"st_setsrid(ST_MakePoint({fromX}, {fromY}, {from.Z.Value.ToString(CultureInfo.InvariantCulture)}), {source_epsg})";
93+
var toBox = $"st_setsrid(ST_MakePoint({toX}, {toY}, {to.Z.Value.ToString(CultureInfo.InvariantCulture)}), {source_epsg})";
14294

14395
var geom = $"st_setsrid(st_makepoint((st_xmin({geometry_column}) + st_xmax({geometry_column}))/2,(st_ymin({geometry_column}) + st_ymax({geometry_column}))/2, (st_zmin({geometry_column}) + st_zmax({geometry_column}))/2), {source_epsg})";
144-
where = keepProjection ?
145-
$"ST_3DIntersects({geom}, ST_3DMakeBox({fromBox}, {toBox})) {query}" :
146-
$"ST_3DIntersects({geom}, st_transform(ST_3DMakeBox({fromBox}, {toBox}), {source_epsg})) {query}";
96+
where = $"ST_3DIntersects({geom}, ST_3DMakeBox({fromBox}, {toBox})) {query}";
14797
}
14898

14999
return where;

‎src/b3dm.tileset/OctreeTiler.cs‎

Lines changed: 4 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,6 @@
11
using System.Collections.Generic;
22
using System.IO;
33
using System.Linq;
4-
using B3dm.Tileset;
54
using Npgsql;
65
using subtree;
76
using Wkx;
@@ -43,7 +42,7 @@ public List<Tile3D> GenerateTiles3D(BoundingBox3D bbox, int level, Tile3D tile,
4342

4443
var where = inputTable.GetQueryClause();
4544

46-
var numberOfFeatures = FeatureCountRepository.CountFeaturesInBox(conn, inputTable.TableName, inputTable.GeometryColumn, new Point(bbox.XMin, bbox.YMin, bbox.ZMin), new Point(bbox.XMax, bbox.YMax, bbox.ZMax), where, inputTable.EPSGCode, tilingSettings.KeepProjection, processedGeometries);
45+
var numberOfFeatures = FeatureCountRepository.CountFeaturesInBox(conn, inputTable.TableName, inputTable.GeometryColumn, new Point(bbox.XMin, bbox.YMin, bbox.ZMin), new Point(bbox.XMax, bbox.YMax, bbox.ZMax), where, inputTable.EPSGCode, processedGeometries);
4746
if (numberOfFeatures == 0) {
4847
var t2 = new Tile3D(level, tile.X, tile.Y, tile.Z);
4948
t2.Available = false;
@@ -76,10 +75,9 @@ public List<Tile3D> GenerateTiles3D(BoundingBox3D bbox, int level, Tile3D tile,
7675
var bbox3d = new BoundingBox3D(xstart, ystart, z_start, xend, yend, zend);
7776

7877
var bboxOctant = new BoundingBox(xstart, ystart, xend, yend);
79-
var filteredProcessedGeometries = GeometryRepository.FilterHashesByEnvelope(conn, inputTable.TableName, inputTable.GeometryColumn, bboxOctant, inputTable.EPSGCode, localProcessedGeometries, tilingSettings.KeepProjection);
8078

8179
var new_tile = new Tile3D(level, tile.X * 2 + x, tile.Y * 2 + y, tile.Z * 2 + z);
82-
GenerateTiles3D(bbox3d, level, new_tile, tiles, tileBounds, filteredProcessedGeometries);
80+
GenerateTiles3D(bbox3d, level, new_tile, tiles, tileBounds, localProcessedGeometries);
8381
}
8482
}
8583
}
@@ -101,7 +99,7 @@ private HashSet<string> CreateTileForLargestGeometries3D(BoundingBox3D bbox, int
10199
int target_srs = tilingSettings.KeepProjection ? inputTable.EPSGCode : 4978;
102100

103101
var bbox1 = new double[] { bbox.XMin, bbox.YMin, bbox.XMax, bbox.YMax, bbox.ZMin, bbox.ZMax };
104-
var geometriesToProcess = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, bbox1, inputTable.EPSGCode, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, tilingSettings.KeepProjection, processedGeometries, tilingSettings.MaxFeaturesPerTile, tilingSettings.SortBy);
102+
var geometriesToProcess = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, bbox1, inputTable.EPSGCode, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, processedGeometries, tilingSettings.MaxFeaturesPerTile, tilingSettings.SortBy);
105103

106104
if (geometriesToProcess.Count > 0) {
107105
foreach (var geom in geometriesToProcess.Where(geom => !string.IsNullOrEmpty(geom.Hash))) {
@@ -127,7 +125,7 @@ private void CreateTile3D(BoundingBox3D bbox, int level, Tile3D tile, List<Tile3
127125
int target_srs = tilingSettings.KeepProjection ? inputTable.EPSGCode : 4978;
128126

129127
var bbox1 = new double[] { bbox.XMin, bbox.YMin, bbox.XMax, bbox.YMax, bbox.ZMin, bbox.ZMax };
130-
var geometries = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, bbox1, inputTable.EPSGCode, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, tilingSettings.KeepProjection, processedGeometries, tilingSettings.MaxFeaturesPerTile, tilingSettings.SortBy);
128+
var geometries = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, bbox1, inputTable.EPSGCode, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, processedGeometries, tilingSettings.MaxFeaturesPerTile, tilingSettings.SortBy);
131129

132130
if (geometries.Count > 0) {
133131
// Collect hashes of processed geometries

‎src/b3dm.tileset/QuadtreeTiler.cs‎

Lines changed: 4 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -55,7 +55,7 @@ public List<Tile> GenerateTiles(BoundingBox bbox, Tile tile, List<Tile> tiles, i
5555
where += $" and {lodquery}";
5656
}
5757

58-
var numberOfFeatures = FeatureCountRepository.CountFeaturesInBox(conn, inputTable.TableName, inputTable.GeometryColumn, new Point(bbox.XMin, bbox.YMin), new Point(bbox.XMax, bbox.YMax), where, source_epsg, keepProjection, processedGeometries);
58+
var numberOfFeatures = FeatureCountRepository.CountFeaturesInBox(conn, inputTable.TableName, inputTable.GeometryColumn, new Point(bbox.XMin, bbox.YMin), new Point(bbox.XMax, bbox.YMax), where, source_epsg, processedGeometries);
5959

6060
if (numberOfFeatures == 0) {
6161
tile.Available = false;
@@ -82,9 +82,7 @@ public List<Tile> GenerateTiles(BoundingBox bbox, Tile tile, List<Tile> tiles, i
8282
var new_tile = new Tile(z, tile.X * 2 + x, tile.Y * 2 + y);
8383
new_tile.BoundingBox = bboxQuad.ToArray();
8484

85-
var filteredProcessedGeometries = GeometryRepository.FilterHashesByEnvelope(conn, inputTable.TableName, inputTable.GeometryColumn, bboxQuad, source_epsg, localProcessedGeometries, keepProjection);
86-
87-
GenerateTiles(bboxQuad, new_tile, tiles, lod, createGltf, keepProjection, filteredProcessedGeometries);
85+
GenerateTiles(bboxQuad, new_tile, tiles, lod, createGltf, keepProjection, localProcessedGeometries);
8886
}
8987
}
9088
}
@@ -107,7 +105,7 @@ private HashSet<string> CreateTileForLargestGeometries(BoundingBox bbox, Tile ti
107105

108106
int target_srs = keepProjection ? source_epsg : 4978;
109107

110-
var geometriesToProcess = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, tile.BoundingBox, source_epsg, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, keepProjection, processedGeometries, maxFeaturesPerTile, sortBy: sortBy);
108+
var geometriesToProcess = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, tile.BoundingBox, source_epsg, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, processedGeometries, maxFeaturesPerTile, sortBy: sortBy);
111109

112110
if (geometriesToProcess.Count > 0) {
113111

@@ -161,7 +159,7 @@ private void CreateTile(BoundingBox bbox, Tile tile, List<Tile> tiles, string wh
161159

162160
int target_srs = keepProjection ? source_epsg : 4978;
163161

164-
var geometries = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, tile.BoundingBox, source_epsg, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, keepProjection, processedGeometries, sortBy: sortBy);
162+
var geometries = GeometryRepository.GetGeometrySubset(conn, inputTable.TableName, inputTable.GeometryColumn, tile.BoundingBox, source_epsg, target_srs, inputTable.ShadersColumn, inputTable.AttributeColumns, where, inputTable.RadiusColumn, processedGeometries, sortBy: sortBy);
165163

166164
if (geometries.Count > 0) {
167165
// Collect hashes of processed geometries

‎src/pg2b3dm.database.tests/UnitTest1.cs‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -43,9 +43,12 @@ public void TestArvieuxBuildingsOctree()
4343
OutputDirectoryCreator.GetFolders("output");
4444
var connectionString = _containerPostgres.GetConnectionString();
4545
var conn = new NpgsqlConnection(connectionString);
46-
var bbox_table = BoundingBoxRepository.GetBoundingBoxForTable(conn, "arvieux_batiments", "geom");
46+
// Get bbox in source projection (EPSG:5698) for spatial queries
47+
var bbox_table = BoundingBoxRepository.GetBoundingBoxForTable(conn, "arvieux_batiments", "geom", true);
48+
// Get WGS84 bbox separately for ECEF translation
49+
var bbox_wgs84 = BoundingBoxRepository.GetBoundingBoxForTable(conn, "arvieux_batiments", "geom", false);
4750

48-
var center_wgs84 = bbox_table.bbox.GetCenter();
51+
var center_wgs84 = bbox_wgs84.bbox.GetCenter();
4952
var translation = SpatialConverter.GeodeticToEcef((double)center_wgs84.X!, (double)center_wgs84.Y!, 0);
5053
var trans = new double[] { translation.X, translation.Y, };
5154

0 commit comments

Comments
 (0)