Mercurial > gemma

changeset 969:a4e003ba0074

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Throw out triangles out of triangulation result which are not covered by the the bounding polygon of the sounding result.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 17 Oct 2018 17:23:46 +0200
parents a4fe07a21ba7
children 166782ebab96
files cmd/tin2octree/main.go
diffstat 1 files changed, 28 insertions(+), 33 deletions(-) [+]
line wrap: on
line diff
--- a/cmd/tin2octree/main.go	Wed Oct 17 16:00:49 2018 +0200
+++ b/cmd/tin2octree/main.go	Wed Oct 17 17:23:46 2018 +0200
@@ -23,7 +23,6 @@
 	date       = flag.String("date", "", "date info")
 	file       = flag.String("file", "", "save to file")
 	insert     = flag.Bool("insert", false, "write as SQL insert statement")
-	utm        = flag.Bool("utm", true, "fetch in matchin UTM zone")
 
 	dbhost     = flag.String("dbhost", "localhost", "database host")
 	dbport     = flag.Uint("dbport", 5432, "database port")
@@ -55,18 +54,23 @@
 }
 
 const (
-	tinSQL = `
-SELECT ST_AsBinary(ST_DelaunayTriangles(point_cloud::geometry, 0, 2))
-FROM waterway.sounding_results
-WHERE bottleneck_id = $1 AND date_info = $2
-`
 	tinUTMSQL = `
-SELECT ST_AsBinary(
-  ST_DelaunayTriangles(
-	  ST_Transform(point_cloud::geometry, $3::int), 0, 2))
-FROM waterway.sounding_results
-WHERE bottleneck_id = $1 AND date_info = $2
-`
+WITH trans AS (
+  SELECT
+	ST_Buffer(ST_Transform(area::geometry, $3::int), 0.001) AS area,
+	ST_Transform(point_cloud::geometry, $3::int) AS point_cloud
+  FROM waterway.sounding_results
+  WHERE bottleneck_id = $1 AND date_info = $2
+),
+triangles AS (
+  SELECT t.geom AS geom, ST_MakePolygon(ST_ExteriorRing(t.geom)) AS poly FROM (
+    SELECT (st_dump(
+      ST_DelaunayTriangles(point_cloud, 0, 2))).geom
+    FROM trans) t
+)
+SELECT ST_AsBinary(ST_Collect(triangles.geom)) FROM triangles, trans
+WHERE st_covers(trans.area, triangles.poly)`
+
 	centroidSQL = `
 SELECT ST_X(ST_Centroid(point_cloud::geometry)), ST_Y(ST_Centroid(point_cloud::geometry))
 FROM waterway.sounding_results
@@ -102,27 +106,20 @@
 	if err := run(func(db *sql.DB) error {
 		var utmZ int
 
-		if *utm {
-			var cx, cy float64
-			err := db.QueryRow(centroidSQL, *bottleneck, dateInfo).Scan(&cx, &cy)
-			switch {
-			case err == sql.ErrNoRows:
-				return nil
-			case err != nil:
-				return err
-			}
-			log.Printf("lat/lon: [%f, %f]\n", cx, cy)
-			utmZ = utmZone(cx, cy)
-			log.Printf("UTM zone: %d\n", utmZ)
+		var cx, cy float64
+		err := db.QueryRow(centroidSQL, *bottleneck, dateInfo).Scan(&cx, &cy)
+		switch {
+		case err == sql.ErrNoRows:
+			return nil
+		case err != nil:
+			return err
 		}
+		log.Printf("lat/lon: [%f, %f]\n", cx, cy)
+		utmZ = utmZone(cx, cy)
+		log.Printf("UTM zone: %d\n", utmZ)
 
 		start := time.Now()
-		var err error
-		if *utm {
-			err = db.QueryRow(tinUTMSQL, *bottleneck, dateInfo, utmZ).Scan(&t)
-		} else {
-			err = db.QueryRow(tinSQL, *bottleneck, dateInfo).Scan(&t)
-		}
+		err = db.QueryRow(tinUTMSQL, *bottleneck, dateInfo, utmZ).Scan(&t)
 		switch {
 		case err == sql.ErrNoRows:
 			return nil
@@ -131,9 +128,7 @@
 		}
 		log.Printf("query took: %s\n", time.Since(start))
 
-		if *utm {
-			t.EPSG = uint32(utmZ)
-		}
+		t.EPSG = uint32(utmZ)
 
 		return nil
 	}); err != nil {