diff --git a/noisemodelling-jdbc/src/main/java/org/noise_planet/noisemodelling/jdbc/BezierContouring.java b/noisemodelling-jdbc/src/main/java/org/noise_planet/noisemodelling/jdbc/BezierContouring.java index d1fd24e03..b96ef6ed8 100644 --- a/noisemodelling-jdbc/src/main/java/org/noise_planet/noisemodelling/jdbc/BezierContouring.java +++ b/noisemodelling-jdbc/src/main/java/org/noise_planet/noisemodelling/jdbc/BezierContouring.java @@ -21,6 +21,7 @@ */ package org.noise_planet.noisemodelling.jdbc; +import org.h2gis.utilities.GeometryTableUtilities; import org.h2gis.utilities.JDBCUtilities; import org.h2gis.utilities.TableLocation; import org.h2gis.utilities.jts_utils.Contouring; @@ -53,6 +54,8 @@ public class BezierContouring { List isoLevels; List isoLabels; boolean smooth = true; + + boolean mergeTriangles = true; double smoothCoefficient = 1.0; double deltaPoints = 0.5; // minimal distance between bezier points double epsilon = 0.05; @@ -96,6 +99,20 @@ public double getSmoothCoefficient() { return smoothCoefficient; } + /** + * @return True if triangles will be merged using isolevel attribute + */ + public boolean isMergeTriangles() { + return mergeTriangles; + } + + /** + * @param mergeTriangles True if triangles will be merged using isolevel attribute, Z ordinate will be lost + */ + public void setMergeTriangles(boolean mergeTriangles) { + this.mergeTriangles = mergeTriangles; + } + /** * @param smoothCoefficient Coefficient of polygons smoothing [0-1] */ @@ -181,11 +198,11 @@ static List curve4(Coordinate anchor1, Coordinate control1, Coordina static Coordinate[] generateBezierCurves(Coordinate[] coordinates, Quadtree segmentTree, double pointsDelta) { ArrayList pts = new ArrayList<>(); - pts.add(coordinates[0]); + pts.add(new Coordinate(coordinates[0].x, coordinates[0].y)); for(int i = 0; i < coordinates.length - 1; i++) { final int i2 = i + 1; - Coordinate p1 = coordinates[i]; - Coordinate p2 = coordinates[i2]; + Coordinate p1 = new Coordinate(coordinates[i].x, coordinates[i].y); + Coordinate p2 = new Coordinate(coordinates[i2].x, coordinates[i2].y); Segment segment = new Segment(p1, p2); List segments = (List)segmentTree.query(segment.getEnvelope()); @@ -409,7 +426,7 @@ void processCell(Connection connection, int cellId, Map> entry : polys.entrySet()) { ArrayList polygons = new ArrayList<>(); - if(!smooth) { + if(!smooth && mergeTriangles) { // Merge triangles try { CascadedPolygonUnion union = new CascadedPolygonUnion(entry.getValue()); @@ -454,12 +471,21 @@ public void createTable(Connection connection) throws SQLException { Map> polyMap = new HashMap<>(); int lastCellId = -1; try(Statement st = connection.createStatement()) { + String geometryType = "GEOMETRY(POLYGONZ,"+srid+")"; + if(smooth || mergeTriangles) { + geometryType = "GEOMETRY(POLYGON,"+srid+")"; + } st.execute("DROP TABLE IF EXISTS " + TableLocation.parse(outputTable)); - st.execute("CREATE TABLE " + TableLocation.parse(outputTable) + "(PK SERIAL, CELL_ID INTEGER, THE_GEOM GEOMETRY, ISOLVL INTEGER, ISOLABEL VARCHAR);"); - String query = "SELECT CELL_ID, ST_X(p1.the_geom) xa,ST_Y(p1.the_geom) ya,ST_X(p2.the_geom) xb,ST_Y(p2.the_geom) yb,ST_X(p3.the_geom) xc,ST_Y(p3.the_geom) yc, p1."+pointTableField+" lvla, p2."+pointTableField+" lvlb, p3."+pointTableField+" lvlc FROM "+triangleTable+" t, "+pointTable+" p1,"+pointTable+" p2,"+pointTable+" p3 WHERE t.PK_1 = p1."+pkField+" and t.PK_2 = p2."+pkField+" AND t.PK_3 = p3."+pkField+" order by cell_id;"; + st.execute("CREATE TABLE " + TableLocation.parse(outputTable) + + "(PK SERIAL, CELL_ID INTEGER, THE_GEOM "+geometryType+", ISOLVL INTEGER, ISOLABEL VARCHAR);"); + String query = "SELECT CELL_ID, ST_X(p1.the_geom) xa,ST_Y(p1.the_geom) ya, ST_Z(p1.the_geom) za," + + "ST_X(p2.the_geom) xb,ST_Y(p2.the_geom) yb, ST_Z(p2.the_geom) zb," + + "ST_X(p3.the_geom) xc,ST_Y(p3.the_geom) yc, ST_Z(p3.the_geom) zc," + + " p1."+pointTableField+" lvla, p2."+pointTableField+" lvlb, p3."+pointTableField+" lvlc FROM "+triangleTable+" t, "+pointTable+" p1,"+pointTable+" p2,"+pointTable+" p3 WHERE t.PK_1 = p1."+pkField+" and t.PK_2 = p2."+pkField+" AND t.PK_3 = p3."+pkField+" order by cell_id;"; try(ResultSet rs = st.executeQuery(query)) { // Cache columns index - int xa = 0, xb = 0, xc = 0, ya = 0, yb = 0, yc = 0, lvla = 0, lvlb = 0, lvlc = 0, cell_id = 0; + int xa = 0, xb = 0, xc = 0, ya = 0, yb = 0, yc = 0, za = 0, zb = 1, zc = 1, lvla = 0, lvlb = 0, + lvlc = 0, cell_id = 0; ResultSetMetaData resultSetMetaData = rs.getMetaData(); for (int columnId = 1; columnId <= resultSetMetaData.getColumnCount(); columnId++) { switch (resultSetMetaData.getColumnLabel(columnId).toUpperCase()) { @@ -481,6 +507,15 @@ public void createTable(Connection connection) throws SQLException { case "YC": yc = columnId; break; + case "ZA": + za = columnId; + break; + case "ZB": + zb = columnId; + break; + case "ZC": + zc = columnId; + break; case "LVLA": lvla = columnId; break; @@ -495,8 +530,8 @@ public void createTable(Connection connection) throws SQLException { break; } } - if (xa == 0 || xb == 0 || xc == 0 || ya == 0 || yb == 0 || yc == 0 || lvla == 0 || lvlb == 0 || - lvlc == 0 || cell_id == 0) { + if (xa == 0 || xb == 0 || xc == 0 || ya == 0 || yb == 0 || yc == 0 || za == 0 || zb == 0 || zc == 0 + || lvla == 0 || lvlb == 0 || lvlc == 0 || cell_id == 0) { throw new SQLException("Missing field in input tables"); } while(rs.next()) { @@ -508,9 +543,9 @@ public void createTable(Connection connection) throws SQLException { } lastCellId = cellId; // Split current triangle - Coordinate a = new Coordinate(rs.getDouble(xa), rs.getDouble(ya)); - Coordinate b = new Coordinate(rs.getDouble(xb), rs.getDouble(yb)); - Coordinate c = new Coordinate(rs.getDouble(xc), rs.getDouble(yc)); + Coordinate a = new Coordinate(rs.getDouble(xa), rs.getDouble(ya), rs.getDouble(za)); + Coordinate b = new Coordinate(rs.getDouble(xb), rs.getDouble(yb), rs.getDouble(zb)); + Coordinate c = new Coordinate(rs.getDouble(xc), rs.getDouble(yc), rs.getDouble(zc)); // Fetch data TriMarkers triMarkers = new TriMarkers(a, b, c, dbaToW(rs.getDouble(lvla)), dbaToW(rs.getDouble(lvlb)), diff --git a/noisemodelling-jdbc/src/main/java/org/noise_planet/noisemodelling/jdbc/TriangleNoiseMap.java b/noisemodelling-jdbc/src/main/java/org/noise_planet/noisemodelling/jdbc/TriangleNoiseMap.java index 3d814f99f..9e4817ad1 100644 --- a/noisemodelling-jdbc/src/main/java/org/noise_planet/noisemodelling/jdbc/TriangleNoiseMap.java +++ b/noisemodelling-jdbc/src/main/java/org/noise_planet/noisemodelling/jdbc/TriangleNoiseMap.java @@ -299,6 +299,60 @@ public void setEpsilon(double epsilon) { this.epsilon = epsilon; } + public static void generateResultTable(Connection connection, String receiverTableName, String trianglesTableName, + AtomicInteger receiverPK, List vertices, + GeometryFactory geometryFactory, List triangles, int cellI, + int cellJ, int gridDim) throws SQLException { + + if(!JDBCUtilities.tableExists(connection, receiverTableName)) { + Statement st = connection.createStatement(); + st.execute("CREATE TABLE "+TableLocation.parse(receiverTableName)+"(pk serial NOT NULL, the_geom geometry not null, PRIMARY KEY (PK))"); + } + if(!JDBCUtilities.tableExists(connection, trianglesTableName)) { + Statement st = connection.createStatement(); + st.execute("CREATE TABLE "+TableLocation.parse(trianglesTableName)+"(pk serial NOT NULL, the_geom geometry , PK_1 integer not null, PK_2 integer not null, PK_3 integer not null, cell_id integer not null, PRIMARY KEY (PK))"); + } + int receiverPkOffset = receiverPK.get(); + // Add vertices to receivers + PreparedStatement ps = connection.prepareStatement("INSERT INTO "+TableLocation.parse(receiverTableName)+" VALUES (?, ?);"); + int batchSize = 0; + for(Coordinate v : vertices) { + ps.setInt(1, receiverPK.getAndAdd(1)); + ps.setObject(2, geometryFactory.createPoint(v)); + ps.addBatch(); + batchSize++; + if (batchSize >= BATCH_MAX_SIZE) { + ps.executeBatch(); + ps.clearBatch(); + batchSize = 0; + } + } + if (batchSize > 0) { + ps.executeBatch(); + } + // Add triangles + ps = connection.prepareStatement("INSERT INTO "+TableLocation.parse(trianglesTableName)+"(the_geom, PK_1, PK_2, PK_3, CELL_ID) VALUES (?, ?, ?, ?, ?);"); + batchSize = 0; + for(Triangle t : triangles) { + ps.setObject(1, geometryFactory.createPolygon(new Coordinate[]{vertices.get(t.getA()), + vertices.get(t.getB()), vertices.get(t.getC()), vertices.get(t.getA())})); + ps.setInt(2, t.getA() + receiverPkOffset); + ps.setInt(3, t.getC() + receiverPkOffset); + ps.setInt(4, t.getB() + receiverPkOffset); + ps.setInt(5, cellI * gridDim + cellJ); + ps.addBatch(); + batchSize++; + if (batchSize >= BATCH_MAX_SIZE) { + ps.executeBatch(); + ps.clearBatch(); + batchSize = 0; + } + } + if (batchSize > 0) { + ps.executeBatch(); + } + } + public void generateReceivers(Connection connection, int cellI, int cellJ, String receiverTableName, String trianglesTableName, AtomicInteger receiverPK) throws SQLException, LayerDelaunayError, IOException { int ij = cellI * gridDim + cellJ + 1; @@ -367,53 +421,8 @@ public void generateReceivers(Connection connection, int cellI, int cellJ, Strin } nbreceivers += vertices.size(); - if(!JDBCUtilities.tableExists(connection, receiverTableName)) { - Statement st = connection.createStatement(); - st.execute("CREATE TABLE "+TableLocation.parse(receiverTableName)+"(pk serial NOT NULL, the_geom geometry not null, PRIMARY KEY (PK))"); - } - if(!JDBCUtilities.tableExists(connection, trianglesTableName)) { - Statement st = connection.createStatement(); - st.execute("CREATE TABLE "+TableLocation.parse(trianglesTableName)+"(pk serial NOT NULL, the_geom geometry , PK_1 integer not null, PK_2 integer not null, PK_3 integer not null, cell_id integer not null, PRIMARY KEY (PK))"); - } - int receiverPkOffset = receiverPK.get(); - // Add vertices to receivers - PreparedStatement ps = connection.prepareStatement("INSERT INTO "+TableLocation.parse(receiverTableName)+" VALUES (?, ?);"); - int batchSize = 0; - for(Coordinate v : vertices) { - ps.setInt(1, receiverPK.getAndAdd(1)); - ps.setObject(2, geometryFactory.createPoint(v)); - ps.addBatch(); - batchSize++; - if (batchSize >= BATCH_MAX_SIZE) { - ps.executeBatch(); - ps.clearBatch(); - batchSize = 0; - } - } - if (batchSize > 0) { - ps.executeBatch(); - } - // Add triangles - ps = connection.prepareStatement("INSERT INTO "+TableLocation.parse(trianglesTableName)+"(the_geom, PK_1, PK_2, PK_3, CELL_ID) VALUES (?, ?, ?, ?, ?);"); - batchSize = 0; - for(Triangle t : triangles) { - ps.setObject(1, geometryFactory.createPolygon(new Coordinate[]{vertices.get(t.getA()), - vertices.get(t.getB()), vertices.get(t.getC()), vertices.get(t.getA())})); - ps.setInt(2, t.getA() + receiverPkOffset); - ps.setInt(3, t.getC() + receiverPkOffset); - ps.setInt(4, t.getB() + receiverPkOffset); - ps.setInt(5, cellI * gridDim + cellJ); - ps.addBatch(); - batchSize++; - if (batchSize >= BATCH_MAX_SIZE) { - ps.executeBatch(); - ps.clearBatch(); - batchSize = 0; - } - } - if (batchSize > 0) { - ps.executeBatch(); - } + generateResultTable(connection, receiverTableName, trianglesTableName, receiverPK, vertices, geometryFactory, + triangles, cellI, cellJ, gridDim); } public double getRoadWidth() { diff --git a/noisemodelling-jdbc/src/test/java/org/noise_planet/noisemodelling/jdbc/BezierContouringJDBCTest.java b/noisemodelling-jdbc/src/test/java/org/noise_planet/noisemodelling/jdbc/BezierContouringJDBCTest.java index 5974980c4..1c18b5274 100644 --- a/noisemodelling-jdbc/src/test/java/org/noise_planet/noisemodelling/jdbc/BezierContouringJDBCTest.java +++ b/noisemodelling-jdbc/src/test/java/org/noise_planet/noisemodelling/jdbc/BezierContouringJDBCTest.java @@ -5,19 +5,28 @@ import org.h2gis.functions.io.geojson.GeoJsonRead; import org.h2gis.functions.io.shp.SHPRead; import org.h2gis.functions.io.shp.SHPWrite; +import org.h2gis.functions.spatial.mesh.DelaunayData; import org.h2gis.utilities.JDBCUtilities; +import org.h2gis.utilities.SpatialResultSet; +import org.h2gis.utilities.TableLocation; +import org.h2gis.utilities.TableUtilities; import org.junit.After; import org.junit.Before; import org.junit.Test; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.GeometryFactory; import org.noise_planet.noisemodelling.pathfinder.LayerDelaunayError; +import org.noise_planet.noisemodelling.pathfinder.LayerTinfour; import java.io.IOException; -import java.sql.Connection; -import java.sql.SQLException; -import java.sql.Statement; +import java.nio.file.Paths; +import java.sql.*; +import java.util.Arrays; +import java.util.Collections; import java.util.List; import java.util.concurrent.atomic.AtomicInteger; +import static org.junit.Assert.assertEquals; import static org.junit.Assert.assertTrue; public class BezierContouringJDBCTest { @@ -70,6 +79,52 @@ public void testBezierContouring() throws SQLException, IOException { assertTrue(fieldValues.contains("8")); assertTrue(fieldValues.contains("9")); - SHPWrite.exportTable(connection, "target/contouring.shp", "CONTOURING_NOISE_MAP","UTF-8",true); + } + + @Test + public void testContouring3D() throws SQLException, IOException, LayerDelaunayError { + // Will create elevation iso from DEM table + GeoJsonRead.importTable(connection, Paths.get(Paths.get(System.getProperty("user.dir")).getParent().toString(), + "wps_scripts/src/test/resources/org/noise_planet/noisemodelling/wps/dem.geojson").toString()); + LayerTinfour delaunayTool = new LayerTinfour(); + try (PreparedStatement st = connection.prepareStatement( + "SELECT the_geom FROM DEM")) { + try (SpatialResultSet rs = st.executeQuery().unwrap(SpatialResultSet.class)) { + while (rs.next()) { + Geometry pt = rs.getGeometry(); + if(pt != null) { + delaunayTool.addVertex(pt.getCoordinate()); + } + } + } + } + delaunayTool.processDelaunay(); + TriangleNoiseMap.generateResultTable(connection, "RECEIVERS", "TRIANGLES", + new AtomicInteger(), delaunayTool.getVertices(), new GeometryFactory(), delaunayTool.getTriangles(), + 0, 0, 1); + try(Statement st = connection.createStatement()) { + st.execute("ALTER TABLE RECEIVERS ADD COLUMN HEIGHT FLOAT"); + st.execute("UPDATE RECEIVERS SET HEIGHT = ST_Z(THE_GEOM)"); + } + long start = System.currentTimeMillis(); + BezierContouring bezierContouring = new BezierContouring(Arrays.asList(0.,5.,10.,15.,20.,25.,30.,35.), 2154); + bezierContouring.setPointTable("RECEIVERS"); + bezierContouring.setPointTableField("HEIGHT"); + bezierContouring.setSmooth(false); + bezierContouring.setMergeTriangles(false); + bezierContouring.createTable(connection); + System.out.println("Contouring done in " + (System.currentTimeMillis() - start) + " ms"); + + assertTrue(JDBCUtilities.tableExists(connection, "CONTOURING_NOISE_MAP")); + + // Check Z values in CONTOURING_NOISE_MAP + try(Statement st = connection.createStatement()) { + try(ResultSet rs = st.executeQuery("SELECT MAX(ST_ZMAX(THE_GEOM)) MAXZ, MIN(ST_ZMIN(THE_GEOM)) MINZ FROM CONTOURING_NOISE_MAP")) { + assertTrue(rs.next()); + assertEquals(33.2, rs.getDouble("MAXZ"), 0.01); + assertEquals(-1.79, rs.getDouble("MINZ"), 0.01); + } + } + } } \ No newline at end of file