Mercurial > gemma

changeset 2529:45d51a49f191

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
SR import: Use own triangulation and clipping when importing sounding results.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 06 Mar 2019 17:51:58 +0100
parents 113912129481
children 3f61b84ae7a6
files pkg/imports/sr.go pkg/octree/contours.go pkg/octree/tin.go pkg/octree/vertex.go schema/gemma.sql
diffstat 5 files changed, 183 insertions(+), 90 deletions(-) [+]
line wrap: on
line diff
--- a/pkg/imports/sr.go	Wed Mar 06 16:26:45 2019 +0100
+++ b/pkg/imports/sr.go	Wed Mar 06 17:51:58 2019 +0100
@@ -104,36 +104,33 @@
   WHERE import_id = $1 AND
         relation = 'waterway.sounding_results'::regclass)`
 
-	insertPointsSQL = `
+	insertHullSQL = `
 INSERT INTO waterway.sounding_results (
   bottleneck_id,
   date_info,
   depth_reference,
-  point_cloud,
   area
-) VALUES (
+) SELECT
   (SELECT id from waterway.bottlenecks where objnam = $1),
   $2::date,
   $3,
-  ST_Transform(ST_GeomFromWKB($4, $6::integer), 4326)::geography,
-  (SELECT CASE WHEN $5::bytea IS NULL THEN
-    ST_Transform(
-      ST_ConcaveHull(
-          ST_Force2D(ST_GeomFromWKB($4, $6::integer)), 0.7), 4326)::geography
-  ELSE
-    ST_Transform(ST_GeomFromWKB($5, $6::integer), 4326)::geography
-  END)
-)
+  (SELECT
+    CASE WHEN $5::bytea IS NULL THEN
+      ST_Transform(ST_ConcaveHull(ST_Force2D(ST_GeomFromWKB($4, $6::integer)), 0.7), 4326)::geography
+    ELSE
+      ST_Transform(ST_GeomFromWKB($5, $6::integer), 4326)::geography
+    END)
 RETURNING
   id,
-  ST_X(ST_Centroid(point_cloud::geometry)),
-  ST_Y(ST_Centroid(point_cloud::geometry)),
-  CASE WHEN ST_Y(ST_Centroid(point_cloud::geometry)) > 0 THEN
-    32600
-  ELSE
-    32700
-  END + floor((ST_X(ST_Centroid(point_cloud::geometry))+180)/6)::int + 1`
+  ST_X(ST_Centroid(area::geometry)),
+  ST_Y(ST_Centroid(area::geometry)),
+  best_utm(area::geometry),
+  ST_AsBinary(ST_Transform(area::geometry, best_utm(area::geometry)))
+`
 
+	reprojectPointsSQL = `
+SELECT ST_AsBinary(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))
+`
 	insertOctreeSQL = `
 UPDATE waterway.sounding_results SET
   octree_checksum = $2, octree_index = $3
