diff pkg/imports/sr.go @ 3655:a6c671abbc35 single-beam

Started with cleaning up the single beam import.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Thu, 13 Jun 2019 18:33:54 +0200
parents 26e02fa5b992
children 1c3df921361d
line wrap: on
line diff
--- a/pkg/imports/sr.go	Thu Jun 13 17:07:28 2019 +0200
+++ b/pkg/imports/sr.go	Thu Jun 13 18:33:54 2019 +0200
@@ -24,15 +24,12 @@
 	"errors"
 	"fmt"
 	"io"
-	"log"
 	"math"
 	"os"
 	"path"
 	"path/filepath"
-	"runtime"
 	"strconv"
 	"strings"
-	"sync"
 	"time"
 
 	shp "github.com/jonas-p/go-shp"
@@ -151,14 +148,6 @@
   ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.0)),
   ST_AsBinary(ST_Buffer(ST_MakeValid(ST_GeomFromWKB($1, $2::integer)), 0.1))`
 
-	whatsWrongSQL = `
-SELECT
-  ST_IsValidReason(
-    CAST(ST_MakeValid(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326))::geography AS geometry)),
-  ST_IsValid(
-    CAST(ST_MakeValid(ST_Transform(ST_GeomFromWKB($1, $2::integer), 4326))::geography AS geometry))
-`
-
 	insertContourSQL = `
 INSERT INTO waterway.sounding_results_contour_lines (
   sounding_result_id,
@@ -384,7 +373,8 @@
 	var clippingPolygon octree.Polygon
 
 	if boundary == nil {
-		feedback.Info("No boundary given.")
+		feedback.Info("No boundary given. Calulating from XYZ data.")
+		feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge)
 
 		polygon, removed := tri.ConcaveHull(tooLongEdge)
 
@@ -396,19 +386,6 @@
 
 		clippingPolygon.FromLineStringZ(polygon)
 
-		var wrong string
-		var valid bool
-		if err := tx.QueryRowContext(
-			ctx,
-			whatsWrongSQL,
-			clippingPolygon.AsWKB(),
-			epsg,
-		).Scan(&wrong, &valid); err != nil {
-			return nil, err
-		}
-
-		log.Printf("!!!!!!!!: whats wrong (%t): %s\n", valid, wrong)
-
 		var repaired, buffered []byte
 
 		if err := tx.QueryRowContext(
@@ -434,115 +411,58 @@
 
 		firstTree := builder.Tree()
 
-		bbox := polygon.BBox()
-		bboxArea := bbox.Area()
-
-		log.Printf("bbox area: %.2f\n", bboxArea)
-		log.Printf("polygon area: %.2f\n", polygonArea)
+		feedback.Info("Boundary area: %.2fm²", polygonArea)
 
 		numPoints := int(math.Ceil(polygonArea * pointsPerSquareMeter))
 
-		log.Printf("random points for density of ~%d points/m²: %d\n",
-			pointsPerSquareMeter, numPoints)
+		feedback.Info("Generating %d random points for an average density of ~%d points/m².",
+			numPoints, pointsPerSquareMeter)
 
 		start := time.Now()
 
-		var wg sync.WaitGroup
-
-		jobs := make(chan int)
-		out := make(chan []octree.Vertex)
-		done := make(chan struct{})
-
-		cpus := runtime.NumCPU()
-
-		free := make(chan []octree.Vertex, cpus)
-
-		for i, n := 0, cpus; i < n; i++ {
-			wg.Add(1)
-			go func() {
-				defer wg.Done()
-				xRange := common.Random(bbox.X1, bbox.X2)
-				yRange := common.Random(bbox.Y1, bbox.Y2)
-
-				for size := range jobs {
-					var vertices []octree.Vertex
-					select {
-					case vertices = <-free:
-					default:
-						vertices = make([]octree.Vertex, 0, 1000)
-					}
-					for len(vertices) < size {
-						x, y := xRange(), yRange()
-						if z, ok := firstTree.Value(x, y); ok {
-							vertices = append(vertices, octree.Vertex{X: x, Y: y, Z: z})
-						}
-					}
-					out <- vertices
-				}
-			}()
-		}
-
 		generated := make(octree.LineStringZ, 0, numPoints+len(polygon)-1)
 
-		go func() {
-			defer close(done)
-			for vertices := range out {
-				generated = append(generated, vertices...)
-				select {
-				case free <- vertices[:0]:
-				default:
-				}
-			}
-		}()
+		firstTree.GenerateRandomVertices(numPoints, func(vertices []octree.Vertex) {
+			generated = append(generated, vertices...)
+		})
 
-		for remain := numPoints; remain > 0; {
-			if remain > 1000 {
-				jobs <- 1000
-				remain -= 1000
-			} else {
-				jobs <- remain
-				remain = 0
-			}
-		}
-		close(jobs)
-		wg.Wait()
-		close(out)
-		<-done
-
-		log.Printf("Generating %d points took %v\n", len(generated), time.Since(start))
+		feedback.Info("Generating %d points took %v.", len(generated), time.Since(start))
 
 		// Add the boundary to new point cloud.
 		generated = append(generated, polygon[:len(polygon)-1]...)
 
+		feedback.Info("Triangulate new point cloud.")
+
+		start = time.Now()
+
 		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)
+		feedback.Info("Second triangulation took %v.", time.Since(start))
+		feedback.Info("Number triangles: %d.", len(tri.Triangles)/3)
+		feedback.Info("Clipping triangles from new mesh.")
 
+		start = time.Now()
 		secondTin := tri.Tin()
-
-		//polygon.Buffer(1.0)
-
+		secondTin.EPSG = uint32(epsg)
 		clippingPolygonBuffered.Indexify()
 
 		var str octree.STRTree
 		str.Build(secondTin)
-		log.Printf("Building STR tree took %v", time.Since(start))
+		feedback.Info("Building STR tree took %v", time.Since(start))
 
 		start = time.Now()
 
 		removed = str.Clip(&clippingPolygonBuffered)
-		log.Printf("Clipping STR tree took %v.", time.Since(start))
-		log.Printf("Number of triangles to clip %d.", len(removed))
+		feedback.Info("Clipping STR tree took %v.", time.Since(start))
+		feedback.Info("Number of triangles to clip %d.", len(removed))
 
 		start = time.Now()
 
-		secondTin.EPSG = uint32(epsg)
+		feedback.Info("Building final octree index")
 
 		builder = octree.NewBuilder(secondTin)
 		builder.Build(removed)
@@ -552,6 +472,8 @@
 		}
 		feedback.Info("Building octree took %v.", time.Since(start))
 
+		start = time.Now()
+
 		var (
 			id       int64
 			dummy    uint32
@@ -579,7 +501,6 @@
 			return nil, err
 		}
 
-		start = time.Now()
 		h := sha1.New()
 		h.Write(octreeIndex)
 		checksum := hex.EncodeToString(h.Sum(nil))