changeset 3650:01ce3ba9b0d0 single-beam

Fixed generating of random points.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 13 Jun 2019 13:14:40 +0200
parents 810b28f59b8b
children c368a9a20478
files pkg/common/random.go pkg/imports/sr.go pkg/octree/polygon.go pkg/octree/vertex.go
diffstat 4 files changed, 183 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
--- a/pkg/common/random.go	Wed Jun 12 17:30:20 2019 +0200
+++ b/pkg/common/random.go	Thu Jun 13 13:14:40 2019 +0200
@@ -18,6 +18,7 @@
 	"crypto/rand"
 	"io"
 	"log"
+	"math"
 	"math/big"
 	mrand "math/rand"
 	"time"
@@ -71,7 +72,15 @@
 	if low > high {
 		low, high = high, low
 	}
-	rnd := mrand.New(mrand.NewSource(time.Now().Unix()))
+
+	var seed int64
+	if seedInt, err := rand.Int(rand.Reader, big.NewInt(math.MaxInt64)); err != nil {
+		log.Printf("warn: Generating good random seed failed: %v\n", err)
+		seed = time.Now().Unix()
+	} else {
+		seed = seedInt.Int64()
+	}
+	rnd := mrand.New(mrand.NewSource(seed))
 	m := high - low
 	return func() float64 { return rnd.Float64()*m + low }
 }
--- a/pkg/imports/sr.go	Wed Jun 12 17:30:20 2019 +0200
+++ b/pkg/imports/sr.go	Thu Jun 13 13:14:40 2019 +0200
@@ -72,7 +72,7 @@
 
 const (
 	tooLongEdge          = 50.0
-	pointsPerSquareMeter = 5
+	pointsPerSquareMeter = 2
 )
 
 // SRJobKind is the unique name of a SoundingResult import job.
@@ -128,7 +128,7 @@
     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
+      ST_Transform(ST_MakeValid(ST_GeomFromWKB($5, $6::integer)), 4326)::geography
     END)
 RETURNING
   id,
@@ -469,6 +469,109 @@
 
 		xyz = octree.MultiPointZ(generated)
 