@@ -174,6 +171,8 @@
 	feedback Feedback,
 ) (interface{}, error) {
 
+	begin := time.Now()
+
 	z, err := zip.OpenReader(filepath.Join(sr.Dir, "sr.zip"))
 	if err != nil {
 		return nil, err
@@ -234,63 +233,116 @@
 	)
 	start := time.Now()
 
-	err = tx.QueryRow(insertPointsSQL,
+	var hull []byte
+
+	xyzWKB := xyz.AsWKB()
+
+	err = tx.QueryRow(insertHullSQL,
 		m.Bottleneck,
 		m.Date.Time,
 		m.DepthReference,
-		xyz.AsWKB(),
+		xyzWKB,
 		polygon.asWKB(),
 		m.EPSG,
-	).Scan(&id, &lat, &lon, &epsg)
+	).Scan(
+		&id,
+		&lat,
+		&lon,
+		&epsg,
+		&hull,
+	)
 	xyz, polygon = nil, nil // not need from now on.
-	feedback.Info("storing points took %s", time.Since(start))
+	feedback.Info("Calculating hull took %s.", time.Since(start))
+	if err != nil {
+		return nil, err
+	}
+	feedback.Info("Best suited UTM EPSG: %d", epsg)
+
+	start = time.Now()
+
+	var clippingPolygon octree.Polygon
+	if err := clippingPolygon.FromWKB(hull); err != nil {
+		return nil, err
+	}
+	clippingPolygon.Indexify()
+
+	feedback.Info("Building clipping polygon took %v.", time.Since(start))
+
+	start = time.Now()
+
+	var reproj []byte
+
+	err = tx.QueryRow(reprojectPointsSQL,
+		xyzWKB,
+		m.EPSG,
+		epsg,
+	).Scan(&reproj)
+
 	if err != nil {
 		return nil, err
 	}
 
-	feedback.Info("Best suited UTM EPSG: %d", epsg)
+	if err := xyz.FromWKB(reproj); err != nil {
+		return nil, err
+	}
 
-	feedback.Info("Triangulate...")
+	feedback.Info("Reprojecting points took %v.", time.Since(start))
+	feedback.Info("Number of reprojected points: %d", len(xyz))
+
 	start = time.Now()
-	tin, err := octree.GenerateTinByID(ctx, conn, id, epsg)
-	feedback.Info("triangulation took %s", time.Since(start))
+
+	tri, err := octree.Triangulate(xyz)
 	if err != nil {
 		return nil, err
 	}
 
-	if tin == nil {
-		return nil, errors.New("cannot load TIN from database")
-	}
+	feedback.Info("Triangulation took %v.", time.Since(start))
+
+	start = time.Now()
+
+	tin := tri.Tin()
+
+	var str octree.STRTree
+	str.Build(tin)
 
-	feedback.Info("Building octree...")
+	feedback.Info("Building STR tree took %v", time.Since(start))
+
+	start = time.Now()
+
+	removed := str.Clip(&clippingPolygon)
+	feedback.Info("Clipping STR tree took %v.", time.Since(start))
+	feedback.Info("Number of triangles to clip %d.", len(removed))
+
+	start = time.Now()
+
 	builder := octree.NewBuilder(tin)
-	start = time.Now()
-	builder.Build(nil)
+	builder.Build(removed)
 	octreeIndex, err := builder.Bytes()
-	tin = nil // not needed from now on
-	feedback.Info("building octree took %s", time.Since(start))
+	feedback.Info("Building octree took %v.",
+		time.Since(start))
 	if err != nil {
 		return nil, err
 	}
 
-	feedback.Info("Store octree...")
 	start = time.Now()
 	h := sha1.New()
 	h.Write(octreeIndex)
 	checksum := hex.EncodeToString(h.Sum(nil))
 	_, err = tx.Exec(insertOctreeSQL, id, checksum, octreeIndex)
-	feedback.Info("storing octree index took %s", time.Since(start))
+	feedback.Info("Storing octree index took %s.",
+		time.Since(start))
 	if err != nil {
 		return nil, err
 	}
 
 	tree := builder.Tree()
+	tree.EPSG = epsg
 	builder = nil // not needed from now on
 
-	feedback.Info("Generate contour lines...")
 	start = time.Now()
 	err = generateContours(tree, tx, id)
-	feedback.Info("generating and storing contour lines took %s", time.Since(start))
+	feedback.Info("Generating and storing contour lines took %s.",
+		time.Since(start))
 	if err != nil {
 		return nil, err
 	}
@@ -300,10 +352,17 @@
 		return nil, err
 	}
 
