Mercurial > gemma
view cmd/octreediff/main.go @ 2576:647a58ee9ae9
Morphological differences: Centralized generation of height values for differences in octree package.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Mon, 11 Mar 2019 15:04:23 +0100 |
parents | 59e7a011d347 |
children | 5c6a7e69b02d |
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 ) ` 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) } } 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 result := first.Diff(second) now = time.Now() log.Printf("setting in took %v\n", now.Sub(last)) last = now log.Printf("num points: %d\n", len(result)) clippingPolygon, err := octree.LoadClippingPolygon( ctx, conn, first.EPSG, bottleneck, firstDate, secondDate, ) if err != nil { return err } 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 tin := tri.Tin() removed := tin.Clip(clippingPolygon) now = time.Now() log.Printf("clipping tin 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) heights := octree.SampleDiffHeights(tree.Min.Z, tree.Max.Z, 0.1) 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)) }