changeset 664:a3d722e1f593 octree

octree: Load point cloud data projected in a suited UTM zone.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 18 Sep 2018 13:41:04 +0200
parents d856e458dd64
children 7345a249c9e1
files cmd/tin2octree/main.go cmd/tin2octree/tin.go
diffstat 2 files changed, 68 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
--- a/cmd/tin2octree/main.go	Tue Sep 18 00:13:51 2018 +0200
+++ b/cmd/tin2octree/main.go	Tue Sep 18 13:41:04 2018 +0200
@@ -6,6 +6,7 @@
 	"flag"
 	"io"
 	"log"
+	"math"
 	"os"
 	"time"
 
@@ -19,7 +20,9 @@
 	date       = flag.String("date", "", "date info")
 	file       = flag.String("file", "", "save to file")
 
-	snap       = flag.Bool("snappy", false, "use snappy compression")
+	utm  = flag.Bool("utm", false, "fetch in matchin UTM zone")
+	snap = flag.Bool("snappy", false, "use snappy compression")
+
 	dbhost     = flag.String("dbhost", "localhost", "database host")
 	dbport     = flag.Uint("dbport", 5432, "database port")
 	dbname     = flag.String("dbname", "gemma", "database user")
@@ -49,11 +52,36 @@
 	return fn(db)
 }
 
-const tinSQL = `
+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
+`
+	centroidSQL = `
+SELECT ST_X(ST_Centroid(point_cloud::geometry)), ST_Y(ST_Centroid(point_cloud::geometry))
+FROM waterway.sounding_results
+WHERE bottleneck_id = $1 AND date_info = $2
+`
+)
+
+func utmZone(x, y float64) int {
+	var pref int
+	if y > 0 {
+		pref = 32600
+	} else {
+		pref = 32700
+	}
+	zone := int(math.Floor((x+180)/6)) + 1
+	return zone + pref
+}
 
 func main() {
 	flag.Parse()
@@ -70,8 +98,29 @@
 	var t tin
 
 	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)
+		}
+
 		start := time.Now()
-		err := db.QueryRow(tinSQL, *bottleneck, dateInfo).Scan(&t)
+		var err error
+		if *utm {
+			err = db.QueryRow(tinUTMSQL, *bottleneck, dateInfo, utmZ).Scan(&t)
+		} else {
+			err = db.QueryRow(tinSQL, *bottleneck, dateInfo).Scan(&t)
+		}
 		switch {
 		case err == sql.ErrNoRows:
 			return nil
@@ -79,6 +128,11 @@
 			return err
 		}
 		log.Printf("query took: %s\n", time.Since(start))
+
+		if *utm {
+			t.epsg = uint32(utmZ)
+		}
+
 		return nil
 	}); err != nil {
 		log.Fatalf("error: %v\n", err)
--- a/cmd/tin2octree/tin.go	Tue Sep 18 00:13:51 2018 +0200
+++ b/cmd/tin2octree/tin.go	Tue Sep 18 13:41:04 2018 +0200
@@ -22,7 +22,10 @@
 	wkbTriangleZ uint32 = 1000 + 17
 )
 
+const wgs84 = 4326
+
 type tin struct {
+	epsg      uint32
 	vertices  []vertex
 	triangles [][]int32
 
@@ -155,7 +158,11 @@
 		triangles = append(triangles, triangle)
 	}
 
+	log.Printf("bbox: [[%f, %f], [%f, %f]]\n",
+		min.x, min.y, max.x, max.y)
+
 	*t = tin{
+		epsg:      wgs84,
 		vertices:  vertices,
 		triangles: triangles,
 		min:       min,
@@ -186,6 +193,10 @@
 
 func (t *tin) Serialize(w io.Writer) error {
 
+	if err := binary.Write(w, binary.LittleEndian, t.epsg); err != nil {
+		return err
+	}
+
 	if err := t.min.write(w); err != nil {
 		return err
 	}