Mercurial > gemma
view cmd/octreediff/main.go @ 2562:ce39e9954e85
Make upload of AGM require only "fk_gauge_id" "measure_date" and "value"
author | Sascha Wilde <wilde@intevation.de> |
---|---|
date | Fri, 08 Mar 2019 18:57:58 +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)) }