Mercurial > gemma

comparison pkg/imports/sr.go @ 4654:3eda5a7215ab stree-experiment

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Generate STRTrees instaead of Octree during sounding result imports.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Mon, 14 Oct 2019 13:12:38 +0200
parents 89a72e0e2f9b
children 6eab3ac0e849
comparison
equal deleted inserted replaced
4653:8efc6b3289f3 4654:3eda5a7215ab
22 "database/sql" 22 "database/sql"
23 "encoding/hex" 23 "encoding/hex"
24 "errors" 24 "errors"
25 "fmt" 25 "fmt"
26 "io" 26 "io"
27 "log"
28 "math" 27 "math"
29 "os" 28 "os"
30 "path" 29 "path"
31 "path/filepath" 30 "path/filepath"
32 "strconv" 31 "strconv"
154 ST_AsBinary(ST_Buffer(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer), 0.1)), 153 ST_AsBinary(ST_Buffer(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer), 0.1)),
155 ST_Area(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))` 154 ST_Area(ST_Transform(ST_GeomFromWKB($1, $2::integer), $3::integer))`
156 155
157 insertOctreeSQL = ` 156 insertOctreeSQL = `
158 UPDATE waterway.sounding_results SET 157 UPDATE waterway.sounding_results SET
159 octree_checksum = $2, octree_index = $3 158 mesh_checksum = $2, mesh_index = $3
160 WHERE id = $1` 159 WHERE id = $1`
161 160
162 repairBoundarySQL = ` 161 repairBoundarySQL = `
163 SELECT 162 SELECT
164 ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.0)), 163 ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.0)),
496 495
497 feedback.Info("Build virtual DEM based on original XYZ data.") 496 feedback.Info("Build virtual DEM based on original XYZ data.")
498 497
499 start = time.Now() 498 start = time.Now()
500 499
501 var tree *octree.Tree 500 tin := tri.Tin()
502 { 501 var virtual octree.STRTree
503 builder := octree.NewBuilder(tri.Tin()) 502 virtual.BuildWithout(tin, removed)
504 builder.Build(removed)
505 tree = builder.Tree()
506 }
507 503
508 feedback.Info("Building took %v", time.Since(start)) 504 feedback.Info("Building took %v", time.Since(start))
509 505
510 feedback.Info("Boundary area: %.2fm²", polygonArea) 506 feedback.Info("Boundary area: %.2fm²", polygonArea)
511 507
516 512
517 start = time.Now() 513 start = time.Now()
518 514
519 generated := make(octree.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0)) 515 generated := make(octree.LineStringZ, 0, numPoints+clippingPolygon.NumVertices(0))
520 516
521 tree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) { 517 octree.GenerateRandomVertices(
522 generated = append(generated, vertices...) 518 numPoints,
523 }) 519 tin.Min, tin.Max,
520 virtual.Value,
521 func(vertices []octree.Vertex) {
522 generated = append(generated, vertices...)
523 })
524 524
525 feedback.Info("Generating %d points took %v.", len(generated), time.Since(start)) 525 feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
526 526
527 // Add the boundary to new point cloud. 527 // Add the boundary to new point cloud.
528 dupes := map[[2]float64]struct{}{} 528 dupes := map[[2]float64]struct{}{}
530 key := [2]float64{x, y} 530 key := [2]float64{x, y}
531 if _, found := dupes[key]; found { 531 if _, found := dupes[key]; found {
532 return 532 return
533 } 533 }
534 dupes[key] = struct{}{} 534 dupes[key] = struct{}{}
535 if z, ok := tree.Value(x, y); ok { 535 if z, ok := virtual.Value(x, y); ok {
536 generated = append(generated, octree.Vertex{X: x, Y: y, Z: z}) 536 generated = append(generated, octree.Vertex{X: x, Y: y, Z: z})
537 } 537 }
538 }) 538 })
539 539
540 feedback.Info("Triangulate new point cloud.") 540 feedback.Info("Triangulate new point cloud.")
546 return nil, err 546 return nil, err
547 } 547 }
548 feedback.Info("Second triangulation took %v.", time.Since(start)) 548 feedback.Info("Second triangulation took %v.", time.Since(start))
549 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3) 549 feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
550 feedback.Info("Clipping triangles from new mesh.") 550 feedback.Info("Clipping triangles from new mesh.")
551
552 } 551 }
553 552
554 start = time.Now() 553 start = time.Now()
555 tin = tri.Tin() 554 tin = tri.Tin()
556 tin.EPSG = epsg 555 tin.EPSG = epsg
565 removed = str.Clip(&clippingPolygonBuffered) 564 removed = str.Clip(&clippingPolygonBuffered)
566 feedback.Info("Clipping STR tree took %v.", time.Since(start)) 565 feedback.Info("Clipping STR tree took %v.", time.Since(start))
567 feedback.Info("Number of triangles to clip %d.", len(removed)) 566 feedback.Info("Number of triangles to clip %d.", len(removed))
568 567
569 start = time.Now() 568 start = time.Now()
570 s := octree.STRTree{Entries: 16} 569 final := octree.STRTree{Entries: 16}
571 s.BuildWithout(tin, removed) 570 final.BuildWithout(tin, removed)
572 if _, err2 := s.Bytes(); err2 != nil { 571
573 log.Printf("serializing STRTree failed: %v\n", err2)
574 }
575 log.Printf("Building strtree took: %v.\n", time.Since(start))
576
577 // return nil, UnchangedError("nothing to do")
578
579 feedback.Info("Build final octree index")
580
581 builder := octree.NewBuilder(tin)
582 builder.Build(removed)
583 octreeIndex, err := builder.Bytes()
584 if err != nil {
585 return nil, err
586 }
587 feedback.Info("Building octree took %v.", time.Since(start)) 572 feedback.Info("Building octree took %v.", time.Since(start))
588 feedback.Info("Store octree.") 573 feedback.Info("Store octree.")
589 574
590 start = time.Now() 575 start = time.Now()
591 576
621 "no matching bottleneck of given name or time available: %v", err) 606 "no matching bottleneck of given name or time available: %v", err)
622 case err != nil: 607 case err != nil:
623 return nil, err 608 return nil, err
624 } 609 }
625 610
611 index, err := final.Bytes()
612 if err != nil {
613 return nil, err
614 }
615
626 h := sha1.New() 616 h := sha1.New()
627 h.Write(octreeIndex) 617 h.Write(index)
628 checksum := hex.EncodeToString(h.Sum(nil)) 618 checksum := hex.EncodeToString(h.Sum(nil))
629 _, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, octreeIndex) 619 _, err = tx.ExecContext(ctx, insertOctreeSQL, id, checksum, index)
630 if err != nil { 620 if err != nil {
631 return nil, err 621 return nil, err
632 } 622 }
633 feedback.Info("Storing octree index took %s.", time.Since(start)) 623 feedback.Info("Storing octree index took %s.", time.Since(start))
634 err = generateIsos(ctx, tx, feedback, builder.Tree(), id) 624 err = generateIsos(ctx, tx, feedback, &final, id)
635 if err != nil { 625 if err != nil {
636 return nil, err 626 return nil, err
637 } 627 }
638 628
639 // Store for potential later removal. 629 // Store for potential later removal.
836 826
837 func generateIsos( 827 func generateIsos(
838 ctx context.Context, 828 ctx context.Context,
839 tx *sql.Tx, 829 tx *sql.Tx,
840 feedback Feedback, 830 feedback Feedback,
841 tree *octree.Tree, 831 tree *octree.STRTree,
842 id int64, 832 id int64,
843 ) error { 833 ) error {
844 834
845 heights, err := octree.LoadClassBreaks( 835 heights, err := octree.LoadClassBreaks(
846 ctx, tx, 836 ctx, tx,
847 "morphology_classbreaks", 837 "morphology_classbreaks",
848 ) 838 )
849 839
840 minZ, maxZ := tree.Min().Z, tree.Max().Z
841
850 if err != nil { 842 if err != nil {
851 feedback.Warn("Loading class breaks failed: %v", err) 843 feedback.Warn("Loading class breaks failed: %v", err)
852 feedback.Info("Using default class breaks") 844 feedback.Info("Using default class breaks")
853 heights = nil 845 heights = nil
854 h := contourStepWidth * math.Ceil(tree.Min.Z/contourStepWidth) 846 h := contourStepWidth * math.Ceil(minZ/contourStepWidth)
855 for ; h <= tree.Max.Z; h += contourStepWidth { 847 for ; h <= maxZ; h += contourStepWidth {
856 heights = append(heights, h) 848 heights = append(heights, h)
857 } 849 }
858 } else { 850 } else {
859 heights = octree.ExtrapolateClassBreaks(heights, tree.Min.Z, tree.Max.Z) 851 heights = octree.ExtrapolateClassBreaks(heights, minZ, maxZ)
860 // We set steps for InBetweenClassBreaks to 1, so it 852 // We set steps for InBetweenClassBreaks to 1, so it
861 // becomes a null operation. The extra class breaks 853 // becomes a null operation. The extra class breaks
862 // were considered unexpected and confusing by the 854 // were considered unexpected and confusing by the
863 // users. Once we get filled polygones the visual will 855 // users. Once we get filled polygones the visual will
864 // be considerably different anyway. -- sw 856 // be considerably different anyway. -- sw
879 871
880 func generateIsoAreas( 872 func generateIsoAreas(
881 ctx context.Context, 873 ctx context.Context,
882 tx *sql.Tx, 874 tx *sql.Tx,
883 feedback Feedback, 875 feedback Feedback,
884 tree *octree.Tree, 876 tree *octree.STRTree,
885 heights []float64, 877 heights []float64,
886 id int64, 878 id int64,
887 ) error { 879 ) error {
888 feedback.Info("Generate iso areas") 880 feedback.Info("Generate iso areas")
889 total := time.Now() 881 total := time.Now()
890 defer func() { 882 defer func() {
891 feedback.Info("Generating iso areas took %s.", 883 feedback.Info("Generating iso areas took %s.",
892 time.Since(total)) 884 time.Since(total))
893 }() 885 }()
894 886
895 areas := octree.TraceAreas(heights, isoCellSize, tree.Min, tree.Max, tree.Value) 887 areas := octree.TraceAreas(heights, isoCellSize, tree.Min(), tree.Max(), tree.Value)
896 888
897 return storeAreas( 889 return storeAreas(
898 ctx, tx, feedback, 890 ctx, tx, feedback,
899 areas, tree.EPSG, heights, id) 891 areas, tree.EPSG(), heights, id)
900 } 892 }
901 893
902 func storeAreas( 894 func storeAreas(
903 ctx context.Context, 895 ctx context.Context,
904 tx *sql.Tx, 896 tx *sql.Tx,