-	if err = tx.Commit(); err == nil {
-		feedback.Info("Storing sounding result was successful.")
+	if err = tx.Commit(); err != nil {
+		feedback.Error(
+			"Storing sounding result failed after %v.",
+			time.Since(begin))
+		return nil, err
 	}
 
+	feedback.Info(
+		"Storing sounding result was successful after %v.",
+		time.Since(begin))
+
 	summary := struct {
 		Bottleneck string      `json:"bottleneck"`
 		Date       models.Date `json:"date"`
@@ -473,7 +532,7 @@
 	}
 
 	octree.DoContours(tree, heights, func(res *octree.ContourResult) {
-		if err == nil {
+		if err == nil && len(res.Lines) > 0 {
 			_, err = stmt.Exec(
 				id, res.Height, tree.EPSG,
 				res.Lines.AsWKB2D(),
--- a/pkg/octree/contours.go	Wed Mar 06 16:26:45 2019 +0100
+++ b/pkg/octree/contours.go	Wed Mar 06 17:51:58 2019 +0100
@@ -15,7 +15,6 @@
 package octree
 
 import (
-	"log"
 	"runtime"
 	"sync"
 )
@@ -71,8 +70,6 @@
 		return
 	}
 
-	log.Printf("num heights: %d\n", len(heights))
-
 	contours := make([]*ContourResult, len(heights))
 
 	for i, h := range heights {
--- a/pkg/octree/tin.go	Wed Mar 06 16:26:45 2019 +0100
+++ b/pkg/octree/tin.go	Wed Mar 06 17:51:58 2019 +0100
@@ -15,8 +15,6 @@
 
 import (
 	"bytes"
-	"context"
-	"database/sql"
 	"encoding/binary"
 	"errors"
 	"fmt"
@@ -190,47 +188,6 @@
 	return nil
 }
 
-const (
-	tinSQLPrefix = `WITH trans AS (
-  SELECT
-	ST_Buffer(ST_Transform(area::geometry, $1::int), 0.001) AS area,
-	ST_Transform(point_cloud::geometry, $1::int) AS point_cloud
-  FROM waterway.sounding_results
-`
-	tinSQLSuffix = `
-),
-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)`
-
-	loadTinByIDSQL = tinSQLPrefix + `WHERE id = $2` + tinSQLSuffix
-)
-
-// GenerateTinByID generated a TIN by triangulating a point cloud
-// from the database.
-func GenerateTinByID(
-	ctx context.Context,
-	conn *sql.Conn,
-	id int64,
-	epsg uint32,
-) (*Tin, error) {
-	var tin Tin
-	err := conn.QueryRowContext(ctx, loadTinByIDSQL, epsg, id).Scan(&tin)
-	switch {
-	case err == sql.ErrNoRows:
-		return nil, nil
-	case err != nil:
-		return nil, err
-	}
-	tin.EPSG = epsg
-	return &tin, nil
-}
-
 // Scan implements the sql.Scanner interface.
 func (t *Tin) Scan(raw interface{}) error {
 	if raw == nil {
--- a/pkg/octree/vertex.go	Wed Mar 06 16:26:45 2019 +0100
+++ b/pkg/octree/vertex.go	Wed Mar 06 17:51:58 2019 +0100
@@ -16,6 +16,7 @@
 import (
 	"bytes"
 	"encoding/binary"
+	"fmt"
 	"io"
 	"log"
 	"math"
@@ -1059,3 +1060,84 @@
 
 	return buf.Bytes()
 }
+
+func (mpz *MultiPointZ) FromWKB(data []byte) error {
+
+	r := bytes.NewReader(data)
+
+	endian, err := r.ReadByte()
+
+	var order binary.ByteOrder
+
+	switch {
+	case err != nil:
+		return err
+	case endian == wkb.NDR:
+		order = binary.LittleEndian
+	case endian == wkb.XDR:
+		order = binary.BigEndian
+	default:
+		return fmt.Errorf("unknown byte order %x", endian)
+	}
+
+	var geomType uint32
+	err = binary.Read(r, order, &geomType)
+
+	switch {
+	case err != nil:
+		return err
+	case geomType != wkb.MultiPointZ:
+		return fmt.Errorf("unknown geometry type %x", geomType)
+	}
+
+	var numPoints uint32
+	if err := binary.Read(r, order, &numPoints); err != nil {
+		return err
+	}
+
+	points := make(MultiPointZ, numPoints)
+
+	for i := range points {
+		endian, err = r.ReadByte()
+
+		switch {
+		case err != nil:
+			return err
+		case endian == wkb.NDR:
+			order = binary.LittleEndian
+		case endian == wkb.XDR:
+			order = binary.BigEndian
+		default:
+			return fmt.Errorf("unknown byte order %x", endian)
+		}
+
+		err = binary.Read(r, order, &geomType)
+
+		switch {
+		case err != nil:
+			return err
+		case geomType != wkb.PointZ:
+			return fmt.Errorf("unknown geometry type %x", geomType)
+		}
+
+		var x, y, z uint64
+		if err = binary.Read(r, order, &x); err != nil {
+			return err
+		}
+		if err = binary.Read(r, order, &y); err != nil {
+			return err
+		}
+		if err = binary.Read(r, order, &z); err != nil {
+			return err
+		}
+		points[i] = Vertex{
+			X: math.Float64frombits(x),
+			Y: math.Float64frombits(y),
+			Z: math.Float64frombits(z),
+		}
+	}
+
+	*mpz = points
+
+	return nil
+}
--- a/schema/gemma.sql	Wed Mar 06 16:26:45 2019 +0100
+++ b/schema/gemma.sql	Wed Mar 06 17:51:58 2019 +0100
@@ -505,8 +505,6 @@
         surtyp varchar REFERENCES survey_types,
         coverage varchar REFERENCES coverage_types,
         depth_reference varchar(4) NOT NULL, -- REFERENCES depth_references,
-        point_cloud geography(MULTIPOINTZ, 4326) NOT NULL
-            CHECK(ST_IsSimple(CAST(point_cloud AS geometry))),
         octree_checksum varchar,
         octree_index bytea,
         staging_done boolean NOT NULL DEFAULT false