# HG changeset patch # User Sascha L. Teichmann # Date 1560424480 -7200 # Node ID 01ce3ba9b0d07ca62ea9aaf903c1afdc5a5a0827 # Parent 810b28f59b8baacd61ea8f73bcfea0bcc7d0de8b Fixed generating of random points. diff -r 810b28f59b8b -r 01ce3ba9b0d0 pkg/common/random.go --- 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 } } diff -r 810b28f59b8b -r 01ce3ba9b0d0 pkg/imports/sr.go --- 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( diff -r 810b28f59b8b -r 01ce3ba9b0d0 pkg/octree/polygon.go --- 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} +} diff -r 810b28f59b8b -r 01ce3ba9b0d0 pkg/octree/vertex.go --- 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