Mercurial > gemma
view cmd/octreediff/main.go @ 2526:6498267096ae
client: critical bottlenecks: use png instead of vectors for marker
The positioning of text is not always precise. The exclamation mark in the warning sign for critical bottlenecks
was not always perfectly centered. Using a png image is more reliable.
author | Markus Kottlaender <markus@intevation.de> |
---|---|
date | Wed, 06 Mar 2019 15:39:59 +0100 |
parents | 1ec4c5633eb6 |
children | 7686c7c23506 |
line wrap: on
line source
// This is Free Software under GNU Affero General Public License v >= 3.0 // without warranty, see README.md and license for details. // // SPDX-License-Identifier: AGPL-3.0-or-later // License-Filename: LICENSES/AGPL-3.0.txt // // Copyright (C) 2018 by via donau // – Österreichische Wasserstraßen-Gesellschaft mbH // Software engineering by Intevation GmbH // // Author(s): // * Sascha L. Teichmann <sascha.teichmann@intevation.de> package main import ( "context" "database/sql" "errors" "flag" "fmt" "log" "runtime" "sync" "time" "gemma.intevation.de/gemma/pkg/common" "gemma.intevation.de/gemma/pkg/octree" ) var ( bottleneck = flag.String("bottleneck", "", "name of the bottleneck") first = flag.String("first", "", "date of the first sounding result") second = flag.String("second", "", "date of the second sounding result") ) const contourTolerance = 0.1 const ( loadOctreeSQL = ` SELECT sr.octree_index FROM waterway.sounding_results sr JOIN waterway.bottlenecks bn ON sr.bottleneck_id = bn.id WHERE bn.bottleneck_id = $1 AND sr.date_info = $2::date AND sr.octree_index IS NOT NULL` insertContourSQL = ` INSERT INTO diff_contour_lines ( height, lines ) SELECT $4, ST_Transform( ST_Multi( ST_CollectionExtract( ST_SimplifyPreserveTopology( ST_Multi(ST_Collectionextract( ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 2)), $3 ), 2 ) ), 4326 ) ` clippingPolygonSQL = ` WITH joined AS ( SELECT sr.area AS area, sr.date_info AS date_info FROM waterway.sounding_results sr JOIN waterway.bottlenecks bn ON sr.bottleneck_id = bn.id WHERE bn.bottleneck_id = $1 ) SELECT ST_AsBinary( ST_intersection( (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date), (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date) ) ) AS area ` insertContourSQLClipped = ` WITH joined AS ( SELECT sr.area AS area, sr.date_info AS date_info FROM waterway.sounding_results sr JOIN waterway.bottlenecks bn ON sr.bottleneck_id = bn.id WHERE bn.bottleneck_id = $1 ), inter AS ( SELECT ST_Buffer( ST_intersection( (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date), (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date) ), 0.001) AS area ) INSERT INTO diff_contour_lines ( height, lines ) SELECT $5, ST_Transform( ST_Multi( ST_intersection( ST_CollectionExtract( ST_SimplifyPreserveTopology( ST_Multi(ST_Collectionextract( ST_MakeValid(ST_GeomFromWKB($6, $2::integer)), 2)), $7 ), 2 ), area ) ), 4326 ) FROM inter` ) func check(err error) { if err != nil { log.Fatalf("error: %v\n", err) } } type point struct { x float64 y float64 } type pointMap map[point]float64 func (pm pointMap) triangulate() (*octree.Triangulation, error) { vertices := make([]octree.Vertex, len(pm)) var i int for p, z := range pm { vertices[i] = octree.Vertex{X: p.x, Y: p.y, Z: z} i++ } return octree.Triangulate(vertices) } func sliceWork( vs []octree.Vertex, dst pointMap, fn func([]octree.Vertex, func([]octree.Vertex) []octree.Vertex), ) { n := runtime.NumCPU() wg := new(sync.WaitGroup) slices := make(chan []octree.Vertex) out := make(chan []octree.Vertex) pool := make(chan []octree.Vertex, n) const pageSize = 2048 turn := func(p []octree.Vertex) []octree.Vertex { if p != nil { out <- p } select { case p = <-pool: default: p = make([]octree.Vertex, 0, pageSize) } return p } for i := 0; i < n; i++ { wg.Add(1) go func() { defer wg.Done() for slice := range slices { fn(slice, turn) } }() } done := make(chan struct{}) go func() { defer close(done) for s := range out { for i := range s { v := &s[i] key := point{x: v.X, y: v.Y} if z, found := dst[key]; found { dst[key] = (z + v.Z) * 0.5 } else { dst[key] = v.Z } } select { case pool <- s[:0:pageSize]: default: } } }() size := len(vs)/n + 1 for len(vs) > 0 { var l int if len(vs) < size { l = len(vs) } else { l = size } slices <- vs[:l] vs = vs[l:] } close(slices) wg.Wait() close(out) <-done } func process(bottleneck string, firstDate, secondDate time.Time) error { start := time.Now() defer func() { log.Printf("processing took %v\n", time.Since(start)) }() ctx := context.Background() return run(func(db *sql.DB) error { conn, err := db.Conn(ctx) if err != nil { return err } defer conn.Close() tx, err := conn.BeginTx(ctx, nil) if err != nil { return err } defer tx.Rollback() type load struct { date time.Time data []byte err *error dst **octree.Tree } out := make(chan *load) wg := new(sync.WaitGroup) n := runtime.NumCPU() if n > 2 { n = 2 } for i := 0; i < n; i++ { wg.Add(1) go func() { defer wg.Done() for l := range out { if *l.err == nil { *l.dst, *l.err = octree.Deserialize(l.data) } } }() } var firstErr, secondErr error var first, second *octree.Tree for _, l := range []*load{ {date: firstDate, dst: &first, err: &firstErr}, {date: secondDate, dst: &second, err: &secondErr}, } { var data []byte if err := tx.QueryRowContext( ctx, loadOctreeSQL, bottleneck, l.date, ).Scan(&data); err != nil { *l.err = err } else { l.data = data } out <- l } close(out) wg.Wait() if firstErr != nil || secondErr != nil { if firstErr != nil && secondErr != nil { return fmt.Errorf("%v, %v", firstErr, secondErr) } if firstErr != nil { return firstErr } return secondErr } if first.EPSG != second.EPSG { return errors.New("EPSG codes mismatch. Needs transformation slow pass.") } now := time.Now() log.Printf("loading took %v\n", now.Sub(start)) last := now firstVs, secondVs := first.Vertices(), second.Vertices() result := make(pointMap, len(firstVs)+len(secondVs)) sliceWork( firstVs, result, func( slice []octree.Vertex, turn func([]octree.Vertex) []octree.Vertex, ) { p := turn(nil) for i := range slice { v := &slice[i] if z, found := second.Value(v.X, v.Y); found { p = append(p, octree.Vertex{v.X, v.Y, v.Z - z}) if len(p) == cap(p) { p = turn(p) } } } if len(p) > 0 { turn(p) } }) sliceWork( secondVs, result, func( slice []octree.Vertex, turn func([]octree.Vertex) []octree.Vertex, ) { p := turn(nil) for i := range slice { v := &slice[i] if z, found := first.Value(v.X, v.Y); found { p = append(p, octree.Vertex{v.X, v.Y, z - v.Z}) if len(p) == cap(p) { p = turn(p) } } } if len(p) > 0 { turn(p) } }) now = time.Now() log.Printf("setting in took %v\n", now.Sub(last)) last = now log.Printf("num points: %d\n", len(result)) var clip []byte if err := tx.QueryRowContext( ctx, clippingPolygonSQL, bottleneck, first.EPSG, firstDate, secondDate, ).Scan(&clip); err != nil { return err } var clippingPolygon octree.Polygon if err := clippingPolygon.FromWKB(clip); err != nil { return err } clippingPolygon.Indexify() now = time.Now() log.Printf("loading clipping polygon took %v\n", now.Sub(last)) last = now tri, err := result.triangulate() if err != nil { return err } now = time.Now() log.Printf("triangulation took %v\n", now.Sub(last)) last = now var str octree.STRTree tin := tri.Tin() str.Build(tin) now = time.Now() log.Printf("building STR tree took %v\n", now.Sub(last)) last = now removed := str.Clip(&clippingPolygon) now = time.Now() log.Printf("clipping STR tree took %v\n", now.Sub(last)) last = now log.Printf("Number of triangles to clip: %d\n", len(removed)) builder := octree.NewBuilder(tin) builder.Build(removed) now = time.Now() log.Printf("building octree took %v\n", now.Sub(last)) last = now tree := builder.Tree() log.Printf("min/max: %f %f\n", tree.Min.Z, tree.Max.Z) var heights []float64 switch { case tree.Min.Z >= 0: // All values positive. for v := 0.0; v <= tree.Max.Z; v += 0.1 { if v >= tree.Min.Z { heights = append(heights, v) } } case tree.Max.Z <= 0: // All values negative. for v := 0.0; v >= tree.Min.Z; v -= 0.1 { if v <= tree.Max.Z { heights = append(heights, v) } } default: // Positive and negative. for v := 0.1; v <= tree.Max.Z; v += 0.1 { heights = append(heights, v) } for i, j := 0, len(heights)-1; i < j; i, j = i+1, j-1 { heights[i], heights[j] = heights[j], heights[i] } for v := 0.0; v >= tree.Min.Z; v -= 0.1 { heights = append(heights, v) } } log.Printf("num heights: %d\n", len(heights)) var dataSize int stmt, err := tx.PrepareContext(ctx, insertContourSQL) if err != nil { return err } defer stmt.Close() log.Println("do contours") octree.DoContours(tree, heights, func(res *octree.ContourResult) { if err == nil && len(res.Lines) > 0 { log.Printf("%.2f: lines: %d\n", res.Height, len(res.Lines)) wkb := res.Lines.AsWKB2D() dataSize += len(wkb) _, err = stmt.ExecContext( ctx, wkb, first.EPSG, contourTolerance, res.Height, ) } }) if err != nil { return err } now = time.Now() log.Printf("Number of iso lines: %d\n", len(heights)) log.Printf("Total WKB size: %.2fMB\n", float64(dataSize)/(1024*1024)) log.Printf("generating iso lines took %v\n", now.Sub(last)) last = now return tx.Commit() }) } func main() { flag.Parse() firstDate, err := time.Parse(common.DateFormat, *first) check(err) secondDate, err := time.Parse(common.DateFormat, *second) check(err) if *bottleneck == "" { log.Fatalln("Missing bottleneck name") } check(process(*bottleneck, firstDate, secondDate)) }