+		start = time.Now()
+		tri, err := octree.Triangulate(xyz)
+		if err != nil {
+			return nil, err
+		}
+		log.Printf("Second triangulation took %v.", time.Since(start))
+		log.Printf("Number triangles: %d.", len(tri.Triangles)/3)
+
+		secondTin := tri.Tin()
+
+		//polygon.Buffer(1.0)
+
+		clippingPolygon.FromLineStringZ(polygon)
+		clippingPolygon.Indexify()
+
+		var str octree.STRTree
+		str.Build(secondTin)
+		log.Printf("Building STR tree took %v", time.Since(start))
+
+		start = time.Now()
+
+		removed = str.Clip(&clippingPolygon)
+		log.Printf("Clipping STR tree took %v.", time.Since(start))
+		log.Printf("Number of triangles to clip %d.", len(removed))
+
+		start = time.Now()
+
+		secondTin.EPSG = uint32(epsg)
+
+		builder = octree.NewBuilder(secondTin)
+		builder.Build(removed)
+		octreeIndex, err := builder.Bytes()
+		if err != nil {
+			return nil, err
+		}
+		feedback.Info("Building octree took %v.", time.Since(start))
+
+		var (
+			id       int64
+			epsg     uint32
+			lat, lon float64
+		)
+
+		var hull []byte
+
+		if err := tx.QueryRowContext(
+			ctx,
+			insertHullSQL,
+			m.Bottleneck,
+			m.Date.Time,
+			m.DepthReference,
+			nil,
+			clippingPolygon.AsWKB(),
+			m.EPSG,
+		).Scan(
+			&id,
+			&lat,
+			&lon,
+			&epsg,
+			&hull,
+		); err != nil {
+			return nil, err
+		}
+
+		start = time.Now()
+		h := sha1.New()
+		h.Write(octreeIndex)
+		checksum := hex.EncodeToString(h.Sum(nil))
+		_, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, octreeIndex)
+		if err != nil {
+			return nil, err
+		}
+		feedback.Info("Storing octree index took %s.", time.Since(start))
+
+		tree := builder.Tree()
+
+		start = time.Now()
+		err = generateContours(ctx, tx, tree, id)
+		if err != nil {
+			return nil, err
+		}
+		feedback.Info("Generating and storing contour lines took %s.",
+			time.Since(start))
+
+		// Store for potential later removal.
+		if err = track(ctx, tx, importID, "waterway.sounding_results", id); err != nil {
+			return nil, err
+		}
+
+		summary := struct {
+			Bottleneck string      `json:"bottleneck"`
+			Date       models.Date `json:"date"`
+			Lat        float64     `json:"lat"`
+			Lon        float64     `json:"lon"`
+		}{
+			Bottleneck: m.Bottleneck,
+			Date:       m.Date,
+			Lat:        lat,
+			Lon:        lon,
+		}
+
+		return &summary, nil
+
 	} else {
 		var hull []byte
 		if err = tx.QueryRowContext(
--- a/pkg/octree/polygon.go	Wed Jun 12 17:30:20 2019 +0200
+++ b/pkg/octree/polygon.go	Thu Jun 13 13:14:40 2019 +0200
@@ -409,6 +409,30 @@
 	return raySlope >= diagSlope
 }
 
+func (p *Polygon) AsWKB() []byte {
+
+	size := 1 + 4 + 4
+	for _, r := range p.rings {
+		size += 4 + len(r)*8
+	}
+
+	buf := bytes.NewBuffer(make([]byte, 0, size))
+
+	binary.Write(buf, binary.LittleEndian, wkb.NDR)
+	binary.Write(buf, binary.LittleEndian, wkb.Polygon)
+	binary.Write(buf, binary.LittleEndian, uint32(len(p.rings)))
+
+	for _, r := range p.rings {
+		binary.Write(buf, binary.LittleEndian, uint32(len(r)/2))
+		for i := 0; i < len(r); i += 2 {
+			binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+0]))
+			binary.Write(buf, binary.LittleEndian, math.Float64bits(r[i+1]))
+		}
+	}
+
+	return buf.Bytes()
+}
+
 func (p *Polygon) FromWKB(data []byte) error {
 
 	r := bytes.NewReader(data)
@@ -477,3 +501,14 @@
 
 	return nil
 }
+
+func (p *Polygon) FromLineStringZ(ls LineStringZ) {
+	r := make([]float64, 2*len(ls))
+	var pos int
+	for i := range ls {
+		r[pos+0] = ls[i].X
+		r[pos+1] = ls[i].Y
+		pos += 2
+	}
+	p.rings = []ring{r}
+}
--- a/pkg/octree/vertex.go	Wed Jun 12 17:30:20 2019 +0200
+++ b/pkg/octree/vertex.go	Thu Jun 13 13:14:40 2019 +0200
@@ -210,6 +210,11 @@
 	return math.Hypot(v.X-w.X, v.Y-w.Y)
 }
 
+func (v Vertex) Distance(w Vertex) float64 {
+	v = v.Sub(w)
+	return math.Sqrt(v.Dot(v))
+}
+
 // Minimize adjust this vertex v to hold the minimum
 // values component-wise of v and w.
 func (v *Vertex) Minimize(w Vertex) {
@@ -441,6 +446,34 @@
 	}
 }
 
+func (ls LineStringZ) Buffer(radius float64) {
+	if len(ls) == 0 {
+		return
+	}
+	var cx, cy, cz float64
+	for i := range ls {
+		cx += ls[i].X
+		cy += ls[i].Y
+		cz += ls[i].Z
+	}
+
+	s := 1.0 / float64(len(ls))
+
+	cx *= s
+	cy *= s
+	cz *= s
+	c := Vertex{X: cx, Y: cy, Z: cz}
+	for i := range ls {
+		n := c.Sub(ls[i]).Normalize()
+		l := ls[i].Distance(c)
+		// l = 1
+		// l + r = s
+		// (l+r)/l = s
+		s := (l + radius) / l
+		ls[i] = ls[i].Add(n.Scale(s))
+	}
+}
+
 func (ls LineStringZ) order(position func(Vertex) float64) {
 	type posVertex struct {
 		pos float64